The interaction of Pu(IV) with the hematite (001) terminations: a periodic boundary condition DFT study

Abstract

The Enhanced Actinide Removal Plant (EARP), at Sellafield in the UK, is tasked with separating waste actinide species from an aqueous waste stream via base-induced hydrolysis of Fe(III). During this flocculation process ferrihydrite forms and the actinides interact strongly with the surface. It has been shown that Pu remains sorbed in the solid state over long periods of time, during which ferrihydrite undergoes transformation into hematite. It is of critical importance to the operations of future Geological Disposal Facility (GDF) technologies that the Pu@hematite system is studied to understand the binding strength and sorption mechanism. Here we present a comprehensive study of Pu(IV) binding to two well-established basal terminations of hematite using periodic DFT+U_eff. First, we outline our methodology and demonstrate correct prediction of the bulk hematite lattice parameters and electronic band gap, then we generate the (001)-Fe and (001)-O₃ terminations and demonstrate reasonable predictions of the surface energies, inter-layer spacings, and work functions. The (001)-Fe termination is then hydrated with a monolayer of water. We show that Pu(IV) forms multiple Pu–O bonds with both terminations at distances consistent with experimental EXAFS measurements. Density of states and charge density difference analysis reveals strong hybridisation between the Pu(5f) and O(2p) states supporting charge transfer as indicated by depleted charge density surrounding the Pu atom on the surface. Quantum Theory of Atoms in Molecules analysis shows that the Pu–O bonds are partially covalent, in agreement with our previous assessment of Pu bound to α-Fe₁₃, a prenucleation cluster to ferrihydrite (Fh), and Pu bound to ferrihydrite surfaces. The reaction energies for surface binding are significantly exothermic, even more so than was found in our previous analysis of the Pu@Fh surfaces, indicating that Pu(IV) should remain immobile and bound to hematite, particularly in oxygenated conditions as may be found in the GDF or subterranean environments.

---

Introduction

Hematite (α-Fe₂O₃, Hem) is the most thermodynamically stable iron (oxyhydr)oxide mineral and is found in many aquatic and terrestrial systems. It is a common redox-active interface in the environment and is responsible for the uptake of contaminants and cycling of iron which is critical to many environmental and biological processes.¹ Besides its environmental importance, hematite is also favoured for use industrially in catalytic and photoelectric technologies due to its non-toxicity, high natural abundance, and electronic properties.²

Hematite is a significant product in UK nuclear clean-up operations at Sellafield on the Northwestern coast in Cumbria. The Enhanced Actinide Removal Plant (EARP) acts to separate radionuclides from post-PUREX waste for further treatment and long-term storage. The EARP feed primarily consists of aqueous Fe(III) and radionuclides produced at the Thermal Oxide Reprocessing Plant (THORP), which closed in 2018, and the Magnox Reprocessing Plant, which closed in 2022. These plants are now undergoing post-operational clean-out (POCO) to remove radioactive contaminants; therefore, the EARP feed is diversifying to include chemicals used during POCO and which were not typically used during regular operation of those plants. Understanding the implication of this is crucial to the continued operation and safety of the EARP.

The EARP operates via a base-induced flocculation process whereby Fe(III) hydrolysis results in the formation of the highly insoluble ferrihydrite (Fh) phase which interacts strongly with radionuclides and precipitates out of solution. After multiple ultrafiltration steps, the solid waste form containing the separated radionuclides is removed and the clean aqueous stream is discharged into the Irish Sea. In recent years there has been a growing interest in the mechanism of Fh formation, and it was shown to occur via the Fe₁₃ Keggin as a prenucleation cluster under EARP conditions.³ In 2019, it was further shown that under EARP conditions a significant amount of Pu(IV) is removed from solution at very low pH, before Fh has formed, indicating that Pu(IV) may bind directly to the Fe₁₃ Keggin cluster.⁴ We studied this using Density Functional Theory (DFT), and showed that Pu(IV) indeed binds very strongly to the Fe₁₃ Keggin, as do all tetravalent actinides from Th–Pu, primarily through strong An⁴⁺–O²⁻ ionic interactions. However, overlap-driven covalent interactions were observed for Np and Pu showing that the 5f orbitals are important to understanding the bonding in these systems.⁵

Ferrihydrite transforms into hematite without impacting the Pu(IV) coordination environment or releasing significant amounts of the radionuclide, which suggests that Pu(IV) should be strongly bound to the surface of hematite. Pu(IV) is not observed to incorporate into the structure of ferrihydrite or hematite, as the Extended X-ray Absorption Fine Structure (EXAFS) fittings report four iron backscatterers at 3.34 Å,⁴ whereas if Pu(IV) were incorporated, the number of iron backscatters would be higher than four. The EXAFS fitting also shows eight oxygen backsatterers at 2.29 Å. We recently also studied the surface interactions of Pu@Fh and found strong ionic binding, reporting partially covalent Pu–O interactions in line with those found in the Pu–Fe₁₃ Keggin system at distances in good agreement with experiments.⁶

The UK is currently working towards long-term disposal of nuclear waste in a Geological Disposal Facility (GDF), where multiple barriers to radionuclide migration will be employed. Stainless steel canisters are one such barrier which, under the right conditions and on geological timescales, may corrode into ferrihydrite and goethite, which can then undergo phase transitions to hematite.^1,7,8 Given that post-EARP actinides are bound to hematite, that corrosion mechanisms can form hematite, and that the ground where the future GDF may be located is likely to contain hematite, it is very important to study the speciation of Pu-bound hematite. However, it is very difficult to do so experimentally owing to the radioactivity and dangers associated with plutonium. To our knowledge, nobody has attempted to study this complexation theoretically, though bulk and surface studies of bare hematite are plentiful.^9–20

In this work, we first establish a suitable theoretical methodology for bulk hematite, and then the Hem(001) surface under two different terminations. Then, having verified our approach, we introduce a hydrated Pu(IV) fragment to the surface to study the binding mode and electronic structure of Pu-bound hematite (Pu@Hem). The results show that Pu(IV) binds strongly to the Hem(001) terminations due to strong ionic interactions, but the Pu–O bonds are found to be partially covalent according to Quantum Theory of Atoms in Molecules (QTAIM) metrics. The multidentate binding shows Pu–O distances similar to those determined experimentally. The binding is found to be stronger than for Pu(IV) on ferrihydrite, supporting the experimental observation that Pu(IV) is retained during the ferrihydrite → hematite recrystallisation process.

Computational details

General considerations

All Density Functional Theory (DFT) calculations were conducted using the Vienna ab initio Simulation Package (VASP) versions 6.1.2, 6.3.0 and 6.4.1.^21–23 The generalised gradient approximation (GGA) functional PBE²⁴ was used in conjunction with Grimme's D3 dispersion correction.²⁵ It is well known that GGA functionals produce over-delocalised 3d states in iron (oxyhydr)oxide materials, resulting in significantly underestimated magnetic moments and band gaps. Therefore, it is common to apply a Hubbard correction to the Fe 3d states; U_eff = 3–5 eV is typically used for hematite and other similar iron (oxyhydr)oxide materials.^10–14 The Hubbard correction was applied according to the Dudarev formalism.²⁶ A comparison between lattice parameters and band gaps while varying U_eff = 0–5 eV can be found in the SI Fig. S1 and Table S1. As U_eff is increased the a and b lattice parameters increase and the c parameter decreases, and good agreement with experiment is found for U_eff = 4 eV. Plane wave basis sets were employed with core electrons modelled using Projector Augmented Wave (PAW) pseudopotentials.^27,28 These pseudopotentials were generated by the VASP developers by solving the all-electron scalar relativistic Schrödinger equation; therefore, results are limited to scalar effects without the inclusion of spin–orbit coupling. All calculations are spin polarised; initial magnetic moment of ±5μ_B are applied to the Fe(III) sites in the well-known antiferromagnetic (AFM) arrangement (↑↓↓↑, Fig. 1). Convergence testing was performed to determine that a plane wave cutoff of 600 eV is sufficient (Fig. S2 and Table S2). The total energy and forces were converged to <10⁻⁶ eV and <10⁻³ eV Å⁻¹ respectively.

Fig. 1 Antiferromagnetic (AFM) hematite bulk unit cell. The direction of the spin polarisation on each Fe centre is shown by yellow arrows. Colour scheme: Fe, orange; O, red.

Bulk

Convergence with respect to the Γ-centred k-point mesh was carried out (Fig. S2 and Table S3), and the data presented in the results and discussion section correspond to geometric optimisation using a 9 × 9 × 3 mesh with Gaussian smearing, whereas results presented for density of states use a larger 13 × 13 × 7 mesh and the tetrahedron method with Blöchl corrections.²⁹ The electronic band structure was calculated in line-mode with k-points defined to represent high symmetry points in the Brillouin zone and the results were analysed using sumo, a Python toolkit for plotting ab initio solid-state data.³⁰ The phonon band structure was calculated using Phonopy³¹ with a 2 × 2 × 1 supercell expansion (120 atoms). The phonon density of states was computed by interpolating the phonon frequencies onto a 16 × 16 × 16 q-point grid. The bulk structure is a stable minimum as indicated by all phonon frequencies being real (Fig. S3).

Surfaces

METADISE³² was employed to cut the relaxed hematite bulk to form dipole free terminations; this method has previously been used successfully on iron minerals.^6,33–35 A dipole correction was applied in the c direction to minimise the effects of surface dipole formation upon relaxation.³⁶ A vacuum region of 20 Å ensures that the slab does not interact with the periodic image of itself, even after molecular adsorption. Convergence with respect to the Γ-centred k-point mesh was carried out (Tables S4 and S5), and results presented for geometric optimisation use a 7 × 7 × 1 mesh with Gaussian smearing, whereas results presented for density of states use a larger 9 × 9 × 1 mesh and the tetrahedron method with Blöchl corrections. The surface models were then relaxed layer by layer, with layers identified using VASPKIT,³⁷ until the energy is converged to <1 meV per atom (Tables S6 and S7).

In this model, one side of the slab is relaxed while the other is kept fixed. The equations below represent energetic differences between relaxed and unrelaxed stochiometric surfaces with respect to the relaxed bulk and account for asymmetric slab relaxation. The unrelaxed surface energy is obtained by performing a single point calculation on the unrelaxed slab and is calculated according to eqn (1),

and the relaxed surface energy is calculated according to eqn (2),where γ_u and γ_r are the surface energies of the unrelaxed and relaxed surfaces, E_slab,u/r and E_bulk are the electronic energies of the (un)relaxed slab and relaxed bulk, and A is the surface area of the slab.^34,35 The hydrated surface energy was then calculated according to eqn (3),where E_slab+nH2O,r is the electronic energy of the relaxed hydrated slab and E_H2O is the energy of an isolated H₂O molecule.³⁴ On the Fe-termination, multiple initial arrangements of water molecules with 100% coverage were tested, and the optimised geometries and energies are presented in the SI (Fig. S4 and Table S8).

After determining stable surface terminations, a hydrated Pu(IV) species with the formula [Pu(H₂O)₅]⁴⁺ was placed in multiple initial configurations above the surfaces and the geometries were optimised. A Hubbard correction of U_eff = 4 eV was applied to the Pu 5f states for the same reasons as given for Fe 3d.^10–13 To our knowledge, our previous work is the only DFT+U_eff study involving Pu(IV) adsorbed to iron (oxyhydr)oxide mineral surfaces.⁶ In that work we applied U_eff = 4 eV to the Pu 5f based on its success in predicting the structural and magnetic properties of bulk and surface Pu oxides,^38–42 and in doing so verified that it also performs well in the context of contaminant adsorption studies.

Substitution reaction energies are calculated according to the following equations,

Hem + [Pu(H₂O)₉]⁴⁺ → Pu@Hem + 4H₂O (4)

ΔE_r = [E_Pu@Hem + 4E_H2O] − [E_Hem + E_{[Pu(H2O)9]⁴⁺}] (5)

where the nine coordinate Pu species represents Pu(IV) in aqueous EARP conditions.^5,43 Single point energy calculations were then performed including implicit solvation using VASPsol (ε = 78.4, water).^44,45

Analysis of the charge density was carried out using charge density difference and the Quantum Theory of Atoms in Molecules (QTAIM). The charge density difference is calculated according to.

which is the difference between the electron density of the Pu complex and the sum of the hydrated surface and adsorbed species, calculated separately while fixed in the optimised positions found in the complexed surface model. The QTAIM calculations were performed using CRITIC2 which analyses the topology of the electron density in periodic solids.⁴⁶ CRITIC2 requires the all-electron density, which was obtained by performing a single-point energy calculation with the setting LAECHG = .TRUE. and then summing the output core density (AECCAR0) and valence electron density (AECCAR2).

Results and discussion

Structural parameters and electronic structure

As mentioned in the methodology section, hematite displays a well-known AFM (↑↓↓↑) spin arrangement, shown in Fig. 1, and so no other spin arrangements were considered. The calculated magnetic moment on each Fe atom was the same ±4.12μ_B, in good agreement with previous calculation using the DFT+U_eff methodology (Table 1) and work using the hybrid functional HSE06, which yields ±4.16μ_B.⁹ With U_eff = 4 eV, the calculated lattice parameters for hematite are a = b = 5.05 Å and c = 13.84 Å, in good agreement with both previous calculation and experimental measurements (Table 1).

Table 1 Calculated lattice parameters, unit cell volume, magnetic moments, and band gap for the relaxed hematite bulk structure compared to previous theoretical and experimental work. The percentage errors are calculated using the average of the two experimental references

Parameter	Theoretical				Experimental
	This work	Ref. 10	Ref. 11	Ref. 12	Ref. 47	Ref. 48
a Range taken from ref. 49–51.
a = b/Å	5.05 (+0.3%)	5.03 (+0.0%)	5.07 (+0.7%)	5.02 (−0.3%)	5.035	5.031
c/Å	13.84 (+0.6%)	13.74 (−0.1%)	13.88 (+0.9%)	13.66 (−0.7%)	13.747	13.766
c/a	2.74	2.73	2.74	2.72	2.731	2.736
Volume/Å^3	305.67 (+1.3%)	301.42 (−0.1%)	308.98 (+2.4%)	298.12 (−1.2%)	301.76	301.75
μ/μB	±4.12	±4.18	±4.11	±4.23	—
Eg/eV	2.16	2.09	2.00	2.10	2.00–2.20^a

---

In this work, we found that the value of U_eff has little impact on the calculated lattice parameters; as U_eff is varied from 0–5 eV, a = b shows a range of 4.99–5.06 Å and c shows a range from 13.81–13.86 Å. In all cases there is only a small difference compared with experiments, and one could argue that any value of U_eff within the studied range is acceptable on the basis of the structural properties. However, the value of U_eff has a significant impact on the calculated magnetic moments and band gap. It is well documented that hematite displays an experimental band gap of approximately 2 eV (Table 1). The calculated band gap with U_eff ≤ 3 eV is too small, and the value with U_eff = 5 eV is too large (Fig. S1 and Table S1). With U_eff = 4 eV the band gap was calculated as 2.16 eV, in good agreement with experiments and previous calculations (Table 1).

The electronic band structure is shown in Fig. 2, and the projected density of states (pDOS) reveals that the valence band is dominated by the O(2p) states showing hybridisation with the Fe(3d) states. The conduction band is dominated by the Fe(3d) states and therefore hematite displays a charge transfer transition across the band gap which agrees with how hematite is generally considered to be a charge transfer semiconductor rather than Mott–Hubbard insulator.^13,14 This is also consistent with O K-edge X-ray absorption and emission spectroscopy which confirms that the valence band consists of strongly hybridised O(2p) and Fe(3d) states with primarily O(2p) character.⁵²

Fig. 2 Calculated electronic band structure and projected density of states for bulk hematite, with U_eff = 4 eV. Solid yellow bands are spin-up and dashed blue bands are spin-down. The valence band maximum, at 0 eV, and conduction band minimum are marked with a green and red dot respectively.

Surface modelling

The (001) basal surface was selected as a good candidate for Pu(IV) sorption for multiple reasons. Firstly, the (001) surface is stable and is often observed to form in natural crystals of hematite and other corundum-type oxides.^17,53–55 Secondly, there are a significant number of theoretical and experimental investigations of the structure and sorption of contaminants onto the (001) surface to compare results with. Finally, the (001) surface is known to display multiple terminations dependent on the conditions,^18,20 allowing for comparison of different surface arrangements to get a wider view of how Pu(IV) may bind under EARP conditions and potentially under conditions altered over GDF timescales.

Two terminations were selected: the first is Fe rich and second is O rich. These terminations are typically referred to as the “Fe–O₃–Fe” and “O₃–Fe–Fe” terminations to denote the top three surface layers, and hereafter we will refer to these as the Fe and O₃ terminations, or “-term” for short (Fig. 3). The Fe-term is non-polar and stable under lower oxygen pressures, and the O₃-term is polar and stable under a wide range of oxygen pressures.^{11–14,17,20}

Fig. 3 Interlayer spacing in the relaxed surface slabs, relative to the bulk, compared to experiments and previous calculations.

As with any theoretical surface study, the surface slab must be relaxed layer by layer until convergence of the surface energy is achieved; this is shown in Tables S6 and S7. The surface energies of the Fe and O₃ terminations were calculated to be 1.60 J m⁻² and 3.46 J m⁻² respectively, which fall within the range of values calculated previously (Table 2). The same can be said for the work function of the two surfaces, which were calculated as 4.33 eV and 8.47 eV respectively (Table 2).

Table 2 Calculated Hem(001) surface energies and work functions compared to previous work

	γ/J m^−2						Φ/eV
Surface termination	This work	Ref. 11	Ref. 12	Ref. 15	Ref. 16	Ref. 56	This Work	Ref. 17	Ref. 18	Ref. 57
Fe	1.60	1.8	1.66	2.0	1.78	1.0	4.33	4.3	4.7	4.0
O3	3.46	2.4	2.59	4.0	2.63	3.7	8.47	7.6	8.5	7.6

---

Our surface properties are in good agreement with previous simulation, but to our knowledge the experimental surface energies have not been measured. However, another way to determine the viability of our surface slab model is to consider the interlayer spacing, which has been measured experimentally.^54,55 Fig. 3 shows the interlayer spacing as a percentage difference to the unrelaxed slab (i.e. the bulk-like distances) compared to experiments and previous calculations alongside ball-and-stick representations of the surface slabs. The absolute values of the interlayer spacings for the two terminations are provided in Tables S9 and S10. One can see that there is some difference between the relaxation in the top layer, but good agreement is found for all other interlayer spacings (Fig. 3). Experimentally, the O₃-term displays a strong contraction of the top layer of −67%, which is not observed theoretically. This may be due to the limited thickness of the polar slab compared to the size of the bulk in a real crystal.

The pDOS for the two terminations are shown in Fig. 4. Much like the bulk, both surfaces act as charge transfer semiconductors with the valence band consisting of hybridised O(2p) and Fe(3d) states and the conduction band being primarily Fe(3d). The band gap at the surfaces is reduced due to states appearing within the gap. These surface states have been observed before,¹⁹ but no explanation was given as to their source. These states are typically attributed to undercoordinated surface atoms, particularly those at the bottom of the slab which are kept frozen during relaxation. By comparing states associated with the surface atoms to the overall pDOS (Fig. S5 and S6) we confirmed that the states around 1.3 eV in the Fe-term are indeed caused by the undercoordinated surface Fe atoms. The same analysis of the O₃-term shows that the O(2p) states at the Fermi level are surface states, but the Fe(3d) states could not be strictly assigned to the surface.

Fig. 4 Calculated projected density of states for the Fe and O₃ terminated Hem(001) surface.

Smith et al. observed a shell of eight oxygen backscatters close to Pu(IV),⁴ which is not attainable on the bare Fe-term. Therefore, as in our work on ferrihydrite,⁶ we chose to hydrate the Fe-term and study the binding of Pu(IV) to the hydrated surface. A monolayer of H₂O in a variety of different coordination positions with 100% coverage was relaxed (Fig. S4 and Table S8) and the lowest energy surface was used to bind Pu(IV). The bare O₃-term surface is sufficient to obtain sensible coordination numbers and so further hydration was not studied. Note, however, that the Pu(IV) fragment is explicitly hydrated, and therefore the impact of explicit H₂O is not completely omitted in the O₃-term calculations.

Hematite surface-Pu(IV) complexation

For each termination, a variety of starting geometries was optimised and only the lowest energy surface complexes are shown in Fig. 5. Higher energy complexes are shown in Fig. S7 and S8 and data for all complexes are provided in Tables S11 and S12. The most stable surface complexes have the lowest energy and therefore the most negative reaction energy ΔE_r associated with the binding of Pu(IV) to the surface (Table 3).

Fig. 5 Lowest energy Pu@Hem(001) complexes viewed down the a-axis (left) and b-axis (right). The closest eight O and four Fe to Pu are coloured and the rest of the atoms are white for clarity. Colour scheme: Fe, orange; O, red; Pu, blue; H and other atoms, white.

Table 3 Energetic, geometric and electronic structure data for the lowest energy Pu@Hem complexes: electronic and implicitly solvated reaction energy, Bader charge q, average distances, and QTAIM Pu–O bond critical point (BCP) properties (a.u.)

	Pu@Fe-term	Pu@O3-term
ΔEr/eV	−7.85	−13.90
ΔEsol/eV	−2.55	−13.53
q(Pu)	2.15	2.54
Average of closest 8 Pu–O/Å	2.65	2.55
Average of closest 4 Pu–Fe/Å	3.37	3.62
ρBCP	0.05	0.07
HBCP	0.00	−0.02
−(G/V)BCP	0.93	0.84

---

The geometries (Fig. 5) are in good agreement with the experimental EXAFS. As discussed earlier, the Fe-term is stabilised under vacuum conditions and thus was hydrated to mimic conditions found at the EARP; Pu(IV) bound directly to the Fe(III) surface seems very unlikely to occur at the EARP, and does not match the EXAFS model of eight oxygen and four iron back scatterers at 2.29 and 3.34 Å.⁴ The geometries shown in Fig. 5 highlight the closest atoms taken into account in the average distances given in Table 3. Pu(IV) forms a pentadentate complex to the hydrated Fe-term and a tetradentate complex with the O₃-term. We observe that higher coordination of oxygen results in more stable complexation (Table S11, Fig. S7 and S8), which is sensible as Pu(IV) is highly charged and may be stabilised by charge transfer to surface oxygen atoms. This was also found in our previous work studying the complexation of Pu(IV) to low-index ferrihydrite surfaces.⁶

The lack of surface water on the O₃-term reduces the structural complexity of the adsorption mode allowing for a clear description of the binding, as well as a clear charge density difference isosurface (Fig. S11). Pu(IV) binds to four surface oxygen atoms, two at shorter ∼2.18 Å and two at longer ∼2.45 Å distances, which is remarkably similar to our Pu–Fe₁₃ Keggin cluster (2.18 Å and 2.40 Å)⁵ and in good agreement with the split-shell EXAFS fitting with two backscatters at 2.22 Å and two at 2.39 Å.⁴

Our sample of complexes has a range of ΔE_r of −6.01 to −7.85 eV for Pu@Fe-term and −12.98 to −13.90 eV for Pu@O₃-term, showing that Pu(IV) is stabilised significantly by surface complexation. These reaction energies are more exothermic than those determined for a series of Fh surfaces in our previous work,⁶ which were in the range of −3 to −6 eV, suggesting that not only is Pu(IV) stabilised by complexation to Fh but that recrystallisation to hematite is further stabilising, which supports the experimental observation that Pu(IV) is retained during the transformation.⁴

Solvated reaction Gibbs energies were calculated under continuum solvation to mimic more realistic aqueous conditions. These ΔE_sol values are reduced compared to the ΔE_r, which is to be expected as the Pu fragment interacting with the continuum reduces the interaction with the surface. This reduction is significant for the non-polar hydrated Fe-term but much less so for the polar O₃-term (Table 3). It is likely that these ionic interactions are so strong in the latter case, with Pu⁴⁺ bound directly to O²⁻, that the inclusion of a solvent model “competing” with the surface interaction has much less impact.

The increased stability of metal ions bound to the O₃-term compared to the Fe-term has been observed before in the case of V(III) on Hem(001).⁵⁷ After considering a variety of adsorption geometries, as we do in our work, values of ΔE_ads ≈ −5 eV on the Fe-term and −10 eV on the O₃-term were obtained, a trend qualitatively similar to our findings for Pu(IV). The absolute values for Pu(IV) are likely more negative because of the increased charge. Also, strictly speaking, we are reporting a reaction energy, which includes a release of water molecules (eqn (4)), which is why our “binding energy” is labelled ΔE_r rather than simply ΔE_ads. Pu(IV) can compensate for any excess negative charge on the O₃-term by donating electrons to the O(2p) valence band, and in the case of V(III) it was suggested that sub-surface Fe may also receive electrons.⁵⁷

To our knowledge only one other study attempts to calculate and compare adsorption energies between a ferrihydrite surface and a hematite surface; in that case it was found that for Mn(II) ΔE_ads ≈ −4 eV when adsorbed to the Fh(001) surface and ΔE_ads ≈ −7 eV to the Hem(104) surface, in agreement with our observation that binding to hematite is more stabilising than ferrihydrite.⁵⁸ It was also found that the Bader charge on the metal ion q(Mn) was higher in the more stable configurations, which agrees with our findings regarding q(Pu) on ferrihydrite and hematite. The Bader charge for the O₃-term Pu is +2.54, which is more cationic than for the Fe-term Pu at +2.15, and the latter is closer to those determined for the Fh terminations (+2.13 to +2.15) and in all cases significantly more cationic than in [Pu(H₂O)₉]⁴⁺ (+1.65). The charge density difference (Fig. S11) shows clear charge accumulation on the surface oxygen atoms and depletion around Pu(IV), in line with the Bader charge analysis.

QTAIM analysis was performed to determine the nature of the interactions underpinning the binding of Pu(IV) to the surface. The application of QTAIM to study actinide covalency is widespread in molecular studies but essentially unknown in periodic surface studies. QTAIM molecular graphs showing the bond critical points (BCPs) are provided in Fig. S11 and the values quoted in Table 3 are the average of the eight closest oxygen atoms highlighted in Fig. 5. Further breakdowns of the QTAIM properties by surface and explicit water molecule for all the systems studied are provided in Table S12, and show that Pu–O_surf interactions are consistently stronger and more covalent than Pu–O_water, as expected. We explain the interpretation of the density ρ_BCP, total energy density H_BCP, and ratio of the kinetic to potential energy −(G/V)_BCP in detail in our previous work.⁵ Here, it suffices to say that a larger ρ_BCP, a negative H_BCP and a value of −(G/V)_BCP close to 0.5 in the range 0.5–1.0 indicate a partial covalent interaction. The interactions with the O₃-term are much stronger than with the Fe-term (Table 3) in line with the more negative reaction energies, shorter Pu–O bonds, and increased Bader charge on Pu. A comparison of all QTAIM-derived covalency metrics across all geometries, including the higher energy complexes, does not show clear trends, suggesting that the primary binding mode is ionic, as we found for Pu–Fe₁₃ Keggin and Pu@Fh surface complexation.^5,6

Fig. 6 shows how the pDOS evolves from the bare to the hydrated and Pu-complexed Fe-term. As mentioned earlier, the states located around 1.3 eV in the band gap disappear upon hydration as the coordination of the surface Fe atoms is restored. After complexation, Pu(f) states hybridise with O(p) states in the valence band around −1 eV indicating a stabilising bonding interaction. The same can be said for the O₃-term, where states at the Fermi level are greatly reduced upon complexation and strong hybridisation is observed deeper in the valence band at around −4 eV. Hybridisation between Pu(f) and Fe(d) in this region is significant compared to the Fe-term, which may be evidence of sub-surface charge transfer as discussed earlier with respect to V(III).⁵⁷ Figures of the extended pDOS down to −20 eV below the Fermi level are provided in the SI Fig. S9 and S10. We can confirm that the mixed Fe(d)–O(p) valence region extends to around −8 eV below the Fermi level and that the semicore O(s) states are located between −16 eV and −20 eV below the Fermi level. These states are all within the expected range determined by various photoemission spectroscopic methods lending confidence to the quality of the electronic structures presented.⁵⁹

Fig. 6 Projected density of states for the lowest energy Pu@Hem surface complexes. Upper; Fe-term, lower; O₃-term.

Conclusions

The complexation of Pu(IV) to two well characterised Hem(001) terminations has been studied using periodic DFT+U_eff. We benchmark bulk electronic and structural properties, finding good agreement with experiment and previous simulations including those performed with hybrid functionals, thereby justifying our methodology and once more finding that DFT+U_eff is a sensible approach. The relaxed Hem(001) terminations were then benchmarked against experimental and previous theoretical studies, verifying that our slab model is suitable to study actinide complexation.

Pu(IV) has very strong affinity to both surface terminations; the reaction energies for surface binding were found to be significantly negative, even more so than Pu(IV) bound to ferrihydrite using the same methodology. This result is significant as it supports the observations that Pu(IV) may be retained upon recrystallisation from ferrihydrite to hematite, and that Pu(IV) is strongly retained by hematite long-term. Pu(IV) binding to the O₃-term was found to be much stronger than to the Fe-term or any ferrihydrite termination studied in our previous work. This is attributed to strong ionic interactions leading to significant charge transfer and a highly cationic Pu with a Bader charge of +2.54 compared to +2.15 in the Fe-term complex. The pDOS show stronger hybridisation between the Pu(f) and O(p) and Fe(d) states in the valence band in the O₃-term compared to the Fe-term. This result agrees with previous simulations, indicating that adsorbed cationic species may transfer charge to the sub-surface Fe layer on the O₃-term.

Pu–O and Pu–Fe bond lengths were found to be similar to those determined experimentally and to those calculated for both the Pu–Fe₁₃ Keggin and Pu@Fh surface complexes, providing a structural link between the three phases. Similarly, the QTAIM analysis of the Pu–O interactions was found to be in line with those determined for the aforementioned phases, again indicating that the binding mode of Pu(IV) to these iron (oxyhydr)oxide phases is primarily ionic with partial covalency to strengthen the interaction.

Hematite is a likely candidate as a barrier to radionuclide migration in the future GDF in the UK, and we are therefore happy to report that Pu(IV) binds very strongly, particularly under oxygenated conditions. As our methodology may in principle be applied to any actinide/iron (oxyhydr)oxide mineral interaction, we hope our work stimulates further theoretical and experimental work in the field.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data supporting this article have been uploaded as part of the supplementary information (SI). Supplementary information: additional information on the convergence testing with respect to k-points and cutoff energy and the effect of U_eff. Phonon band structure and pDOS for hematite bulk and surfaces. Structural data for higher energy hydrated and Pu-complexed surfaces. Charge density plots of Pu on hematite surfaces. Coordinates of the lowest energy structures in the VASP POSCAR format. See DOI: https://doi.org/10.1039/d6dt00183a.

Acknowledgements

We are grateful for an EPSRC funded PhD studentship to R. L. D. We also thank the University of Manchester for its Computational Shared Facility (CSF3 and CSF4) and associated support services. Via our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/X035859), this work used the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk). We also extend our thanks to Dr Joseph Flitcroft at the University of Manchester for support with the Phonopy calculations.

References

1. U. Schwertmann and R. M. Cornell, Iron oxides in the laboratory: preparation and characterization, Wiley-VCH, Weinheim; New York, 2nd completely rev. and extended edn, 2000.
2. K. Sivula, F. Le Formal and M. Grätzel, ChemSusChem, 2011, 4, 432–449.
3. J. S. Weatherill, K. Morris, P. Bots, T. M. Stawski, A. Janssen, L. Abrahamsen, R. Blackham and S. Shaw, Environ. Sci. Technol., 2016, 50, 9333–9342.
4. K. F. Smith, K. Morris, G. T. W. Law, E. H. Winstanley, F. R. Livens, J. S. Weatherill, L. G. Abrahamsen-Mills, N. D. Bryan, J. F. W. Mosselmans, G. Cibin, S. Parry, R. Blackham, K. A. Law and S. Shaw, ACS Earth Space Chem., 2019, 3, 2437–2442.
5. R. L. Dempsey and N. Kaltsoyannis, Dalton Trans., 2024, 53, 5947–5956.
6. R. L. Dempsey and N. Kaltsoyannis, Environ. Sci.: Processes Impacts, 2025, 27, 2318–2328.
7. O. Stagg, K. Morris, L. T. Townsend, E. S. Ilton, L. Abrahamsen-Mills and S. Shaw, Appl. Geochem., 2023, 159, 105830.
8. C. H. Hatton, A. Christodoulidou, L. S. Natrajan and N. Kaltsoyannis, ACS Omega, 2025, 10, 17717–17726.
9. Z. D. Pozun and G. Henkelman, J. Chem. Phys., 2011, 134, 224706.
10. L. Rassouli and M. Dupuis, J. Phys. Chem. C, 2024, 128, 743–758.
11. A. Rohrbach, J. Hafner and G. Kresse, Phys. Rev. B:Condens. Matter Mater. Phys., 2004, 70, 125426.
12. N. Y. Dzade, A. Roldan and N. H. De Leeuw, Minerals, 2014, 4, 89–115.
13. E. Voloshina, in Encyclopedia of Interfacial Chemistry, Elsevier, 2018, pp. 115–121.
14. J. Goniakowski, F. Finocchi and C. Noguera, Rep. Prog. Phys., 2008, 71, 016501.
15. S. Bandaru, G. Saranya, W. W. Liu and N. J. English, Catal. Sci. Technol., 2020, 10, 1376–1384.
16. N. H. De Leeuw and T. G. Cooper, Geochim. Cosmochim. Acta, 2007, 71, 1655–1673.
17. X.-G. Wang, W. Weiss, Sh. K. Shaikhutdinov, M. Ritter, M. Petersen, F. Wagner and R. Schlögl, Phys. Rev. Lett., 1991, 81, 1038–1041.
18. T. Pabisiak and A. Kiejna, J. Chem. Phys., 2014, 141, 134707.
19. A. Kiejna and T. Pabisiak, J. Phys.: Condens. Matter, 2012, 24, 095003.
20. T. Stirner, D. Scholz and J. Sun, Surf. Sci., 2018, 671, 11–16.
21. G. Kresse and J. Hafner, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558–561.
22. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
23. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
24. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
25. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
26. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B:Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
27. P. E. Blöchl, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
28. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
29. P. E. Blöchl, O. Jepsen and O. K. Andersen, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 49, 16223–16233.
30. A. M. Ganose, A. J. Jackson and D. O. Scanlon, J. Open Source Softw., 2018, 3, 717.
31. A. Togo, J. Phys. Soc. Jpn., 2023, 92, 012001.
32. G. W. Watson, E. T. Kelsey, N. H. De Leeuw, D. J. Harris and S. C. Parker, Faraday Trans., 1996, 92, 433.
33. N. Y. Dzade and N. H. De Leeuw, Environ. Sci.: Processes Impacts, 2018, 20, 977–987.
34. N. Y. Dzade, A. Roldan and N. H. De Leeuw, Environ. Sci. Technol., 2017, 51, 3461–3470.
35. D. Santos-Carballal, A. Roldan, R. Grau-Crespo and N. H. De Leeuw, Phys. Chem. Chem. Phys., 2014, 16, 21082–21097.
36. G. Makov and M. C. Payne, Phys. Rev. B:Condens. Matter Mater. Phys., 1995, 51, 4014–4022.
37. V. Wang, N. Xu, J.-C. Liu, G. Tang and W.-T. Geng, Comput. Phys. Commun., 2021, 267, 108033.
38. B. E. Tegner, M. Molinari, A. Kerridge, S. C. Parker and N. Kaltsoyannis, J. Phys. Chem. C, 2017, 121, 1675–1682.
39. J.-L. Chen and N. Kaltsoyannis, J. Phys. Chem. C, 2019, 123, 15540–15550.
40. J.-L. Chen and N. Kaltsoyannis, J. Phys. Chem. C, 2022, 126, 11426–11435.
41. W. M. Bender and U. Becker, ACS Earth Space Chem., 2019, 3, 637–651.
42. M. S. Talla Noutack, G. Geneste, G. Jomard and M. Freyss, J. Appl. Phys., 2022, 131, 225106.
43. N. L. Banik, V. Vallet, F. Réal, R. M. Belmecheri, B. Schimmelpfennig, J. Rothe, R. Marsac, P. Lindqvist-Reis, C. Walther, M. A. Denecke and C. M. Marquardt, Dalton Trans., 2016, 45, 453–457.
44. K. Mathew, R. Sundararaman, K. Letchworth-Weaver, T. A. Arias and R. G. Hennig, J. Chem. Phys., 2014, 140, 084106.
45. K. Mathew, V. S. C. Kolluru, S. Mula, S. N. Steinmann and R. G. Hennig, J. Chem. Phys., 2019, 151, 234101.
46. A. Otero-de-la-Roza, E. R. Johnson and V. Luaña, Comput. Phys. Commun., 2014, 185, 1007–1018.
47. L. W. Finger and R. M. Hazen, J. Appl. Phys., 1980, 51, 5362–5367.
48. J. S. Olsen, C. S. G. Cousins, L. Gerward, H. Jhans and B. J. Sheldon, Phys. Scr., 1991, 43, 327–330.
49. B. Gilbert, C. Frandsen, E. R. Maxey and D. M. Sherman, Phys. Rev. B:Condens. Matter Mater. Phys., 2009, 79, 035108.
50. S. Mochizuki, Phys. Status Solidi A, 1977, 41, 591–594.
51. N. Beermann, L. Vayssieres, S.-E. Lindquist and A. Hagfeldt, J. Electrochem. Soc., 2000, 147, 2456.
52. Y. Ma, P. D. Johnson, N. Wassdahl, J. Guo, P. Skytt, J. Nordgren, S. D. Kevan, J.-E. Rubensson, T. Böske and W. Eberhardt, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 48, 2109–2111.
53. S. A. Chambers and S. I. Yi, Surf. Sci., 1999, 439, L785–L791.
54. M. Lübbe and W. Moritz, J. Phys.: Condens. Matter, 2009, 21, 134010.
55. A. Barbier, A. Stierle, N. Kasper, M.-J. Guittet and J. Jupille, Phys. Rev. B:Condens. Matter Mater. Phys., 2007, 75, 233406.
56. R. B. Wang and A. Hellman, J. Chem. Phys., 2018, 148, 094705.
57. J. Jin, X. Ma, C.-Y. Kim, D. E. Ellis and M. J. Bedzyk, Surf. Sci., 2007, 601, 3082–3098.
58. S. Tong, H. Wei, J. Zhou, Y. Yang, R. Zhu, Q. Chen, X. Xie, Q. Hu, M. F. Hochella Jr. and J. Liu, Environ. Sci. Technol., 2025, 59, 3961–3971.
59. A. Fujimori, M. Saeki, N. Kimizuka, M. Taniguchi and S. Suga, Phys. Rev. B:Condens. Matter Mater. Phys., 1986, 34, 7318–7328.

---