William D.
Neilson
a,
James T.
Pegg
bc,
Helen
Steele
d and
Samuel T.
Murphy
*a
aEngineering Department, Lancaster University, Lancaster, LA1 4YW, UK. E-mail: samuel.murphy@lancaster.ac.uk
b1QB Information Technologies (1QBit), Vancouver, British Columbia V6E 4B1, Canada
cDepartment of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, UK
dSellafield Ltd, Sellafield, Cumbria, CA20 1PG, UK
First published on 11th February 2021
An increased knowledge of the chemistry of PuO2 is imperative for the design of procedures to store, dispose, or make use of PuO2. In this work, point defect concentrations in PuO2 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 PuO2±x. The results show that the defect chemistry of PuO2 is dominated by oxygen vacancies and interstitials. Hypo-stoichiometric PuO2−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 (O2−i) is found to accommodate hyper-stoichiometry (PuO2+x), but reluctance to form hyper-stoichiometric PuO2+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.
The potential pressurisation of some of the containers used to store PuO2 is of major concern, with multiple pressurisation mechanisms proposed.3 The formation of hyper-stoichiometric PuO2+x may play a role in the chemical reactions producing gas. It is reported that PuO2+x participates in moisture enhanced corrosion of Pu, producing H2,4–6 although this is less of a concern in the UK, where the oxide is stored. Through experimental investigation of PuO2 at conditions relevant to intermediate storage, Haschke et al.4,5 have reported the existence of the hyper-stoichiometric PuO2+x phases, while subsequent work has shown PuO2+x may be best described as a compound containing OH− and H2O species, with excess oxygen accommodated in molecular complexes of Pu(V).7,8 The defect chemistry of PuO2 helps describe the formation and properties of non-stoichiometric PuO2 and the equilibrium of species (e.g. oxygen) between the PuO2 solid and the surrounding gases.
PuO2 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.9 Many of the challenges involved in theoretical calculation arise from trying to correctly describe the 5f electrons. PuO2 is classified as a strongly correlated material, with a ground state of localised 5f electrons.10 In conventional density functional theory (DFT), an over de-localisation of these electrons results in the self-interaction error that causes PuO2 to be described as conducting when in fact PuO2 is classified as a charge-transfer insulator.10
A method used extensively to overcome this shortcoming is to adopt a Hubbard-like correction to give the DFT+U method.10 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 Ueff = U − J in the Dudarev et al. formulism.12 The U and J parameters are obtained by fitting to structural and electronic properties of PuO2.11,13–17 Reproduction of the electronic band gap (Eg) is seen as a strong metric of success. Most previous DFT+U studies of PuO2 fit U and J to reproduce the bandgap of McNeilly et al.18 (Eg = 1.8 eV), however the much more recent experimental study of McCleskey et al.19 reports a higher Eg of 2.8 eV. This large discrepancy in the experimental bandgap for PuO2 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).20 Compared to the DFT+U method, hybrid functionals require significantly more computational effort to utilise, but several studies have applied hybrid functionals to study PuO2 with good experimental replication achieved.14,21
An important aspect in correctly simulating PuO2 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 observations22 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 PuO2.13,17,21,24–26 In their comprehensive study of magnetic order in PuO2, Pegg et al.14 found that 1k AFM states produce an incorrect crystal structure (the correct structure as indicated by experimental results is Fmm crystal symmetry)27 and that the choice of AFM domain (longitudinal or transverse) had a significant impact on the electronic and structural properties of PuO2. This is a key finding, as most studies do not differentiate between collinear magnetic structures. Pegg et al.14 propose a longitudinal 3k AFM ground-state. Of all magnetic configurations tested, this was the lowest in energy, was able to retain Fmm crystal symmetry and replicated well the physical properties of PuO2.
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.14 The experimental measurement of isostructural NpO2 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 PuO2.24,26,28–33 Thermochemical models have been produced that study the defect chemistry of PuO2−x, over a large stoichiometry range,34,35 whilst using DFT+U and a point defect model, Lu et al.24 determined the stability of charged defects within PuO2±x. Lu et al.24 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, Svib. Recent work by Cooper et al.36 and Soulie et al.37 on UO2 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 PuO2 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.
(1) |
(2) |
Following Cooper et al.36 and Soulie et al.,37 two differing atomistic simulation approaches were applied to obtain the parameters of eqn (2): ΔE is found using DFT and ΔSvib 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.
To study the properties of the bulk PuO2, we replicate the study of Pegg et al.14 and use the hybrid Heyd–Scuseria–Ernzerhof (HSE06) functional:46–49
(3) |
For calculation of defect energies in PuO2 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.11 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.52 When using PBEsol+U, the iteration threshold for electronic convergence is 1 × 10−6 eV and for ionic convergence is 1 × 10−2 eV Å−1 and 2 × 10−2 eV Å−1 when simulating PuO2 unit cells and defect containing PuO2 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.20 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 PuO2.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.11 to detrimentally affect the reproduction of the band gap for PuO2. The resulting equilibrium properties found for longitudinal 3k AFM PuO2 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 PuO2 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. |
Defects were inserted into 2 × 2 × 2 expansions of the 12-atom PuO2 unit cell. Only one unique site exists for the vacancies and interstitials of plutonium and oxygen considered, due to the symmetry of the PuO2 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: ×):54
Oxygen interstitials: Oi×, O1−i and O2−i
Oxygen vacancies: , V1+O and V2+O
Plutonium interstitials: , Pu1+i, Pu2+i, Pu3+i and Pu4+i
Plutonium vacancies: , V1−Pu, V2−Pu, V3−Pu and V4−Pu
(4) |
Makov and Payne (MP) extended the PC correction to include a term with L−3 order.62 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 (Vdefect,q & Vbulk respectively).63 The FNV scheme correction energy is summarised following ref. 64 as:
Ecorr = EisoPC − qΔVPC,q/b|far | (5) |
Vq/b = Vdefect,q − Vbulk | (6) |
ΔVPC,q/b = Vq/b − VPC,q | (7) |
Fig. 3 V q/b, VPC,q and ΔVPC,q/b calculated at the atomic positions of a PuO2 supercell with a V4−Pu defect. ΔVPC,q/b|far calculated by averaging ΔVPC,q/b across sampling region. |
(8) |
Fig. 4 The phonon density of states of PuO2 calculated with the CRG empirical potential (blue), compared to experimental data68 (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 PuO2 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 (ΔSvib) and presented in Table 2. The same value of ΔSvib is given to all charge states of a given defect.
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 |
μPu(PO2,T) + μO2(PO2,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(PO2,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.71 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 Å3) matched closely the atomic volume obtained by Söderlind et al.71 when using PBE+U. To determine the chemical potential of oxygen the approach of Finnis et al.72 is adopted. This method uses the known experimental formation energy of the oxide 73 to obtain the chemical potential of oxygen at standard temperature and pressure:
(11) |
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.:74
(12) |
(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 = EVBM + εF, is expressed as the sum of the energy of the valence band maximum (VBM), EVBM, 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:75
(14) |
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 gv(E) and gc(E),the density of electronic states in the valence band and conduction band per formula unit of PuO2, respectively. For calculation of the electron population, ECBM 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 Package76 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 PuO2±x, is calculated from the defect concentrations summed over all charge states:
(15) |
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.77 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.77 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.
Process | Charge states | Reaction energy (eV) | ||||||||
---|---|---|---|---|---|---|---|---|---|---|
DFT+U | Empirical pair potential | Experiment | ||||||||
This work | LSDA+U24 | LDA+U30 | GGA+U31 | GGA+U32 | Ref. 26 | Ref. 27 | Ref. 32 | O diffusion77 | ||
O Frenkel pair | V2+O + O2−i | 4.89 | 3.48 | 4.58 | 9.78 | 3.90 | 5.52 | 2.66 | 4.90 | 2.72–2.92 |
Pu Frenkel pair | V4−Pu + Pu4+i | 13.86 | 15.19 | 10.02 | 19.78 | 11.90 | 13.58 | 10.03 | 24.00 | |
Schottky defect | V4−Pu + 2V2+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−5 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 PuO2.77 All plutonium vacancy and interstitial defects are predicted to have such low concentrations that do not appear in Fig. 6.
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 [Oi] = [VO], meaning that vibrational entropy inclusion further promotes Oi defects. From the point of perfect stoichiometry, increasing or decreasing temperature results in hypo-stoichiometric or hyper-stoichiometric PuO2±x, respectively. The diagrams display the reluctance for PuO2 to be hyper-stoichiometric in oxygen, with x in PuO2+x peaking at approx. 10−7 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 PuO2+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 PuO2.78
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 Oi defects dominate and are charge compensated by holes. The [Oi]:[p+] ratio is approximately 1:2, due to the dominance of the O2−i defect. Increasing the temperature past the point of peak-hyper-stoichiometry, which coincides with peak [Oi], the hole concentration is no longer set such that it provides charge compensation to Oi defects. Instead, the concentration of holes matches the concentration of electrons and the Fermi level remains fixed at the midpoint of the bandgap. [VO] increases with temperature, resulting in increasing hypo-stoichiometry. At high temperatures, [VO] reaches levels close to the concentration of holes and electrons. To accommodate any further increases in [VO], 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 [VO]:[e−] ratio varies in order to provide charge compensation to both V2+O and V1+O defects.
The defect chemistry of PuO2 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 PuO2±x is at its smallest, O2−i is the dominant oxygen interstitial. The dominant oxygen vacancies are found to be V2+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:
(16) |
(17) |
(18) |
Application of the law of mass action to eqn (16), (17) and (18) finds equilibrium constants of k1, k2 and k3, respectfully:
(19) |
(20) |
(21) |
Given that [e−] = [p+], it can be seen that [O2−i] is proportional to , whilst [V2+O] and are proportional to in this region. Increasing the oxygen partial pressure beyond the region of near-stoichiometry sees the concentration holes, O2−i defects and x in PuO2+x transition to an oxygen partial pressure dependence of . With too few electrons available, two holes are created for charge compensation of the O2−i defect. Eqn (22) and (23) show the defect reaction responsible and the corresponding equilibrium constant, respectively. Given that [O2−i] = 2[p+], it can be shown using eqn (23) that [O2−i] is proportional to in the hyper-stoichiometric region.
(22) |
(23) |
Reducing the oxygen partial pressure below the region of near-stoichiometry sees the concentration of V2+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 PuO2−x. The dependence of x in PuO2−x has been shown experimentally to have a dependence on the oxygen partial pressure of ,35 indicative of the V1+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 PuO2 significantly.
The finding that the defect chemistry of PuO2 is dominated by oxygen interstitials and vacancies agrees with the predictions in the PuO2 point defect study of Lu et al.24 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.24 found [O2−i] and [V2+O] defects to be multiple orders of magnitude greater than any non-formally charged defect, here is the dominant oxygen vacancy and concentration of O1−i is predicted to be much closer to that of O2−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−5 atm Lu et al.24 report a value of x in PuO2+x of ≈10−3, several orders of magnitude greater than the value of ≈10−6 found in this study.
Fig. 8 Planar averaged Vq/b, VPC,q and ΔVPC,q/b calculated at the atomic positions of a PuO2 supercell with a VPu4−defect. |
This journal is © the Owner Societies 2021 |