A high throughput computational investigation of the solid solution mechanisms of actinides and lanthanides in zirconolite

Abstract

In this work, we perform a theoretical investigation of the actinide and lanthanide solid solution mechanisms of zirconolite-2M, prototypically CaZrTi₂O₇, as a candidate immobilisation matrix for plutonium. Solid solution energies were calculated using static atomistic simulations by means of the General Utility Lattice Program, for formulations of relevance to ceramic wasteform deployment, with substitution on the Ca²⁺ and Zr⁴⁺ sites by Ce⁴⁺, Pu⁴⁺, Th⁴⁺, and U⁴⁺, and appropriate charge balance by substitution of Al³⁺ or Fe³⁺ on Ti⁴⁺ sites. In simple solid solutions involving substitution on the Zr⁴⁺ site, we found that whereas substitution of Ce⁴⁺, U⁴⁺ and Pu⁴⁺ were energetically favoured, substitution of Th⁴⁺ was not energetically favoured. For more complex solid solutions involving Ce⁴⁺, Pu⁴⁺, Th⁴⁺, and U⁴⁺ substitution on the Ca²⁺ site, we found the most energetically favoured scheme involved co-substitution of Al³⁺ or Fe³⁺ on the five-fold co-ordinate Ti⁴⁺ site in the zirconolite-2M structure. Comparison of these computational data with experimental evidence, where available, demonstrated broad agreement. Consequently, this study provides useful insight into formulation design and the efficacy of Ce⁴⁺, U⁴⁺ and Th⁴⁺ as Pu⁴⁺ surrogates in zirconolite-2M ceramic wasteforms for plutonium disposition.

---

1 Introduction

The UK holds the world's largest inventory of separated civil plutonium, forecast to reach 140 tons at the end of reprocessing operations.¹ UK government policy is to manage this material to a safe and secure end point, the preferred strategy for which is reuse in mixed oxide (MOX) fuel in light water reactors. However, should this strategy not prove implementable, immobilisation of the inventory will be required, along with the fraction of material known to be unsuitable for MOX fuel manufacture.

Numerous natural and synthetic materials have been proposed as wasteforms for the immobilisation of actinides, these including ceramics, glasses, and glass-ceramics.^2–6 Geological disposal of actinides places greater emphasis on the performance of the wasteform and near field barriers, so as to assure adequate containment of fissile material over the required timescales, which, in the geological context, are comparably short.^2–5,7 Zirconolite, prototypically CaZrTi₂O₇, is a naturally occurring mineral and the dominant actinide bearing phase in the SYNROC C ceramic wasteform;^8,9 it is known to be highly resistant to alteration and dissolution.^10,11 As a result, zirconolite is an attractive material for plutonium immobilisation and its potential as a wasteform has been well established.^8,12–15 The zirconolite-2M polytype structure (space group C2/c), adopted by the prototypical composition CaZrTi₂O₇, comprises alternating layers of CaO₈ and ZrO₇ polyhedra aligned parallel to (110); parallel to the [001] direction, these polyhedra are interleaved 1 : 1 with hexagonal tungsten bronze motifs formed by corner sharing TiO₆ and TiO₅ polyhedra.^16,17 The 2M nomenclature thus signifies a monoclinic unit cell with a two layer repeat sequence along [001]; other zirconolite polytype structures with different interlayer relationships are known, as discussed below.

The use of Pu in laboratory based studies is hazardous, challenging and expensive. Consequently, Ce, U and Th are frequently used as inactive or low active surrogates to emulate the behaviour of Pu in laboratory based studies.^18,19 This is due to the similarity of the ionic size of Ce⁴⁺, U⁴⁺, Th⁴⁺ and Pu⁴⁺ and to CeO₂, PuO₂, ThO₂, and UO₂ all having a common fluorite crystal structure and exhibiting solid solution at any ratio, implying similar solid state chemistry.^20–23

In this work we aim to investigate the plutonium immobilisation potential of zirconolite-2M by atomistic simulations. Previous simulation based studies of zirconolite-2M at the atomistic level have focused on studying the defect chemistry of actinide additions,²⁴ and within the regime of molecular dynamics for investigations the radiation damaged structure, and crystalline to amorphous phase transition, arising from α-decay of Pu.^25–27 More recently, Ce and actinide solid solution mechanisms in zirconolite-2M were studied at the electronic structure level, within the density functional theory (DFT) regime.²⁸ Importantly, the DFT investigation of Tanti et al. broadly agreed with findings of the atomistic simulations of Gilbert et al. for Ce-substituted zirconolite-2M. The choice of methodology and accuracy level is a critical consideration in such investigations. The broad agreement between DFT and atomistic simulations shows that we can obtain accurate insight with the atomistic approach. Further, given the low computational cost, it is feasible to simulate relatively large lattices at the atomistic level employing an innovative high-throughput workflow described herein.

Our investigation develops and extends a previous computational study of the defect chemistry of zirconolite-2M,²⁴ with regard to incorporation of Ce^3+/4+ and Pu^3+/4+. In this contribution we expand the previous study by examining the solid solution of Pu⁴⁺ and its typical surrogates, Ce⁴⁺, U⁴⁺, Th⁴⁺, on the Ca²⁺ and Zr⁴⁺ sites, at concentrations greater than point defects, with necessary charge compensation provided by Al³⁺ and Fe³⁺ substitution on the Ti⁴⁺ sites.

We focus on the Pu⁴⁺ oxidation state which has been shown to be the dominant species in fluorite related zirconolite-2M and pyrochlore structured ceramics synthesised under conditions relevant to wasteform manufacture.^29–33 Under conditions of hot isostatic pressing with PuO₂ as a feedstock, synthesis of the zirconolite ceramic wasteform will be under the redox control of the Fe/FeO buffer imposed by the stainless steel can. Consideration of Ellingham diagrams shows that this will not be sufficient to effect PuO₂/Pu₂O₃ reduction.³⁴ Indeed, Pu³⁺ is stabilised by annealing only under strongly reducing 5% H₂/N₂, or 5% H₂/Ar, which is not relevant to the technological focus of wasteform manufacture by hot isostatic pressing.^31,32

2 Theory

The work presented in this paper examines three substitution schemes to investigate Pu and surrogate incorporation in zirconolite-2M. The substitution schemes are based on compositions relevant to optimisation of zirconolite ceramic formulations. We used the supercell approach to study the defects of the system. Here, the defects are added as absolute concentrations in a solid solution, and our concentration values can therefore be directly compared to experimental compositions. This differs from previous work on this system²⁴ that used a Mott–Littleton method where the substitution defects were effectively at infinite dilution.

In the first substitution scheme, we replaced Zr⁴⁺ sites in prototypical zirconolite-2M with Ce⁴⁺, Pu⁴⁺, Th⁴⁺ and U⁴⁺. The chemical reaction for the substitution scheme was as follows,

CaZrTi₂O₇ + xMO₂ → CaZr_1−xM_xTi₂O₇ + xZrO₂

where M = Ce, Pu, Th, U. Here, 2M denotes the polytype structure of monoclinic symmetry.

The second substitution scheme targeted the substitution of Ce⁴⁺, Pu⁴⁺, Th⁴⁺ and U⁴⁺ on the Ca²⁺ site, with charge balance provided by replacing 2Ti⁴⁺ sites with Al³⁺, for every Ca²⁺ ion replaced.^35–37 The second substitution scheme followed the reaction,

CaZrTi₂O₇ + xMO₂ + xAl₂O₃ → Ca_1−xM_xZrTi_2−2xAl_2xO₇ + xCaO + 2xTiO₂

where M = Ce, Pu, Th, U.

There are three unique Ti⁴⁺ sites in the zirconolite-2M structure that may accommodate charge balancing cations such as Al³⁺. An illustration of the Ti⁴⁺ site orientations and coordination is shown in Fig. 1, where example Ti(1), Ti(2) and Ti(3) sites are coloured in green, yellow and fuschia respectively, with Ca and Zr omitted for clarity. The Ti(1) and Ti(3) sites adopt octahedral co-ordination by O²⁻, whereas the Ti(2) site adopts a trigonal bipyramidal configuration by O²⁻. The Ti(2) site is partially occupied, with a 50% probability of lying either side of the site axis.¹⁶

Fig. 1 Ti⁴⁺ site orientation in our base zirconolite-2M system (Ca²⁺, Zr⁴⁺ hidden). Ti sites are labelled, example Ti(1), Ti(2) and Ti(3) sites are coloured green, yellow and fuchsia, respectively. The Ti(1) and Ti(3) sites have octahedral coordination, whereas the Ti(2) site has trigonal bipyramidal co-ordination and is 50% occupied. Figure generated in VESTA.⁴⁰

Experimental studies have shown that charge balancing ions are generally preferentially accommodated in the 5-fold coordinate Ti(2) site.^36–38 Although, in some instances, charge balancing species such as Cr³⁺ have been shown to preferentially adopt 6-fold Ti⁴⁺ sites as may be expected from consideration of crystal field stabilisation energy.³⁴ Therefore, this substitution scheme needs to consider potential preferential substitution of the charge balancing ions for particular Ti⁴⁺ sites. To address this question, we considered 6 different Al³⁺ site combinations: two Ti(1) sites; two Ti(2) sites; two Ti(3) sites; one Ti(1) and one Ti(2); one Ti(1) and one Ti(3); and one Ti(2) and one Ti(3).

The third substitution scheme was identical to the second scheme, however, the reaction was charge balanced with Fe³⁺.³⁹ The third substitution scheme followed the reaction,

CaZrTi₂O₇ + xMO₂ + xFe₂O₃ → Ca_1−xM_xZrTi_2−2xFe_2xO₇ + xCaO + 2xTiO₂

where M = Ce, Pu, Th, U.

3 Method

All the calculations in this work were performed with the General Utility Lattice Program (GULP).⁴¹ Ions are treated as charged spheres represented by their formal charge⁴² with a coulombic attraction/repulsion. Short range interactions between ions are described with Buckingham potentials of the form,
where r_ij is the distance between two ions i and j, and A, ρ and C₆ are parametrised constants specific to each interaction pair, as summarised in Table 1.

Table 1 Force field parameters for Buckingham potentials used in this work

Interaction	A (eV)	ρ (Å)	C6 (eV Å^6)	Ref.
O^2−–O^2−	25.410	0.6937	32.320	48
Ca^2+–O^2−	2272.741	0.2986	0.000	48
Zr^4+–O^2−	7290.347	0.2610	0.000	49
Ti^4+–O^2−	4545.823	0.2610	0.000	49
Al^3+–O^2−	2409.505	0.2649	0.000	48
Fe^3+–O^2−	3219.335	0.2641	0.000	48
Ce^4+–O^2−	2409.505	0.3260	0.000	49
Pu^4+–O^2−	752.224	0.4007	0.000	24
Th^4+–O^2−	8638.5	0.2856	70.000	50
U^4+–O^2−	9296.65	0.2796	90.000	50

---

The polarisability of the system is described by the shell model⁴³ where the charged core interacts with a massless “shell” via a spring constant, k. Only O²⁻ was polarised in this study so the interaction potentials presented in Table 1 are cation–anion core–shell interaction potentials. The shell model data for O²⁻ are presented in Table 2.

Table 2 Shell parameters for O²⁻ used in this work

Effective charge (core/shell)	k (eV Å^−2)	Ref.
0.513/−2.513	20.53	51

---

The model zirconolite lattice is the stoichiometric 2M polytype taken from the work of Gilbert et al.,²⁴ based on the structure published by Rossell.¹⁷ For our simulations the structure was expanded to a 2 × 2 × 2 supercell (704 atoms). The calculations were performed at constant pressure and the structure and atomic positions were optimised using a Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.^44–47

Substitution sites were randomly chosen. Our base zirconolite-2M system was composed of 704 atoms, 64 of those were Zr⁴⁺, 64 were Ca²⁺, 128 were Ti⁴⁺ (64 Ti(1), 32 Ti(2), 32 Ti(3)). For the first solid solution scheme, we chose to produce 30 random site substitution configurations. For the second solid solution scheme we produced 120 random configurations per site combination, e.g. CeTi(1)Ti(1). That is, for every substitution introduced to the system we replaced one Ca²⁺ site at random and two of the chosen type of Ti⁴⁺ sites at random with Al³⁺. For the third solid solution scheme, our aim was to directly compare Al³⁺ and Fe³⁺ as charge balancing species. Therefore, we made a direct substitution of Fe³⁺ on the sites that were occupied by Al³⁺ in the second solid solution scheme. The structures remained identical, otherwise. In our simulations we considered substitution concentrations of 3, 6, 9, 12, 15, 18 and 21%. In practice, this was the percentage of the number of atoms of the element in the cell to be replaced by the substitution rounded to the nearest integer, which was 2, 4, 6, 8, 10, 12 and 13 sites for Zr⁴⁺ and Ca²⁺, and 4, 8, 12, 20, 24 and 26 sites for Ti⁴⁺, when considering Al³⁺ and Fe³⁺ charge compensation. For a cell of 704 atoms, there were 64 Zr⁴⁺ and 64 Ca²⁺ sites that were potential substitution sites. Therefore, we were able to sample a wide range of configurations without artificial symmetry restrictions and provide data that could be compared to experimental investigation. Importantly, the upper substitution limit in this study concentration approximates that required for wasteform deployment.

Once the groundstate energy for each randomly substituted lattice was obtained we calculated the solution energy, that is the enthalpy of solid solution, for each substitution scheme. The calculation was based on the reactions presented in Section 2. We obtained the energy of each component by performing a geometry optimisation calculation on each structure, we assumed the polymorph of TiO₂ to be rutile.

Configurations that failed to optimise were removed from the spread of data. The solution energies were averaged. For an optimised spread the ground state energy differed by 2 eV to 3 eV. Each point on the graphs presented in Fig. 3 is the average solution energy for the denoted M⁴⁺ concentration. Each line presented in Fig. 4 and 5 is the trend in the average solution energy over all substitutions for each scheme. Where the trend in the solution energy of a substitution did not follow the global trend in the scheme, it is presented separately on the same graph with a dashed line.

To perform the above randomisation of site substitution in structures for each solid solution scheme, we wrote Python based software to enable the automation of the generation of the structures that are then passed to a high-performance computer server, enabling us to rapidly perform simulations and, from the following analysis, provide meaningful suggestions for material synthesis. An example time frame would be about 1–12 hours per simulation on a single core, with simulations queued as an array job; this resulted in about 2 days for the 120 simulation cell analyses from queueing to data clean-up and analysis. The above workflow is illustrated in Fig. 2.

Fig. 2 Flowchart demonstrating the high throughput methodology described in the Method section.

Importantly, Fig. 2 extends the work performed in this study to a pipeline workflow for high-throughput materials discovery. Enclosed within the box is the work performed within this paper. Outside of the box are future steps that can be performed where more expensive computational methods are applied. This demonstrates the method of use computationally cheaper methods – force fields – to scan the search space of crystal structure to quickly provide targets for investigation by more expensive methods. Within the context of the work here, this is beneficial when it is impractical to physically synthesise all possible solid solution stoichiometries. Within the wider context of materials science, our method has application to the exploration of similar materials e.g. high entropy alloys, capacitor ceramics, perovskites. Furthermore, the use of this initial screening can then direct subsequent, more costly, ab initio simulations, where the initial screening method has narrowed down the range of target compositions. The full calculation outputs are presented in the associated ESI.†

4 Results

4.1 CaZr_1−xM_xTi₂O₇ (M = Ce, Pu, Th, U)

Fig. 3 shows the mean solution energy for the following scheme, plotted against Ce, Pu, Th, and U concentration.

CaZrTi₂O₇ + xMO₂ → CaZr_1−xM_xTi₂O₇ + xZrO₂

Fig. 3 Comparison of solution energy as a function of M⁴⁺ substitution on the Zr⁴⁺ site (M = Ce, Pu, Th, U).

The solution energies of Ce⁴⁺ and Pu⁴⁺ substitution on the Zr⁴⁺ site become progressively more negative with increasing Ce⁴⁺ and Pu⁴⁺ concentration; the slope for Pu⁴⁺ is greater than for Ce⁴⁺. In contrast, for Th⁴⁺ substitution, an increase in solution energy is observed with increasing concentration. The solution energy for U⁴⁺ substitution initially follows a downward trend, however, beyond a U concentration of 12%, there is an abrupt increase in solution energy. Thereafter, the solution energies for U⁴⁺ substitution follow a downward trend despite the discontinuity in solution energy from 12–15% concentration.

4.2 Ca_1−xM_xZrTi_2−2xAl_2xO₇ (M = Ce, Pu, Th, U)

Results for Ce⁴⁺, Pu⁴⁺, Th⁴⁺ and U⁴⁺ substitution on the Ca²⁺ site and charge balance by replacement of Ti⁴⁺ with Al³⁺, are presented in Fig. 4. The mean solution energies are plotted against M⁴⁺ substitution concentration, with Al³⁺ as the charge balancing species substituted on sites Ti(1)Ti(1), Ti(2)Ti(2), Ti(3)Ti(3), Ti(1)Ti(2), Ti(1)Ti(3), and Ti(2)Ti(3). The solution energies of particular substitution schemes which do not follow the general trends are plotted separately with their own dashed lines. There is an upward trend in solution energy with increasing Ce, Pu, Th and U concentration for all schemes involving substitution on either the Ti(1) or Ti(3) sites, and the solution energy is always positive above an M⁴⁺ concentration of 6%. The solution energy only decreases with increasing M⁴⁺ concentration when charge balancing with Al³⁺ on the Ti(2)Ti(2) site combination, for which the solution energy is always negative.

Fig. 4 Comparison of solution energy as a function of M⁴⁺ substitution on the Ca²⁺ site with charge balance of Al³⁺ on Ti⁴⁺ site combinations, Ti(1)Ti(1), Ti(2)Ti(2), Ti(3)Ti(3), Ti(1)Ti(2), Ti(1)Ti(3), and Ti(2)Ti(3) sites (M = Ce, Pu, Th, U).

While Pu⁴⁺ substitution with Al³⁺ charge balance on two Ti(1) sites follows the general trend of the other substitutions, where the solution energy increases with increasing concentration, we observe maxima in the solution energy at 6%, 12% and 21% M⁴⁺ concentration. In the case of Al³⁺ charge balance on the combination of Ti(1) and Ti(3) sites, Ce⁴⁺ substitution did not follow the trend of Pu⁴⁺, Th⁴⁺ and U⁴⁺ so the data for it is plotted individually on the graph.

4.3 Ca_1−xM_xZrTi_2−2xFe_2xO₇ (M = Ce, Pu, Th, U)

Fig. 5 shows the results for substitution of Ce⁴⁺, Pu⁴⁺, Th⁴⁺ and U⁴⁺ substitution on the Ca²⁺ site of the zirconolite-2M structure with charge balance by replacement of Ti⁴⁺ with Fe³⁺. The mean solution energies are plotted against M⁴⁺ substitution concentration, with Fe³⁺ as the charge balancing species substituted on sites Ti(1)Ti(1), Ti(2)Ti(2), Ti(3)Ti(3), Ti(1)Ti(2), Ti(1)Ti(3), and Ti(2)Ti(3). Again, the solution energies of particular substitution schemes which do not follow the general trends are plotted separately with their own dashed lines.

Fig. 5 Comparison of solution energy as a function of M⁴⁺ substitution on the Ca²⁺ site and charge balance of Fe³⁺ on Ti⁴⁺ site combinations, Ti(1)Ti(1), Ti(2)Ti(2), Ti(3)Ti(3), Ti(1)Ti(2), Ti(1)Ti(3), and Ti(2)Ti(3) sites (M = Ce, Pu, Th, U).

The general trends in solution energy with Fe³⁺ as a charge balancing species are similar to those of Al³⁺. All substitution schemes showed an increase in solution energy with increasing M⁴⁺ concentration, with two exceptions. In the case of Ce⁴⁺ substitution with Fe³⁺ charge balancing on the Ti(1)Ti(3) site, the solution energy decreased with increased Ce⁴⁺ concentration, however, we observe a large increase in solution energy in the compositional interval between 18–21% Ce⁴⁺ incorporation. Whereas, in the case of U⁴⁺ substitution with Fe³⁺ charge balancing on Ti(2)Ti(3) sites, the solution energies for each compositional interval are much lower than for counterpart M⁴⁺ substitutions.

5 Discussion

The negative solution energies presented in Fig. 3 suggest that zirconolite-2M may fully accommodate Ce⁴⁺, U⁴⁺, Th⁴⁺ and Pu⁴⁺ on the Zr⁴⁺ site at low to moderate concentrations, which is in broad agreement with experimental validations for corresponding CaZr_1−xM_xTi₂O₇ solid solutions (M = Ce, U, Th, Pu). Furthermore, as the concentration of substitution is increased, the mixing of Ce⁴⁺, U⁴⁺ and Pu⁴⁺ is increasingly favoured, but only up to a value of 15% in the case of U⁴⁺, where we observe a discontinuity in solution energy. The observed discontinuity in solution energy for U⁴⁺ is consistent with the apparent solid solution limit of U in the zirconolite-2M structure as reported by Vance et al.¹³ Transformation to the zirconolite-4M polytype structure was reported in excess of approximately 15% U⁴⁺ substitution in the Zr⁴⁺ site of the zirconolite-2M structure.¹³ The 4M polytype also crystallises in the space group C2/c and is commonly described as an intergrowth of zirconolite-2M and pyrochlore-type layers, parallel to the [001] axis, resulting in a doubling of the unit cell.⁵² The zirconolite-4M phase remains the dominant structure in the CaZr_1−xU_xTi₂O₇ system up to a value of approximately 40% substitution, after which the cubic pyrochlore CaUTi₂O₇-type structure is preferred. Similar solid solution limits for Ce in the corresponding CaZr_1−xCe_xTi₂O₇ system were reported by Blackburn et al.⁵³ and Begg et al.⁵⁴ with the Ce inventory preferentially accommodated in the zirconolite-4M structure above 20% incorporation. However, it must be recognised that the tendency of Ce⁴⁺ to undergo reduction to Ce³⁺, when processing under inert or reducing conditions, does not permit formation of zirconolite-4M in the same solid solution, rather a Ce-rich CaTiO₃ phase is preferentially formed.⁵³ Nevertheless, targeting equimolar Ce³⁺ substitution between Ca²⁺ and Zr⁴⁺ sites, i.e. Ca_1−xZr_1−xCe_2xTi₂O₇ was observed to yield a transformation to zirconolite-4M.⁵⁵ Begg et al. fabricated the CaZr_1−xPu_xTi₂O₇ solid solution confirming that Pu⁴⁺ was preferentially accommodated in the 4M structure at around 15% incorporation, consistent with data for Ce⁴⁺ and U⁴⁺.⁵⁶ Consequently, simulation studies of substituent reduction and polytype transitions of zirconolite are necessary.

It follows that a similar trend would be expected for Pu in the data presented in Fig. 3, however, a continuous trend of negative solution energy was observed. It should be noted that a number of configurations did fail to optimise in our simulations, suggesting that certain defect arrangements are highly unfavourable. This may correspond to experimental observations, in which the 2M structure becomes less favourable towards high substitution concentrations, possibly due to substituent proximity within a lattice. Our observation that Ce⁴⁺ and Pu⁴⁺ substitute favourably for Zr⁴⁺ in zirconolite is further supported by the observations of Gilbert et al.²⁴ These data indicate that, whilst Ce remains a safe and practical analogue for Pu in wasteform development trials, it cannot fully replicate the substitution behaviour of Pu in zirconolite. Despite a similar trend to Pu⁴⁺ and U⁴⁺ at low concentrations, a clear variation in the solution energy, as a function of substitution, was observed. Nevertheless, the limitations of Ce–Pu surrogacy have been previously discussed in the context of Pu immobilisation in ceramic materials.^57,58

In contrast, the substitution of Th⁴⁺ for Zr⁴⁺ produces a continuous positive upward trend in solution energy, which becomes positive above 9% Th substitution, suggesting that Th⁴⁺ may have a narrow solid solution range in the zirconolite-2M structure. These data are consistent with recent observations by Blackburn et al.⁵⁹ in which it was confirmed that the solubility of Th⁴⁺ in the CaZr_1−xTh_xTi₂O₇ solid solution was limited to 10% substitution for Zr⁴⁺, with Th⁴⁺ preferentially accommodated in a pyrochlore-structured phase between 0.10 ≤ x ≤ 0.50. The single phase pyrochlore compound CaZr_0.40Th_0.60Ti₂O₇ was produced when targeting x = 0.60. Interestingly, a phase transition to the zirconolite-4M structure, as reported in analogue Ce and U solid solutions, was not observed.

The data presented in Fig. 4 and 5 demonstrate that the substitution of Ce⁴⁺, U⁴⁺, Th⁴⁺ and Pu⁴⁺ in the Ca²⁺ site, with charge balance provided by Al³⁺ and Fe³⁺ is favoured at M⁴⁺ concentrations around 3%. Yet, these solid solutions become rapidly unfavourable, tending towards positive solution energy with the exception of the substitution scheme in which charge compensators were accommodated in the Ti(2) site. This substitution scheme follows a downward trend in solution energy, suggesting that surrogate species may be accommodated in the Ca²⁺ site up to a 21% substitution, with Al³⁺ and Fe³⁺ preferentially accommodated in the Ti(2) site, consistent with some experimental observations. Loiseau et al.⁶⁰ fabricated the Ca_1−xNd_xZrTi_2−xAl_xO₇ solid solution, confirming that zirconolite-2M was produced as a single phase in the compositional range x ≤ 0.60, with Nd³⁺ deployed as a trivalent actinide surrogate, and further substitution resulted in the formation of the orthorhombic 3O polytype. Rietveld refinement of a zirconolite-2M structural model, in which Al³⁺ was constrained in the Ti(2) site, was refined against powder X-ray diffraction data for Ca_0.7Nd_0.3ZrTi_1.7Al_0.3O₇, confirmed that Al³⁺ preferentially occupied this site relative to Ti(1) and Ti(3). Similarly, Fe³⁺ was demonstrated by Whittle et al.³⁷ to substitute for Ti(2) in the CaZrTi_2−2xNb_xFe_xO₇ solid solution, however, it must be recognised that no surrogate as targeted to replace Ca²⁺. Conversely, Fe K-edge XANES has failed to resolve any preferential occupation of Fe³⁺ between Ti(1)/Ti(3) and Ti(2) sites in the Ca_1−xHo_xZrTi_2−xFe_xO₇.⁶¹ Similarly, Forder et al.⁶² resolved Fe³⁺ coordination in single phase zirconolite-2M (targeting Ca_1−xCe_xZrTi_2−2xFe_2xO₇) using ⁵⁷Fe Mössbauer spectroscopy, confirming that whilst Fe³⁺ occupied both 5- and 6-fold coordination Ti⁴⁺ sites, occupation of the Ti(2) site was preferred at low Fe³⁺ concentration. Previous simulation studies did not identify a significant preference of Fe³⁺ in any particular Ti⁴⁺ site of zirconolite-2M.²⁴ Cr³⁺ coordination was probed by Blackburn et al. in the Ca_1−xCe_xZrTi_2−2xCr_2xO₇ system, with deconvolution of the pre-edge Cr K-edge XANES region consistent with Cr³⁺ arranged in 6-fold coordination, inferring occupation in the octahedral Ti(1) and Ti(3) sites, as expected from consideration of crystal field stabilisation energy.⁶³

The variation in dominant charge compensation mechanism and preferential site occupancy of charge compensation species in zirconolite-2M may be attributed to a several factors, not limited to: the choice of Pu surrogate deployed, the valence state, electronic structure and ionic radius of the charge compensation cation, and the partial oxygen pressure imposed during the fabrication route, which has been systematically demonstrated to influence both surrogate oxidation state and partitioning in the zirconolite-2M structure. This work also provided some evidence that, in line with experimental observations, Pu⁴⁺ may be favourably accommodated on the Ca²⁺ site in the zirconolite structure, whereas substitution of Ce⁴⁺ may be unfavourable in some instances. However, it is important to note that the underlying mechanism constraining the site occupancy of Ce within zirconolite-2M is controlled by the prevailing redox conditions imposed during synthesis, and is not entirely governed by the chosen solid solution regime. Whereas Pu⁴⁺ is readily incorporated in the zirconolite-2M phase under oxidising conditions, Ce⁴⁺ has a tendency to partially reduce to Ce³⁺ regardless of sintering environment. For example, synthesis of the CaZr_1−xCe_xTi₂O₇ and Ca_1−xCe_xZrTi_2−2xCr_2xO₇ solid solutions in air consistently resulted in partial reduction of the Ce⁴⁺ inventory to Ce³⁺. It is this underlying auto-reduction tendency, that does not present itself with Pu under such conditions, that is the limiting factor in Ce–Pu surrogacy. Near single phase zirconolite-2M materials with nominal composition Ca_0.8Pu_0.2HfTi_1.6Al_0.4O₇ and Ca_0.8Pu_0.2ZrTi_1.8Al_0.2O₇ have been previously reported by both conventional sintering and hot pressing techniques, targeting Pu⁴⁺ and Pu³⁺ respectively.^20,64 Analysis of the compositionally analogous Ca_1−xCe_xZrTi_2−2xAl_2xO₇ solid solution, at x = 0.10, has been observed to produce a minor perovskite phase, attributed to partial Ce³⁺ speciation-2M.⁶⁵ The differences between simulated Ce⁴⁺ and Pu⁴⁺ substitution behaviour suggests that clustering and localised substituent–substituent interactions, which would be present in the current work yet excluded from previous data,²⁴ may be key for stabilising the substituents. The combined presence of numerous defects may relieve the localised stress they create, as opposed to lone, or few, defects.²⁸ This may also explain the data presented by Ji et al.⁶¹ where the incorporation of Ln³⁺ species in the Ca²⁺ site would lead to a weaker binding energy with Fe³⁺ defects distributed across the Ti⁴⁺ sites, and thus preferential occupation of Fe³⁺ within any specific Ti⁴⁺ site was not reported. It remains clear that in the present study, and as has been confirmed in a selection of laboratory investigations, that low valence charge balancing cations preferentially occupy the trigonal biprymidal TiO₅ site. It follows that the partially occupied nature of this site (50% occupied), relative to Ti(1) and Ti(3), permits the accommodation of cations of varied size.

From examination of the collective substituent behaviour we establish that at low concentrations (e.g. 3%) a selection of charge compensated substitution schemes are viable within the zirconolite-2M phase, and it has been experimentally determined that environmental conditions, chemical activity and ion mobility may dictate the solid solution mechanisms that occur. As the nominal concentration of substitution is increased, it is to be expected that clustering of defects will occur in either the Zr⁴⁺ site or Ca²⁺ site, facilitated by charge compensation on the Ti(2) site. The investigation presented herein shows that substitution schemes involving the Ca²⁺ site reach lower solution energies when charge balancing with Al³⁺ on the Ti(2) site, this requires further computational work to elucidate, an interpretation is that Al³⁺ being smaller than Fe³⁺ can be accommodated more easily. Nevertheless, these data support the deployment of zirconolite-2M as a potential single host phase for the immobilisation of Pu oxides, and are in general agreement with a selection of recent publications concerning the solid solution behaviour of Ce, U, Th and Pu.

6 Conclusions

We have demonstrated trends in the energetics of zirconolite-2M solid solutions with Ce⁴⁺, Pu⁴⁺, Th⁴⁺, and U⁴⁺ cations that are in general agreement with published experimental data concerning the deployment of zirconolite-2M as a host for actinides. Consequently, we have shown that atomistic simulations can effectively guide the formulation development of these materials, and inform experimental validation. For example, using the high throughput methodology developed and reported here, it is possible to rapidly screen many solid solution schemes in silico and evaluate their relative stability, to guide more resource intensive ab initio simulations and laboratory investigation for validation. Indeed, this method should have much wider utility in the exploration in the optimisation of functional materials such as high entropy alloys, capacitor ceramics, and perovskite catalysts. Our investigations have also shown that Ce is not a direct analogue for the actinide cations such as Pu, as has been validated in a number of wasteform development trials.

Author contributions

S. D. – data curation, investigation, formal analysis, methodology, software, visualization, writing – original draft, writing – review & editing. C. M. H. – conceptualization, formal analysis, investigation, methodology, software, writing – review & editing. L. R. B. – investigation, validation, writing – review & editing. C. L. F. – project administration, resources, supervision, writing – review & editing. N. C. H. – conceptualization, funding acquisition, project administration, resources, supervision, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We are grateful for financial support from the Nuclear Decommissioning Authority and EPSRC under grant reference EP/S01019X/1 and EP/R511754/1. This research utilised the HADES/MIDAS facility at the University of Sheffield established with financial support from EPSRC and BEIS, under grant EP/T011424/1.⁶⁶ The computational work was made possible by the resources provided by Research-IT at the University of Sheffield. We would also like to extend our gratitude to Ondřej Krejčí for invaluable input on the presentation of this work.

Notes and references

1. N. C. Hyatt, npj Mater. Degrad., 2020, 4, 28.
2. W. J. Weber, R. C. Ewing, C. A. Angell, G. W. Arnold, A. N. Cormack, J. M. Delaye, D. L. Griscom, L. W. Hobbs, A. Navrotsky, D. L. Price, A. M. Stoneham and M. C. Weinberg, J. Mater. Res., 1997, 12, 1948–1978.
3. W. J. Weber, R. C. Ewing, C. R. A. Catlow, T. D. de la Rubia, L. W. Hobbs, C. Kinoshita, H. Matzke, A. T. Motta, M. Nastasi, E. K. H. Salje, E. R. Vance and S. J. Zinkle, J. Mater. Res., 1998, 13, 1434–1484.
4. R. C. Ewing, Earth Planet. Sci. Lett., 2005, 229, 165–181.
5. R. C. Ewing, Prog. Nucl. Energy, 2007, 49, 635–643.
6. N. C. Hyatt and M. I. Ojovan, Materials, 2019, 12, 3611.
7. R. C. Ewing, C. R. Geosci., 2011, 343, 219–229.
8. A. E. Ringwood, S. E. Kesson, N. G. Ware, W. Hibberson and A. Major, Nature, 1979, 278, 219–223.
9. W. Lutze and R. C. Ewing, Radioactive waste forms for the future, North-Holland, 1988.
10. K. P. Hart, E. R. Vance, M. W. A. Stewart, J. Weir, M. L. Carter, M. Hambley, A. Brownscombe, R. A. Day, S. Leung, C. J. Ball, B. Ebbinghaus, L. Gray and T. Kan, MRS Proceedings, 1997, 506, 161.
11. K. Zhang, B. Luo and H. Zhang, Mater. Res. Express, 2019, 6, 115526.
12. R. F. Haaker and R. C. Ewing, Naturally occurring crystalline phases: analogues for radioactive waste forms, Richland, WA (United States), 1981.
13. E. R. Vance, G. R. Lumpkin, M. L. Carter, D. J. Cassidy, C. J. Ball, R. A. Day and B. D. Begg, J. Am. Ceram. Soc., 2002, 85, 1853–1859.
14. X. Deschanels, V. Picot, B. Glorieux, F. Jorion, S. Peuget, D. Roudil, C. Jégou, V. Broudic, J. N. Cachia, T. Advocat, C. Den Auwer, C. Fillet, J. P. Coutures, C. Hennig and A. Scheinost, J. Nucl. Mater., 2006, 352, 233–240.
15. L. R. Blackburn, D. J. Bailey, S.-K. Sun, L. J. Gardner, M. C. Stennett, C. L. Corkhill and N. C. Hyatt, Adv. Appl. Ceram., 2021, 120, 69–83.
16. B. M. Gatehouse, I. E. Grey, R. J. Hill and H. J. Rossell, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1981, 37, 306–312.
17. H. J. Rossell, Nature, 1980, 283, 282–283.
18. M. C. Stennett, C. L. Corkhill, L. A. Marshall and N. C. Hyatt, J. Nucl. Mater., 2013, 432, 182–188.
19. B. D. Begg, E. R. Vance and G. R. Lumpkin, MRS Proceedings, 1997, 506, 79.
20. B. D. Begg and E. R. Vance, MRS Proceedings, 1996, 465, 333.
21. C. Dickson and F. Glasser, Cem. Concr. Res., 2000, 30, 1619–1623.
22. W. A. Lambertson, M. H. Mueller and F. H. Gunzel, J. Am. Ceram. Soc., 1953, 36, 397–399.
23. P. W. Tasker, Le Journal de Physique Colloques, 1980, 41, C6-488–C6-491.
24. M. Gilbert and J. H. Harding, Phys. Chem. Chem. Phys., 2011, 13, 13021.
25. L. Veiller, J.-P. Crocombette and D. Ghaleb, J. Nucl. Mater., 2002, 306, 61–72.
26. H. R. Foxhall, K. P. Travis, L. W. Hobbs, S. C. Rich and S. L. Owens, Philos. Mag., 2013, 93, 328–355.
27. H. R. Foxhall, K. P. Travis and S. L. Owens, J. Nucl. Mater., 2014, 444, 220–228.
28. J. Tanti and N. Kaltsoyannis, J. Nucl. Mater., 2021, 543, 152525.
29. J. A. Fortner, A. J. Kropf, A. J. Bakel, M. C. Hash, S. B. Aase, E. C. Buck and D. B. Chamberlain, MRS Proceedings, 1999, 608, 401.
30. J. A. Fortner, A. J. Kropf, R. J. Finch, A. J. Bakel, M. C. Hash and D. B. Chamberlain, J. Nucl. Mater., 2002, 304, 56–62.
31. C. Nästren, R. Jardin, J. Somers, M. Walter and B. Brendebach, J. Solid State Chem., 2009, 182, 1–7.
32. D. J. Gregg, Y. Zhang, S. C. Middleburgh, S. D. Conradson, G. Triani, G. R. Lumpkin and E. R. Vance, J. Nucl. Mater., 2013, 443, 444–451.
33. S. Finkeldei, M. C. Stennett, P. M. Kowalski, Y. Ji, E. de Visser-Týnová, N. C. Hyatt, D. Bosbach and F. Brandt, J. Mater. Chem. A, 2020, 8, 2387–2403.
34. R. L. Miller and G. A. Reimann, Thermodynamics of gas-metal-slag equilibria for applications in in situ and ex situ vitrification melts, 1993.
35. P. Loiseau, D. Caurant, N. Baffier, L. Mazerolles and C. Fillet, J. Nucl. Mater., 2004, 335, 14–32.
36. G. R. Lumpkin, K. R. Whittle, C. J. Howard, Z. Zhang, F. J. Berry, G. Oates, C. T. Williams and A. N. Zaitsev, MRS Proceedings, 2006, 932, 53.1.
37. K. R. Whittle, N. C. Hyatt, K. L. Smith, I. Margiolaki, F. J. Berry, K. S. Knight and G. R. Lumpkin, Am. Mineral., 2012, 97, 291–298.
38. F. J. Berry, G. R. Lumpkin, G. Oates and K. R. Whittle, Hyperfine Interact., 2006, 166, 363–366.
39. M. R. Gilbert, C. Selfslag, M. Walter, M. C. Stennett, J. Somers, N. C. Hyatt and F. R. Livens, IOP Conf. Ser.: Mater. Sci. Eng., 2010, 9, 012007.
40. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
41. J. D. Gale and A. L. Rohl, Mol. Simul., 2003, 29, 291–341.
42. M. Born and K. Huang, Dynamical Theory of Crystal Lattices, Clarendon Press, 1954.
43. B. G. Dick and A. W. Overhauser, Phys. Rev., 1958, 112, 90–103.
44. C. G. Broyden, Mathematics of Computation, 1970, 24, 365.
45. R. Fletcher, Comput. J., 1970, 13, 317–322.
46. D. Goldfarb, Mathematics of Computation, 1970, 24, 23.
47. D. F. Shanno, Mathematics of Computation, 1970, 24, 647.
48. T. S. Bush, J. D. Gale, C. R. A. Catlow and P. D. Battle, J. Mater. Chem., 1994, 4, 831.
49. S. M. Woodley, P. D. Battle, J. D. Gale and C. R. A. Catlow, Phys. Chem. Chem. Phys., 1999, 1, 2535–2542.
50. M. Nadeem, M. J. Akhtar, R. Shaheen, M. N. Haque and a Y. Khan, J. Mater. Sci. Technol., 2001, 17, 638–642.
51. C. R. A. Catlow, in Solid State Chemistry Techniques, Clarendon Press, 1987, ch. 7.
52. A. A. Coelho, R. W. Cheary and K. L. Smith, J. Solid State Chem., 1997, 129, 346–359.
53. L. R. Blackburn, S. Sun, L. J. Gardner, E. R. Maddrell, M. C. Stennett and N. C. Hyatt, J. Nucl. Mater., 2020, 535, 152137.
54. B. D. Begg, E. R. Vance, R. Day, M. Hambley and S. D. Conradson, MRS Proceedings, 1996, 465, 325.
55. C. Meng, X. Ding, W. Li, J. Zhao and H. Yang, J. Mater. Sci., 2016, 51, 5207–5215.
56. B. D. Begg, R. A. Day and A. Brownscombe, MRS Proceedings, 2000, 663, 259.
57. P. A. Bingham, R. J. Hand, M. C. Stennett, N. C. Hyatt and M. T. Harrison, MRS Proceedings, 2008, 1107, 421.
58. M. V. Zamoryanskaya and B. E. Burakov, MRS Proceedings, 2000, 663, 301.
59. L. R. Blackburn, S. Sun, L. J. Gardner, E. R. Maddrell, M. C. Stennett, C. L. Corkhill and N. C. Hyatt, J. Am. Ceram. Soc., 2020, 2.
60. P. Loiseau, D. Caurant, N. Baffier and C. Fillet, MRS Proceedings, 2002, 757, II6.3.
61. S. Ji, M. Su, C. Liao, S. Ma, Z. Wang, K. Shih, C. Chang, J. Lee, T. Chan and Y. Li, J. Am. Ceram. Soc., 2020, 103, 1463–1475.
62. S. D. Forder, P. A. Bingham, O. J. McGann, M. C. Stennett and N. C. Hyatt, Hyperfine Interact., 2013, 217, 83–90.
63. L. R. Blackburn, S.-K. Sun, S. M. Lawson, L. J. Gardner, H. Ding, C. L. Corkhill, E. R. Maddrell, M. C. Stennett and N. C. Hyatt, J. Eur. Ceram. Soc., 2020, 40, 5909–5919.
64. E. R. Vance, C. J. Ball, R. A. Day, K. L. Smith, M. G. Blackford, B. D. Begg and P. J. Angel, J. Alloys Compd., 1994, 213–214, 406–409.
65. L. R. Blackburn, S. K. Sun, L. J. Gardner, E. R. Maddrell, M. C. Stennett and N. C. Hyatt, MRS Adv., 2020, 5, 93–101.
66. N. C. Hyatt, C. L. Corkhill, M. C. Stennett, R. J. Hand, L. J. Gardner and C. L. Thorpe, IOP Conf. Ser.: Mater. Sci. Eng., 2020, 818, 012022.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1ra02914b

---