The defect chemistry of non-stoichiometric PuO_2±x

Abstract

An increased knowledge of the chemistry of PuO₂ is imperative for the design of procedures to store, dispose, or make use of PuO₂. In this work, point defect concentrations in PuO₂ are determined by combining density functional theory (DFT) defect energies and empirical potential calculations of vibrational entropies. The obtained defect concentrations are expressed as a function of temperature and oxygen partial pressure and used to calculate non-stoichiometry in PuO_2±x. The results show that the defect chemistry of PuO₂ is dominated by oxygen vacancies and interstitials. Hypo-stoichiometric PuO_2−x is accommodated by both the uncharged oxygen vacancy and positively charged oxygen vacancy at small values of x, with increasingly dominant with increasing x. The negatively charged oxygen interstitial (O²⁻_i) is found to accommodate hyper-stoichiometry (PuO_2+x), but reluctance to form hyper-stoichiometric PuO_2+x is observed, with oxygen interstitials present only in very low concentrations irrespective of conditions. The small degree of hyper-stoichiometry found is favoured by low temperatures.

---

1 Introduction

The handling, storage and disposal of PuO₂ (the form civilian plutonium stockpiles are maintained in)¹ presents considerable safety and security challenges. Large stockpiles of PuO₂ currently exist in interim storage awaiting decisions regarding its long-term disposal or re-use; PuO₂ attracts interest as a component of future nuclear reactor fuels and ²³⁸Pu, in the form ²³⁸PuO₂, is the most commonly used isotope in radioisotope thermoelectric generators and heating units for space applications.² To enable safe storage and/or application a thorough understanding of how the material evolves in the short to medium term is required, which can be aided by studying the materials defect chemistry.

The potential pressurisation of some of the containers used to store PuO₂ is of major concern, with multiple pressurisation mechanisms proposed.³ The formation of hyper-stoichiometric PuO_2+x may play a role in the chemical reactions producing gas. It is reported that PuO_2+x participates in moisture enhanced corrosion of Pu, producing H₂,^4–6 although this is less of a concern in the UK, where the oxide is stored. Through experimental investigation of PuO₂ at conditions relevant to intermediate storage, Haschke et al.^4,5 have reported the existence of the hyper-stoichiometric PuO_2+x phases, while subsequent work has shown PuO_2+x may be best described as a compound containing OH⁻ and H₂O species, with excess oxygen accommodated in molecular complexes of Pu(V).^7,8 The defect chemistry of PuO₂ helps describe the formation and properties of non-stoichiometric PuO₂ and the equilibrium of species (e.g. oxygen) between the PuO₂ solid and the surrounding gases.

PuO₂ is very challenging to study both experimentally and theoretically. Experimental challenges include strict regulation, the difficulty in preparing and handling samples, and continuous crystal lattice damage and transmutation of Pu due to self-irradiation.⁹ Many of the challenges involved in theoretical calculation arise from trying to correctly describe the 5f electrons. PuO₂ is classified as a strongly correlated material, with a ground state of localised 5f electrons.¹⁰ In conventional density functional theory (DFT), an over de-localisation of these electrons results in the self-interaction error that causes PuO₂ to be described as conducting when in fact PuO₂ is classified as a charge-transfer insulator.¹⁰

A method used extensively to overcome this shortcoming is to adopt a Hubbard-like correction to give the DFT+U method.¹⁰ The U parameter modifies the electron repulsion of the 5f electrons and addresses the self-interaction error, yielding good reproduction of the ground-state electronic properties seen in actinide dioxides.^10,11 Two adjustable variables help construct the U correction parameter. These are the on-site Coulombic (U) and exchange interaction (J) parameters, which are used to calculate U_eff = U − J in the Dudarev et al. formulism.¹² The U and J parameters are obtained by fitting to structural and electronic properties of PuO₂.^11,13–17 Reproduction of the electronic band gap (E_g) is seen as a strong metric of success. Most previous DFT+U studies of PuO₂ fit U and J to reproduce the bandgap of McNeilly et al.¹⁸ (E_g = 1.8 eV), however the much more recent experimental study of McCleskey et al.¹⁹ reports a higher E_g of 2.8 eV. This large discrepancy in the experimental bandgap for PuO₂ poses significant issues for determination of the appropriate values for U and J.

An alternative method to DFT+U used to avoid the problems encountered with the self-interaction error is to use hybrid functionals. Hybrid functionals, which blend a portion of the Hartree–Fock (HF) exchange into a part of a density functional, have been shown to offer significantly improved descriptions of band gaps, especially in small- to medium- gap systems (∼<5 eV).²⁰ Compared to the DFT+U method, hybrid functionals require significantly more computational effort to utilise, but several studies have applied hybrid functionals to study PuO₂ with good experimental replication achieved.^14,21

An important aspect in correctly simulating PuO₂ is selection of an accurate description of its magnetic ground state. Experimental and theoretical studies have indicated differing ground states: a singlet diamagnetic (DM) ground-state is indicated from experimental observations²² whereas ferromagnetic (FM) or anti-ferromagnetic (AFM) ground states have in the past been predicted from theory.^14,17,21,23 Many studies select 1k AFM order to describe PuO₂.^{13,17,21,24–26} In their comprehensive study of magnetic order in PuO₂, Pegg et al.¹⁴ found that 1k AFM states produce an incorrect crystal structure (the correct structure as indicated by experimental results is Fm3̄m crystal symmetry)²⁷ and that the choice of AFM domain (longitudinal or transverse) had a significant impact on the electronic and structural properties of PuO₂. This is a key finding, as most studies do not differentiate between collinear magnetic structures. Pegg et al.¹⁴ propose a longitudinal 3k AFM ground-state. Of all magnetic configurations tested, this was the lowest in energy, was able to retain Fm3̄m crystal symmetry and replicated well the physical properties of PuO₂.

To obtain this ground-state, spin orbit interaction (SOI) was included. Not including SOI resulted in a different magnetic ground state being obtained, highlighting the importance of its inclusion. The different observations found in experiment are perhaps a consequence of the aforementioned experimental challenges. Possible theories include that current experimental resolution cannot identify an ordered magnetic state and/or an AFM-DM transition could be occurring outside the temperature ranges investigated.¹⁴ The experimental measurement of isostructural NpO₂ has also proven difficult, where an odd number of electrons should result in a magnetic moment.

Several studies have investigated the formation energies of intrinsic and extrinsic defects in PuO₂.^{24,26,28–33} Thermochemical models have been produced that study the defect chemistry of PuO_2−x, over a large stoichiometry range,^34,35 whilst using DFT+U and a point defect model, Lu et al.²⁴ determined the stability of charged defects within PuO_2±x. Lu et al.²⁴ find that, at 1000 K, oxygen vacancies dominate in the region of hypo-stoichiometry, whilst oxygen interstitials dominate the region of hyper-stoichiometry.

It is common when creating a point defect model from first principles to assume that the Gibbs free energy to form a defect can be well approximated by the defect formation energy. However, this approximation neglects the contribution of the vibrational entropy, S_vib. Recent work by Cooper et al.³⁶ and Soulie et al.³⁷ on UO₂ has shown that it is only by including the change in vibrational entropy due to defect incorporation that it is possible to reproduce the defect chemistry observed in experiment.

Therefore, the objective of this paper is to understand the defect chemistry of PuO₂ and how it might evolve under storage conditions taking into account the complex electronic and magnetic structures and incorporating the vibrational entropy into an updated point defect model.

2 Method

To determine how the concentrations of defects in PuO₂ varies with environmental conditions requires the combination of DFT defect energies and simple thermodynamics, starting from eqn (1). Using Boltzmann statistics, this equation shows that the concentration of defect i, c_i, can be related to the change in Gibbs free energy, ΔGⁱ_f, to form defect i:where m_i is the multiplicity of equivalent sites, k_B is the Boltzmann constant and T is the temperature. ΔGⁱ_f is defined as:where ΔE is the defect energy: the difference in the total DFT-derived energies of the defective and perfect supercells. ΔS_vib is the difference in vibrational entropy between defective and perfect supercells. Within the third term, n_α is the number of atoms of species, α, added or removed from the system to make defect i, where μ_α is the chemical potential of species α. q_i is the charge of the defect i, and μ_e is the electron chemical potential. Finally, E_corr is a correction term necessary to correct for finite-size effects arising due to the use of relatively small simulation supercells.

Following Cooper et al.³⁶ and Soulie et al.,³⁷ two differing atomistic simulation approaches were applied to obtain the parameters of eqn (2): ΔE is found using DFT and ΔS_vib is obtained using empirical pair potentials. The vibrational entropy of a lattice is a function of the lattices phonon frequencies which are determined through force calculations associated with the displacements of atoms from their ground state positions. This becomes a very large calculation when defects are introduced, due to the removal of symmetry. Therefore, empirical potentials are used to calculate the vibrational entropies. The calculation and use of all the terms that contribute to eqn (2) are discussed in more detail throughout the rest of this section.

2.1 Defect energies

DFT calculations were performed using the Vienna ab initio Simulation Package (VASP).^38–41 The projector augmented wave (PAW)^42,43 method has been implemented with the frozen-core approximation. 5f⁴6s²6p⁶6d²7s² of Pu and 1s²2s²2p⁶ of O are treated as valence electrons. Following convergence testing, the cut-off energy for the planewave basis set was selected to be 500 eV and a 4 × 4 × 4 Monkhurst-Pack k-point mesh⁴⁴ was used for the 12-atom PuO₂ unit cells. For 2 × 2 × 2 supercells used in defect simulations, the k-point mesh was adjusted to 2 × 2 × 2. SOI⁴⁵ was considered in the calculations with a longitudinal 3k AFM magnetic configuration adopted for PuO₂, described in Fig. 1.

Fig. 1 The crystal (left) and magnetic (centre) structure of Fm3̄m PuO₂, with plutonium and oxygen atoms depicted with blue and red spheres respectively. The direction of the effective magnetic moment (μ) on each plutonium ion is shown by an arrow and the μ = (a,b,c) wave-vector matrix (right).

To study the properties of the bulk PuO₂, we replicate the study of Pegg et al.¹⁴ and use the hybrid Heyd–Scuseria–Ernzerhof (HSE06) functional:^46–49

where E^HSE_XC is the exchange–correlation HSE06 energy, E^HF,SR_X is the short-range HF exchange, E^PBE,SR_X and E^PBE,LR_X are the short-range and long-range components of the Perdew–Burke–Ernzerhof (PBE) exchange functional respectively and E^PBE_C is the PBE correlation energy. A quarter of the short-range HF exchange is used in this functional and the screening parameter (μ) is 0.207 Å⁻¹. The iteration threshold for electronic and ionic convergence when using this functional is 1 × 10⁻⁶ eV and 1 × 10⁻² eV Å⁻¹, respectively.

For calculation of defect energies in PuO₂ supercells hybrid functionals are computationally too expensive; DFT calculations are therefore performed with the generalised gradient approximation (GGA) using the PBE functional revised for solids (PBEsol).^50,51 The PBEsol functional was shown by Pegg et al.¹¹ to best reproduce experimental properties in actinide dioxides after the HSE06 functional. The strong correlations due to f-electrons in Pu were accounted for by applying the DFT+U method using the Liechtenstein formulism.⁵² When using PBEsol+U, the iteration threshold for electronic convergence is 1 × 10⁻⁶ eV and for ionic convergence is 1 × 10⁻² eV Å⁻¹ and 2 × 10⁻² eV Å⁻¹ when simulating PuO₂ unit cells and defect containing PuO₂ supercells, respectively. In this study, the U parameter within PBEsol+U was selected such that it reproduced the band gap obtained from the HSE06 functional (3.04 eV). It was chosen to reproduce the HSE06 bandgap as the experimental data shows a large variation and this functional has been proven to replicate experimental bandgaps.²⁰ This resulted in a U parameter of 7.0 eV being used for this study (see Fig. 2). The use of a high U parameter is not without precedence, with several previous studies adopting U values exceeding 6.0 eV when studying actinide dioxides, including PuO₂.^11,53 The J parameter was fixed at a value of 0.0 eV throughout this study, as any introduction of J was shown by Pegg et al.¹¹ to detrimentally affect the reproduction of the band gap for PuO₂. The resulting equilibrium properties found for longitudinal 3k AFM PuO₂ simulated using either the HSE06 functional or PBEsol+U (U = 7.0 eV) are presented in Table 1. The discrepancy between the theoretical and experimental description of the magnetic properties of PuO₂ is discussed in the introduction.

Fig. 2 Variation of band gap as a function of the Coulomb modifier (U). Band gap value calculated using the HSE06 functional shown with horizontal dotted line.

Table 1 The lattice constant (Å), magnetic moment (μ_B/Pu ion), band gap energy (eV) and crystal symmetry for longitudinal 3k AFM PuO₂ calculated using either the HSE06 functional or PBEsol+U (U = 7.0 eV). In addition, we show reported experimental values

		Lattice constant (Å)	Magnetic moment (μB/Pu ion)	Band gap (eV)	Crystal symmetry	Magnetic configuration
This work	HSE06	5.379	3.65	3.04	Fm3̄m	Longitudinal 3k AFM
	PBEsol+U	5.415	3.82	3.03	Fm3̄m	Longitudinal 3k AFM
Experiment		5.3995	0.0022	1.80^18/2.80^19	Fm3̄m^27	DM^22

---

Defects were inserted into 2 × 2 × 2 expansions of the 12-atom PuO₂ unit cell. Only one unique site exists for the vacancies and interstitials of plutonium and oxygen considered, due to the symmetry of the PuO₂ lattice. Different charge states were considered for each defect, modelled by adding or removing electrons from the supercell. Defect-containing supercells were relaxed under constant volume, using the lattice constants obtained from defect-free simulations. The defects considered in this study are denoted below, using Kröger–Vink notation, modified to display charge as an integer value (no charge indicated by: ×):⁵⁴

Oxygen interstitials: O_i^×, O¹⁻_i and O²⁻_i

Oxygen vacancies: , V¹⁺_O and V²⁺_O

Plutonium interstitials: , Pu¹⁺_i, Pu²⁺_i, Pu³⁺_i and Pu⁴⁺_i

Plutonium vacancies: , V¹⁻_Pu, V²⁻_Pu, V³⁻_Pu and V⁴⁻_Pu

2.2 Finite size correction

The introduction of charged defects into the small simulation supercells accessible using DFT introduces a number of finite size effects as discussed extensively in ref. 55. These include Coulombic interactions between the defect and its periodic image as well as with the neutralising background charge. The result is that defect formation energies exhibit a strong dependence on the size of the supercell used and this must be corrected for. The most basic correction that can be applied to remove these interactions is the point charge (PC) correction. For a cubic system, such as PuO₂, the PC correction energy can be written as:⁵⁶where α = 2.837 is the lattice-type dependent Madelung constant, L is the supper cell lattice constant and ε is the static dielectric constant. The static dielectric constant of PuO₂ was calculated using density functional perturbation theory (DFPT)^57,58 as implemented in VASP. We use the previously described longitudinal 3k AFM PuO₂ unit cell with PBEsol+U (U = 7.0 eV). The static dielectric tensor of PuO₂, ε_ij is calculated by summing its electronic and ionic components. The average of the static dielectric tensor, is calculated and used as the static dielectric constant of PuO₂. We obtain a static dielectric constant of 19.66. This compares favourably to the experimental value (18.9⁵⁹) and falls between the experimental value for UO₂ (23.8⁶⁰) and the DFT-obtained value for AmO₂ (17.40⁶¹), thus fits within the emerging trend.

Makov and Payne (MP) extended the PC correction to include a term with L⁻³ order.⁶² More recently Freysoldt, Neugebauer and Van de Walle (FNV) developed a correction that compares the planar-averaged electrostatic potentials of supercells with and without defects (V_defect,q & V_bulk respectively).⁶³ The FNV scheme correction energy is summarised following ref. 64 as:

E_corr = E^iso_PC − qΔV_PC,q/b|_far (5)

V_q/b = V_defect,q − V_bulk (6)

ΔV_PC,q/b = V_q/b − V_PC,q (7)

ΔV_PC,q/b is the potential difference between the defect induced potential (V_q/b) and the PC potential V_PC,q.⁶³ ΔV_PC,q/b|_far is ΔV_PC,q/b at a position far from the defect site but still within the supercell. The FNV scheme has been demonstrated to be an effective correction in many systems,⁶⁵ however the use of planar-averaged electrostatic potentials can provide unreliable corrections when the atomic positions of the defective supercell are allowed to relax. This was observed in the system investigated in this study: with the example of a V⁴⁻_Pu supercell, V_PC,q/b|_far could not be determined accurately as V_q/b was seen to strongly fluctuate whilst V_PC,q remains parabolic (illustration in appendix). Consequently, this study employs the correction of Kumagai and Oba which uses atomic site electronic potentials as opposed to planar averaged electrostatic potentials.⁶⁴ Using the same example of a V⁴⁻_Pu supercell, Fig. 3 shows that this scheme is able to produce ΔV_PC,q/b|_far that is close to constant across the sampling region. Kumagai and Oba also demonstrate that a potential alignment correction, sometimes adopted to shift the valence band maximum of a charged defect to that of a pristine host, is not required.⁶⁴ The difference in ΔE between a V⁴⁻_Pu 3 × 3 × 3 supercell and a V⁴⁻_Pu 2 × 2 × 2 supercell was found to decrease from 1.29 eV when no correction was applied to 0.06 eV after correction following the Kumagai and Oba scheme. This gives confidence that the correction works well for the system studied here.

Fig. 3 V _q/b, V_PC,q and ΔV_PC,q/b calculated at the atomic positions of a PuO₂ supercell with a V⁴⁻_Pu defect. ΔV_PC,q/b|_far calculated by averaging ΔV_PC,q/b across sampling region.

2.3 Vibrational entropies. Vibrational entropies are obtained using empirical potentials, with the General Utility Lattice Program (GULP)⁶⁶ together with the Cooper, Rushton and Grimes (CRG)^33,67 potential. This is a many-body potential model used to describe actinide oxide systems which achieves good reproduction of thermodynamic and mechanical properties. The phonon density of states (DOS) for PuO₂, calculated with the CRG potential and compared to the experimental data of Manley et al.⁶⁸ is shown in Fig. 4. The calculated phonon DOS compares reasonably well to the experimental data. Many of the characteristic peaks and troughs are replicated, although at high frequency the magnitude of these features is smaller for the calculated data. This is a result of adopting a large k-point mesh that smooths the phonon DOS, and at high frequencies the experimental data has wide error bars.⁶⁸ The calculation of vibrational entropies, S_vib, closely follows the approach described in ref. 36, 69 and 70, where using eqn (8), defect entropies are calculated from the normal vibrational frequencies, v_n, which are themselves calculated by diagonalizing the dynamical matrix of the system.

Fig. 4 The phonon density of states of PuO₂ calculated with the CRG empirical potential (blue), compared to experimental data⁶⁸ (red).

In this formula, h is Plancks's constant and N is the number of atoms in the crystal. In this study, the system to calculate vibrational entropies is a 4 × 4 × 4 expansion of the PuO₂ unit cell. Defective supercells are created by adding or removing atoms and are relaxed under constant volume in the same way as in the DFT simulations. Defect vibrational entropies are found by calculating the difference in vibrational entropies between the defective and perfect supercell (ΔS_vib) and presented in Table 2. The same value of ΔS_vib is given to all charge states of a given defect.

Table 2 The difference in vibrational entropy of PuO₂ due to the addition of a defect, calculated using the CRG potential

T (K)	Defect entropy, ΔSvib/kB
	VO	Oi	VPu	Pui
400	0.545	4.572	−3.783	2.518
600	−0.453	5.373	−5.245	3.110
800	−1.230	6.058	−6.208	3.725
1000	−1.868	6.638	−6.940	4.270
1200	−2.391	7.137	−7.508	4.746
1400	−2.843	7.578	−7.996	5.164
1600	−3.238	7.961	−8.402	5.547
1800	−3.574	8.309	−8.761	5.883
2000	−3.899	8.599	−9.086	6.174

---

2.4 Chemical potentials

The starting point for determining the chemical potentials is to define the chemical potential of solid PuO₂ (μ_PuO2(s)) in terms of the chemical potential per formula unit of the constituent species, namely plutonium, μ_Pu(P_O2,T), and oxygen, μ_O2(P_O2,T):

μ_Pu(P_O2,T) + μ_O2(P_O2,T) = μ_PuO2(s) (9)

For a solid, , therefore the temperature and pressure dependencies have been dropped. In equilibrium conditions, the chemical potential of Pu cannot exceed that of solid Pu, otherwise a Pu precipitate would form. This upper bound is the Gibbs free energy of Pu in its natural state. It can therefore be said that in Pu rich conditions:

μ_Pu(P_O2,T) = μ_Pu(s) (10)

To find μ_Pu(s) we simulate the α phase of Pu with DFT using PBEsol+U. We encompass the recommendation of the review by Söderlind et al.⁷¹ to use small U and J values, setting U and J parameters at 2.2 eV and 0.58 eV respectively. The atomic volume obtained with these values (18.27 ų) matched closely the atomic volume obtained by Söderlind et al.⁷¹ when using PBE+U. To determine the chemical potential of oxygen the approach of Finnis et al.⁷² is adopted. This method uses the known experimental formation energy of the oxide ⁷³ to obtain the chemical potential of oxygen at standard temperature and pressure:

Unlike the solid species in eqn (11), the temperature and pressure dependence of the oxygen chemical potential cannot be neglected and is extrapolated from using formulas constructed by Johnston et al.:⁷⁴

where the temperature dependent Gibbs free energy per mole is fitted to a polynomial derived from experimental data (coefficients listed in Table 3):⁷⁴

Table 3 Coefficients for Gibbs free energy expression in eqn (13)⁷⁴

A	29.659 × 10^−3 kJ mol^−1 K^−1
B	6.137261 × 10^−6 kJ mol^−1 K^−1
C	−1.186521 × 10^−9 kJ mol^−1 K^−1
D	0.095780 × 10^−12 kJ mol^−1 K^−1
E	−0.219663 × 10^−3 kJ mol^−1 K^−1
F	−9.861391 kJ mol^−1 K^−1
G	237.948 × 10^−3 kJ mol^−1 K^−1

---

The electron chemical potential, μ_e = E_VBM + ε_F, is expressed as the sum of the energy of the valence band maximum (VBM), E_VBM, and the electron chemical potential above the VBM, ε_F. As overall charge neutrality of the system must be maintained, the concentrations of ionic and electronic defects must be such that at any given temperature and oxygen partial pressure the following criteria is met:⁷⁵

The first term is the sum of the charges the of ionic defects. The second and third terms are given by applying Fermi–Dirac statistics to the electronic DOS to obtain the concentrations of electrons (e⁻) in the conduction band and concentration of holes (p⁻) in the valence band, respectively. Within these two integrals are g_v(E) and g_c(E),the density of electronic states in the valence band and conduction band per formula unit of PuO₂, respectively. For calculation of the electron population, E_CBM is the energy of the conduction band minimum.

With the above expressions it is possible to construct diagrams showing the variation of defect concentrations with environmental conditions; at each condition there exists a unique value of ε_F that satisfies the charge neutrality condition. The Defect Analysis Package⁷⁶ employs a linear bisection to find the value of ε_F that ensures charge neutrality for any given oxygen partial pressure and temperature. Additionally, the deviation from stoichiometry, x in PuO_2±x, is calculated from the defect concentrations summed over all charge states:

3 Results & discussions

The formation energy of all defects considered are plotted as a function of the Fermi level in Fig. 5. Fig. 5 shows the formation energies at 1000 K and O-rich (10⁻² atm)/O-poor (10⁻³⁰ atm) conditions, displaying only the charge state of each defect that corresponds to the lowest formation energy at a given position in the bandgap. Under O-rich conditions Pu_i defects have very high formation energy and so are not shown in Fig. 5. For each defect it can be seen that, to varying degrees, the charge state with the lowest energy varies across the bandgap. V_Pu defects have the least variation, its full ionic charge state (V⁴⁻_Pu) has the lowest energy across the majority of the bandgap. Conversely, under halfway across the bandgap, the neutral defect is seen to overtake the fully charged V²⁺_O defect as the oxygen vacancy with the lowest formation energy. Prodan et al.²¹ report that O¹⁻_i is the most energetically favourable charge state for the oxygen interstitial. We find that although the O²⁻_i defect has the lowest energy for much of the bandgap, O¹⁻_i becomes favourable at low Fermi level values, with O^×_i most favourable at the lowest Fermi levels close to the valence band. Therefore, under certain conditions we also find O¹⁻_i to be the preferred interstitial. Across the band gap the favourable charge state of Pu_i defects are seen to vary, with each charge state except Pu²⁺_i possessing the lowest formation energy over fairly large ranges in Fermi level.

Fig. 5 Defect formation energies for O vacancies & interstitials as well as Pu vacancies & interstitials as a function of Fermi energy. Calculated at 1000 K and oxygen partial pressure of 10⁻² atm (left) and 10⁻³⁰ atm (right). Only the charge state with the lowest formation energy for a given Fermi level is shown for each defect, represented with numeric label. At an oxygen partial pressure of 10⁻² atm Pu interstitials are omitted from the plot as their formation energy is very high.

Table 4 displays the calculated reaction energy of unbound oxygen and plutonium Frenkel pairs (FP) and the unbound Schottky tri-vacancy (STV), with defects formally charged. Each reaction energy is found to be within the range of previous theoretical calculations, which are provided for comparison in Table 4. The ranges are seen to be very wide; in the DFT+U studies, the U and J parameters adopted varied significantly indicating the importance of electron repulsion of the 5f electrons on the formation of point defects. The review of Murch et al.⁷⁷ provides the only experimental comparison available, indicating that the oxygen Frenkel pair has a much lower formation energy than calculated here. Caution is however attached to the oxygen Frenkel pair experimental result, which is described by Murch et al.⁷⁷ as “on the low side”. Our results match the established pattern in all techniques, where O-FP < STV < Pu-FP, showing that oxygen-type defects are significantly more favourable than plutonium-type defects.

Table 4 Reaction energies for the intrinsic processes in PuO₂ with comparison to experimental and previous theoretical studies. Theoretical studies are separated further, based on if results are obtained using DFT or empirically

Process	Charge states	Reaction energy (eV)
		DFT+U					Empirical pair potential			Experiment
		This work	LSDA+U^24	LDA+U^30	GGA+U^31	GGA+U^32	Ref. 26	Ref. 27	Ref. 32	O diffusion^77
O Frenkel pair	V^2+O + O^2−i	4.89	3.48	4.58	9.78	3.90	5.52	2.66	4.90	2.72–2.92
Pu Frenkel pair	V^4−Pu + Pu^4+i	13.86	15.19	10.02	19.78	11.90	13.58	10.03	24.00
Schottky defect	V^4−Pu + 2V^2+O	7.65	7.51	6.09	14.92		10.40	3.54

---

The resulting point defect concentrations, calculated as a function of temperature at a fixed oxygen partial pressure of 10⁻⁵ atm, are presented in Fig. 6. Concentrations of oxygen defects are observed to be multiple orders of magnitude greater than cation defect concentrations across all temperatures. This agrees with previous experimental work that metal defects are minority defects in PuO₂.⁷⁷ All plutonium vacancy and interstitial defects are predicted to have such low concentrations that do not appear in Fig. 6.

Fig. 6 The defect chemistry of PuO_2±x as a function of temperature at an oxygen partial pressure of 10⁻⁵ atm. The diagrams include (top) and exclude (bottom) vibrational entropy contributions when calculating defect concentrations.

Fig. 6 also displays the impact of including vibrational entropy contributions. Inclusion of vibrational entropy is seen to increase the temperature at which perfect stoichiometry occurs: ≈1200 K when vibrational entropy is omitted versus ≈1425 K when included. Perfect stoichiometry occurs when [O_i] = [V_O], meaning that vibrational entropy inclusion further promotes O_i defects. From the point of perfect stoichiometry, increasing or decreasing temperature results in hypo-stoichiometric or hyper-stoichiometric PuO_2±x, respectively. The diagrams display the reluctance for PuO₂ to be hyper-stoichiometric in oxygen, with x in PuO_2+x peaking at approx. 10⁻⁷ at 1210 K (when vibrational entropy included), representing negligible hyper-stoichiometry. As temperature decreases from this peak, the concentration of defects and the value of x in PuO_2+x falls. The concentration of defects becomes increasingly negligible; this continues beyond the range of Fig. 6, at lower temperatures. It can be seen hypo-stoichiometry is much more favourable, which correlates with the wide hypo-stoichiometric region evident in the phase diagram of PuO₂.⁷⁸

From this point onwards, all plots include vibrational entropy. Broadly speaking, the diagram in Fig. 6 can be broken down into three regions based upon electron and hole concentrations: these are a hole-dominant region, a region of equal concentration and an electron dominant region. At low temperatures the O_i defects dominate and are charge compensated by holes. The [O_i] : [p⁺] ratio is approximately 1 : 2, due to the dominance of the O²⁻_i defect. Increasing the temperature past the point of peak-hyper-stoichiometry, which coincides with peak [O_i], the hole concentration is no longer set such that it provides charge compensation to O_i defects. Instead, the concentration of holes matches the concentration of electrons and the Fermi level remains fixed at the midpoint of the bandgap. [V_O] increases with temperature, resulting in increasing hypo-stoichiometry. At high temperatures, [V_O] reaches levels close to the concentration of holes and electrons. To accommodate any further increases in [V_O], the concentration of electrons increases to maintain charge neutrality. The concentration of holes then decreases as the temperature increases resulting in the Fermi level increasing. In this region, the [V_O] : [e⁻] ratio varies in order to provide charge compensation to both V²⁺_O and V¹⁺_O defects.

The defect chemistry of PuO₂ as a function of oxygen partial pressure is presented in the Brouwer diagrams of Fig. 7, at temperatures 1000 K and 2000 K. At 1000 K, in the region close to perfect stoichiometry, where x in PuO_2±x is at its smallest, O²⁻_i is the dominant oxygen interstitial. The dominant oxygen vacancies are found to be V²⁺_O and , which have concentrations of similar magnitude, with at narrowly higher concentrations. In this region of near-stoichiometry, the defect chemistry and dependence on oxygen partial pressure can be attributed to the following defect reactions occurring:

Fig. 7 Calculated Brouwer diagram showing the defect chemistry of PuO_2±x as a function of the oxygen partial pressure at 1000 K (top) and 2000 K (bottom). Both diagrams include vibrational entropy contributions when calculating defect concentrations.

Application of the law of mass action to eqn (16), (17) and (18) finds equilibrium constants of k₁, k₂ and k₃, respectfully:

Given that [e⁻] = [p⁺], it can be seen that [O²⁻_i] is proportional to , whilst [V²⁺_O] and are proportional to in this region. Increasing the oxygen partial pressure beyond the region of near-stoichiometry sees the concentration holes, O²⁻_i defects and x in PuO_2+x transition to an oxygen partial pressure dependence of . With too few electrons available, two holes are created for charge compensation of the O²⁻_i defect. Eqn (22) and (23) show the defect reaction responsible and the corresponding equilibrium constant, respectively. Given that [O²⁻_i] = 2[p⁺], it can be shown using eqn (23) that [O²⁻_i] is proportional to in the hyper-stoichiometric region.

Reducing the oxygen partial pressure below the region of near-stoichiometry sees the concentration of V²⁺_O defects transition to a dependence on the oxygen partial pressure of , whilst concentration of defects maintain a dependence on the oxygen partial pressure of . Consequently, it is eqn (18) that is the dominant defect reaction in the hypo-stoichiometric region, causing defects to become the dominant defect and the source of hypo-stochiometric PuO_2−x. The dependence of x in PuO_2−x has been shown experimentally to have a dependence on the oxygen partial pressure of ,³⁵ indicative of the V¹⁺_O vacancy being dominant. This result reflects experimental data at large non-stoichiometry, where the point defect model breaks down due to the occurrence of complex processes involving defect clusters. Our result, that reflects very small non-stoichiometry, can therefore not be directly compared.

Comparing the two Brouwer diagrams in Fig. 7, increasing temperature creates a more reducing environment, increasing oxygen vacancy concentrations. Perfect stoichiometry therefore occurs at higher oxygen partial pressures. Aside from this, in the 1000 K temperature range investigated here, temperature is found not to alter the defect chemistry of PuO₂ significantly.

The finding that the defect chemistry of PuO₂ is dominated by oxygen interstitials and vacancies agrees with the predictions in the PuO₂ point defect study of Lu et al.²⁴ However, the present study has found a significant increase in the concentrations of non-formally charged defects relative to their formally charged counterparts. Whereas Lu et al.²⁴ found [O²⁻_i] and [V²⁺_O] defects to be multiple orders of magnitude greater than any non-formally charged defect, here is the dominant oxygen vacancy and concentration of O¹⁻_i is predicted to be much closer to that of O²⁻_i. This can be partially attributed to the larger bandgap used in this work (3.04 eV vs. 1.70 eV). Electrons and holes, which are required to provide charge compensation to charged defects, have an increased formation energy with increased band gap. Therefore, increasing the bandgap increases the favourability of defects with a smaller charge magnitude. Additionally, we report a much lower degree of hyper-stoichiometry. For example, at a temperature of 1000 K and oxygen partial pressure of 10⁻⁵ atm Lu et al.²⁴ report a value of x in PuO_2+x of ≈10⁻³, several orders of magnitude greater than the value of ≈10⁻⁶ found in this study.

4 Conclusions

In summary, using DFT+U with a longitudinal 3k AFM magnetic structure and SOI, the point defect concentrations in PuO_2±x as a function of temperature and oxygen partial pressure have been determined. The results predict that the defect chemistry of PuO₂ is dominated by oxygen vacancies and interstitials, but that even at high oxygen partial pressures, PuO₂ is found reluctant to form hyper-stoichiometric PuO_2+x, with oxygen interstitials present only in very low concentrations irrespective of conditions. Hypo-stoichiometric PuO_2−x is predicted, in agreement with the PuO₂ phase diagram,⁷⁸ and is accommodated by both the and V²⁺_O defects at small values of x, with increasingly dominant with increasing x. The little hyper-stoichiometry is accommodated by the O²⁻_i defect. Aside from increasing the temperature at which perfect stoichiometry occurs, inclusion of vibrational entropy into the defect formation energy calculation was not found to alter which defects accommodate non-stoichiometry in PuO_2±x, as observed in UO₂.^36,37

Conflicts of interest

There are no conflicts to declare.

Appendix

See Fig. 8.

Fig. 8 Planar averaged V_q/b, V_PC,q and ΔV_PC,q/b calculated at the atomic positions of a PuO₂ supercell with a V_Pu⁴⁻defect.

Acknowledgements

We thank D. O. Scanlon and M. Pacey for their valuable help. This work made use of the High End Computing facility at Lancaster University. The visualisation of crystal structures was performed with VESTA,⁷⁹ finite-size corrections were calculated using pymatgen.⁸⁰ This project was funded as part of EPSRCs TRANSCEND project (EP/S01019X/1).

References

1. J. M. Haschke and J. L. Stakebake, Handling, Storage, and Disposition of Plutonium and Uranium, in The Chemistry of the Actinide and Transactinide Elements, ed. L. R.Morss, N. M. Edelstein and J. Fuger, Springer, Dordrecht, 3rd edn, 2010, vol. 5, ch. 29, pp. 3199–3272..
2. R. C. O’Brien, R. M. Ambrosi, N. P. Bannister, S. D. Howe and H. V. Atkinson, J. Nucl. Mater., 2008, 377, 506–521.
3. G. Bailey, E. Bluhm, J. Lyman, R. Mason, M. Paffet, G. Polansky, G. D. Roberson, M. Sherman, K. Veirs and L. Worl, Gas Generation from Actinide Oxide Materials, LA-13781-MS, Los Alamos National Laboratory, Los Alamos, NM, 2000.
4. J. M. Haschke, T. H. Allen and L. A. Morales, Science, 2000, 287, 285–287.
5. J. M. Haschke and T. H. Allen, J. Alloys Compd., 2002, 336, 124–131.
6. J. M. Haschke, T. H. Allen and L. A. Morales, J. Alloys Compd., 2001, 314, 78–91.
7. J. D. Farr, R. K. Schulze and M. P. Neu, J. Nucl. Mater., 2004, 328, 124–136.
8. S. T. Conradson, B. D. Begg, D. L. Clark, C. den Auwer, M. Ding, P. K. Dorhout, F. J. Espinosa-Faller, P. L. Gordan, R. G. Haire, N. J. Hess, R. F. Hess, D. Webster Keogh, L. A. Morales, M. P. Neu, P. Paviet-Hartmann, W. Runde, C. Drew Tait, D. Kirk Veirs and P. M. Villella, J. Am. Chem. Soc., 2004, 126, 13443–13458.
9. Los Alamos Science, ed. N. G. Cooper, Los Alamos National Laboratory, Los Alamos, NM, 2000, vol. 26.
10. X.-D. Wen, R. L. Martin, T. M. Henderson and G. E. Scuseria, Chem. Rev., 2013, 113, 1063–1096.
11. J. T. Pegg, X. Aparicio-Anglès, M. Storr and N. H. de Leeuw, J. Nucl. Mater., 2017, 492, 269–278.
12. 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.
13. P. S. Ghosh and A. Arya, Phys. Chem. Chem. Phys., 2019, 21, 16818–16829.
14. J. T. Pegg, A. E. Shields, M. T. Storr, A. S. Wills, D. O. Scanlon and N. H. de Leeuw, Phys. Chem. Chem. Phys., 2018, 20, 20943–20951.
15. H. Nakamura, M. Machida and M. Kato, Prog. Nucl. Sci. Technol, 2011, 2, 16–19.
16. P. Zhang, B.-T. Wang and X.-G. Zhao, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 144110.
17. B. Sun, P. Zhang and X.-G. Zhao, J. Chem. Phys., 2008, 128, 084705.
18. C. E. McNeilly, J. Nucl. Mater., 1964, 11, 53–58.
19. T. M. McCleskey, E. Bauer, Q. Jia, A. K. Burrell, B. L. Scott, S. D. Conradson, A. Mueller, L. Roy, X. Wen, G. E. Scuseria and R. L. Martin, J. Appl. Phys., 2013, 113, 013515.
20. A. J. Garza and G. E. Scuseria, J. Phys. Chem. Lett., 2016, 7, 4165–4170.
21. I. D. Prodan and G. E. Scuseria, J. Chem. Phys., 2005, 123, 014703.
22. Y. Tokunaga, H. Sakai, T. Fujimoto, S. Kambe, R. E. Walstedt, K. Ikushima, H. Yasuoka, D. Aoki, Y. Homma, Y. Haga, T. D. Matsuda, S. Ikeda, E. Yamamoto, A. Nakamura, Y. Shiokawa, K. Nakajima, Y. Arai and Y. Ōnuki, J. Alloys Compd., 2007, 444-445, 241–245.
23. D. Gryaznov, E. Heifets and D. Sedmidubsky, Phys. Chem. Chem. Phys., 2010, 12, 12273–12278.
24. Y. Lu, Y. Yang and P. Zhang, J. Alloys Compd., 2015, 649, 544–552.
25. G. Jomard, B. Amadon, F. Bottin and M. Torrent, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 075125.
26. S. C. Hernandez and E. F. Holby, J. Phys. Chem. C, 2016, 120, 13095–13102.
27. D. L. Clark, S. S. Hecker, G. D. Jarvinen and M. P. Neu, Plutonium, in The Chemistry of the Actinide and Transactinide Elements, ed. L. R. Morss, N. M. Edelstein and J. Fuger, Springer, Dordrecht, 3rd edn, 2008, vol. 2, ch. 7, pp. 813–1264.
28. L.-F. Wang, B. Sun, H.-F. Liu, D.-Y. Lin and H.-F. Song, J. Nucl. Mater., 2019, 526, 151762.
29. M. S. D. Read, S. R. Walker and R. A. Jackson, J. Nucl. Mater., 2014, 448, 20–25.
30. H. Nakamura and M. Machida, Prog. Nucl. Sci. Technol, 2018, 5, 132–135.
31. X. Tian, T. Gao, C. Lu, J. Shang and H. Xiao, Eur. Phys. J. B, 2013, 86, 179.
32. P. Tiwary, A. van de Walle, B. Jeon and N. Grønbech-Jensen, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 094104.
33. M. W. D. Cooper, M. J. D. Rushton and R. W. Grimes, J. Phys.: Condens. Matter, 2014, 26, 105401.
34. M. Stan and P. Cristea, Defects and oxygen diffusion in PuO_2−x, J. Nucl. Mater., 2005, 344, 213–218.
35. M. Kato, H. Nakamura, M. Watanabe, T. Matsumoto and M. Machida, Defect Diffus. Forum, 2017, 375, 57–70.
36. M. W. D. Cooper, S. T. Murthy and D. A. Andersson, J. Nucl. Mater., 2018, 504, 251–260.
37. A. Soulie, F. Bruneval, M. Marinica, S. T. Murphy and J. Crocombette, Phy. Rev. Mater., 2018, 2, 083607.
38. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
39. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251.
40. G. Kresse and J. Furthmüller, Comput. Mater. Sci, 1996, 6, 15.
41. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
42. P. E. Blochl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
43. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
44. H. J. Monkhurst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 49, 16223.
45. S. Steiner, S. Khmelevskyi, M. Marsmann and G. Kresse, Phys. Rev. B, 2016, 93, 224425.
46. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207.
47. J. Heyd and G. E. Scuseria, J. Chem. Phys., 2004, 121, 1187.
48. J. Heyd Scuseria and G. E. Scuseria, J. Chem. Phys., 2006, 124, 219906.
49. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
50. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2009, 100, 136406.
51. G. I. Csonka, J. P. Perdew, A. Ruzsinszky, P. H. T. Philipsen, S. Lebègue, J. Paier, O. A. Vydrov and J. G. Ángyán, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 155107.
52. A. I. Liechtenstein, V. I. Anisimov and J. Zaanen, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 52, R5467(R).
53. A. Shick, J. Kolorenc, L. Havela, T. Gouder and R. Caciuffo, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 041109.
54. F. A. Kröger and H. J. Vink, Solid State Phys., 1956, 3, 307–435.
55. T. Durrant, S. Murphy, M. B. Watkins and A. L. Shluger, J. Chem. Phys., 2018, 149, 024103.
56. M. Leslie and N. J. Gillan, J. Phys. C: Solid State Phys., 1985, 18, 973.
57. S. Baroni and R. Resta, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 7017.
58. M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller and F. Bechstedt, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 045112.
59. L. Manes and A. Barisic, Phys. Status Solidi, 1970, 3, 971.
60. N. Hampton, G. Saunders, J. Harding and A. Stoneham, J. Nucl. Mater., 1987, 149, 18–20.
61. M. S. Talla Noutack, G. Geneste, G. Jomard and M. Freyss, Phys. Rev. Mater., 2019, 3, 035001.
62. G. Makov and M. C. Payne, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 4014.
63. C. Freysoldt, J. Neugebauer and C. G. Van, de Walle, Phys. Rev. Lett., 2009, 102, 016402.
64. Y. Kumagai and F. Oba, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 195205.
65. H.-P. Komsa, T. T. Rantala and A. Pasquarello, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 045112.
66. J. D. Gale, J. Chem. Soc., Faraday Trans., 1997, 93, 629–637.
67. M. W. D. Cooper, S. T. Murphy, M. J. D. Rushton and R. W. Grimes, J. Nucl. Mater., 2015, 461, 206–214.
68. M. Manley, J. Jeffries, A. Said, C. A. Marianetti, H. Cynn, B. M. Leu and M. A. Wall, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 132301.
69. Y. Mishin, M. R. Sørensen and A. F. Voter, Philos. Mag. A, 2001, 81, 2591–2612.
70. D. A. Andersson, P. Garcia, X.-Y. Liu, G. Pastore, M. Tonks, P. Millett, B. Dorado, D. Gaston, D. Andrs, R. L. Williamson, R. C. Martineau, B. P. Uberuaga and C. R. Stanek, J. Nucl. Mater., 2014, 451, 225–242.
71. P. Söderlind, A. Landa and B. Sadigh, Adv. Phys., 2019, 68, 1–47.
72. M. W. Finnis, A. Y. Lozovoi and A. Alavi, Annu. Rev. Mater. Res., 2005, 35, 167–207.
73. G. K. Johnson, E. H. Deventer, O. L. Kruger and W. N. Hubbard, J. Chem. Thermodyn., 1969, 1, 89–98.
74. K. Johnston, M. R. Castell, A. T. Paxton and M. W. Finnis, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 085415.
75. M. Youseff and B. Yildiz, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 144109.
76. S. T. Murphy and N. D. M. Hine, Chem. Mater., 2014, 26, 1629–1638.
77. G. E. Murch, C. Richard and A. Catlow, J. Chem. Soc., Faraday Trans. 2, 1987, 83, 1157–1169.
78. C. Guéneau, C. Chatillon and B. Sundman, J. Nucl. Mater., 2008, 378, 257–272.
79. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
80. S. Ong, W. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. Chevrier, K. Persson and G. Ceder, Comput. Mater. Sci., 2013, 68, 314–319.

---