Investigation of oxygen self-diffusion in PuO₂ by combining molecular dynamics with thermodynamic calculations

Abstract

In the present work, the defect properties of oxygen self-diffusion in PuO₂ are investigated over a wide temperature (300–1900 K) and pressure (0–10 GPa) range, by combining molecular dynamics simulations and thermodynamic calculations. Based on the well-established cBΩ thermodynamic model which connects the activation Gibbs free energy of diffusion with the bulk elastic and expansion properties, various point defect parameters such as activation enthalpy, activation entropy, and activation volume were calculated as a function of T and P. Molecular dynamics calculations provided the necessary bulk properties for the proper implementation of the thermodynamic model, in the lack of any relevant experimental data. The estimated compressibility and the thermal expansion coefficient of activation volume are found to be more than one order of magnitude greater than the corresponding values of the bulk plutonia. The diffusion mechanism is discussed in the context of the temperature and pressure dependence of the activation volume.

---

1. Introduction

Nuclear energy is being considered as an efficient way to address future global energy demands, in conjunction with a reduction in greenhouse gas emissions. The primary component of commercial nuclear fuel is UO₂, however this can be mixed with other actinide oxides (i.e. ThO₂ and PuO₂) to form mixed-oxide fuel (MO_x).^1–6 The future of the nuclear industry will necessitate the employment of alternative nuclear fuel cycles and the sustainable use of reprocessed spent fuel in so-called “closed fuel cycles”. Additionally, the use of PuO₂, directly or as a MO_x fuel form, can transform legacy stockpiles of plutonium-rich material to a “strategic fuel reserve” rather than waste. Finally, the technological importance of PuO₂ can be further extended by considering nice applications such as its use as a fuel in thermo-electric converters for space and planetary exploration.⁷

Defect processes, including self-diffusion mechanisms and the tentative formation of more complex defects,^8,9 in actinide oxides nuclear fuels (UO₂, ThO₂ and PuO₂, and their solid solutions) are important in determining the physical properties of the fuel. In particular, oxygen self-diffusion in nuclear fuel is important as it is linked with radiation damage tolerance,¹⁰ clustering of oxygen point defects that can drive the formation of high-burn-up microstructures during operation,¹¹ and finally the migration and solubility of fission products.¹² The use of atomistic simulations has the advantage that it can overcome the limitations and challenges of interacting directly with nuclear materials in reactor conditions. This has been realized by the community for decades and there is now comprehensive data on nuclear fuel materials that can advance our understanding of their properties and/or complement experimental studies.^13–15

The high migration energies of nuclear fuel materials constitute the use of thermodynamic models in conjunction with atomistic simulations necessary to address defect processes at a wide temperature and pressure range (in particular low temperatures).¹⁶ In recent studies, molecular dynamics (MD) calculations have been employed with the so-called cBΩ thermodynamic model. In this model,^17–20 Varotsos and Alexopoulos described Gibbs energy gⁱ as proportional to the isothermal bulk modulus B and the mean volume per atom Ω, with c being the constant of proportionality. The cBΩ model has been applied to describe the defect processes of numerous materials^21–26 including nuclear fuels such as UO₂, PuO₂, ThO₂ and MO_x.^16,27–31 However, the studies of diffusion kinetics in stoichiometric PuO₂ are rather limited as compared to other fluorite-structured nuclear oxides, such as UO₂.

In the present study we investigate the self-diffusion of oxygen in PuO₂ at elevated temperatures (300–1900 K) and pressures (0–10 GPa), by combining the cBΩ thermodynamic model with MD simulations. For the first time, a systematic study of self-diffusion in a nuclear fuel (PuO₂) is carried out in the framework of the cBΩ model and various point defect parameters are estimated as a function of temperature and pressure. The results are discussed in the context of the related diffusion processes and are compared with other similar systems.

2. Methodology

2.1 Molecular dynamics calculations

Classical MD is a powerful tool to simulate the thermo-physical and chemical properties, as well as the kinetics evolution of point defects in various actinide oxide nuclear fuels and MO_x.^15,32–34 In the present study, MD simulations were carried out by means of the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) classical MD code³⁵ using the Cooper–Rushton–Grimes (CRG) interatomic potential set³⁶ which reproduces the thermophysical and thermomechanical properties of numerous actinide oxides, over a wide temperature range (300–3000 K). The CRG potential model is described in detail elsewhere.^36,37 Briefly, this model incorporates many-body interactions using the embedded atom method³⁸ and investigates the variation of bulk modulus as a function of pressure. A PuO₂ supercell of 10 × 10 × 10 fluorite unit cells is used with ions slightly displaced from their lattice sites in random directions. Oxygen point defects were introduced into the simulation supercells, either as vacancies or as interstitials, by removing an oxygen ion from its initial lattice site to the centre of the fluorite cell. The latter, is not actually a stable site, but the time the oxygen spends there, is sufficient to allow another oxygen anion to jump into the space left. In this way, oxygen vacancies are produced by short-lived closely bound Frenkel-pair formation (also known as pseudo-Frenkel pairs due to their spontaneous recombination).^36,37,39 MD calculations were performed in a constant NPT ensemble (pressure, temperature and particle number remain constant) where pressure and temperature were controlled through a Nosé–Hoover thermostat and barostat with relaxation times of 0.1 ps and 0.5 ps, respectively. At a fixed temperature, the pressure is varied from 0 to 10 GPa at intervals of 0.1 GPa and the supercell is equilibrated for 10 ps, with the volume and pressure averaged over the final 5 ps. In this way, a P(V) isothermal curve is obtained which is subsequently fitted with a third-order polynomial. The bulk modulus is then determined through the well known definition, B(T, P) = −V(∂P/∂V), where the first derivative ∂P/∂V is calculated from the analytical polynomial expression of P(V). The calculated data of the pressure and temperature dependence of the atomic volume and the bulk modulus are necessary for the implementation of the cBΩ model, in the lack of any corresponding experimental data.

Our MD calculations are restricted only to these bulk properties since, at the low temperature range of our investigation (300–1900 K) the calculation of the diffusion coefficients is a very time-consuming process with large uncertainties, due to the high migration energy of oxygen defects. Notably, the pressure range (0–10 GPa) used in our MD calculations is significantly higher that the pressure of the fuels under operating conditions inside a nuclear reactor. However, apart from the theoretical interest as concern the pressure dependence of activation volume which is related to the mechanisms of diffusion, high pressures could be developed in crack tips and microstructural defects such as dislocations.⁴⁰

2.2 The cBΩ thermodynamic model

The so-called cBΩ thermodynamic model^{17–20,41,42} relates the activation Gibbs free energy g^act due to the formation and migration of a point defect, with the bulk elastic and expansion properties of the host material, through the relation

g^act = c^actBΩ (1)

where B is the isothermal bulk modulus, Ω is the mean atomic volume and c^act is a dimensionless coefficient which depends on the process and is temperature and pressure independent under certain experimental limitations. Thus, in the case of a single diffusion mechanism, the diffusion process is described by an Arrhenius equation of the diffusion coefficients D, which according to eqn (1) is modified to

D(T, P) = fga₀²νe^{−cactBΩ/k_BT} (2)

In the above equation, f is the diffusion correlation factor which depends on the diffusion mechanism and the structure, g is a geometry factor, a₀ is the lattice parameter, ν is the attempt frequency and k_B is Boltzmann's constant.

According to eqn (1), various point defect parameters such as the activation entropy s^act, activation enthalpy h^act and activation volume υ^act can be derived from the elastic and expansion properties of the bulk material, at any temperature and/or pressure, through the following relations²⁰

where β is the volume thermal expansion coefficient which is also temperature and pressure dependent.

Eqn (3) and (4) provide us with an alternative way of calculating the experimentally determined activation enthalpy h^act_exp and activation entropy s^act_exp in self- or hetero-diffusion processes or predicting them over a broad temperature and pressure range, in the case of limited diffusion data. Activation entropy is indirectly extracted in diffusion experiments from the pre-exponential factor D_o which is derived after the modification of eqn (2). If we replace g^act (=c^actBΩ) in eqn (2) with h^act − Ts^act, then D(T, P) = fga₀²ν exp(s^act/k_B)exp(−h^act/k_BT) and thus, D_o = fga₀²ν exp(s^act/k_B), which is experimentally determined, assuming that the parameters f, g and ν are known. Furthermore, the value and sign of the activation volume (refer to eqn (5)) is important in order to distinguish the kind of the diffusion mechanism, i.e., vacancies or self-interstitials.⁴³

In order to estimate all the previous point defect parameters by means of the cBΩ model (see eqn (3)–(7)), the constant c^act should be determined. In principle, c^act equals to h^act_o/B_oΩ_o where the subscript refers to T = 0 K, since according to eqn (4), g^act = h^act_o at T = 0 K.²⁰ The calculation of c^act can be also achieved in different ways, based on the available experimental diffusion data. The mean value method is reliable when diffusivities are experimentally known over a broad temperature or pressure range.^20,22,44 In the case where only a single diffusion measurement D_i is available at a given temperature T_i (or pressure), the parameter c^act is derived⁴⁵ from eqn (2), under the condition that the pre-exponential factor fga₀ν is roughly known. Subsequently, the diffusion coefficient D can be calculated at any desired temperature T (or pressure), provided that the elastic and expansion data are known at this temperature (or pressure).

3. Results and discussion

The bulk modulus and the mean atomic volume data that were derived from the MD calculations were fitted with similar multi-variable polynomial expressions of temperature and pressure as follows

χ_i(T, P) = χ_o,i + b_iT + c_iT² + d_iP + e_iP² + f_i(TP) + g_i(TP)² (6)

where i = 1, 2 stands for Ω(T, P) and B(T, P), respectively. The parameters χ_o,i denote the mean atomic volume (Ω_o) and the bulk modulus (B_o) at zero pressure and temperature. The parameters b_i, c_i and d_i, e_i express the 1^st and 2^nd order temperature and pressure dependences of χ_i respectively, while, the parameters f_i and g_i account for the interdependency between temperature and pressure. All the above parameters were fitted to the modeling data of bulk modulus and mean atomic volume with very good correlation coefficients in both cases (R = 0.9997) and their values are listed in Table 1. Based on this fitting procedure, the pressure and temperature dependence of Ω(T, P) and B(T, P) in the ranges 300–1900 K and 0–10 GPa are shown as 3D surface plots in Fig. 1(a) and (b), respectively.

Table 1 Fitting parameters of the multi-variable polynomial expressions (see eqn (6)) that describe the atomic mean volume Ω and the bulk modulus B of PuO₂ as a function of pressure and temperature. The parameter d₂ is the pressure derivative of B (dB/dP), at T = 0 and P = 0

	Mean atomic volume, Ω (i = 1)	Bulk modulus, B (i = 2)
χo,i	(12.9738 ± 0.0014) Å^3	(228.36 ± 0.08) GPa
bi	(3.67 ± 0.02) × 10^−4 Å^3 K^−1	(−4.98 ± 0.01) × 10^−2 K^−1 GPa
ci	(4.75 ± 0.07) × 10^−8 Å^3 K^−2	(−7.34 ± 0.50) × 10^−7 K^−2 GPa
di	(−5.84 ± 0.04) × 10^−2 Å^3 GPa^−1	4.69 ± 0.03
ei	(1.04 ± 0.03) × 10^−3 Å^3 GPa^−2	(−0.0750 ± 0.0018) GPa^−1
fi	(−1.53 ± 0.03) × 10^−5 Å^3 GPa^−1 K^−1	(2.05 ± 0.02) × 10^−3 K^−1
gi	(−2.8 ± 1.3) × 10^−11 Å^3 GPa^−2 K^−2	(−2.39 ± 0.08) × 10^−8 K^−2 GPa^−1

---

Fig. 1 3D surface plot of (a) the mean atomic volume Ω and (b) the bulk modulus B of plutonia as a function of pressure (0–10 GPa) and temperature (300–1900 K), as derived from multi-variable fittings of the MD calculated data.

To proceed further to the implementation of the cBΩ model by using the simulated data of Ω(T, P) and B(T, P), the estimation of c^act is necessary. In a recent study of oxygen self-diffusion in PuO₂ in the high temperature range (1800–3000 K), Chroneos et al.²⁷ reported the validity of the cBΩ model and estimated by the mean value method at temperatures below the superionic transition temperature (∼2400 K) the value c^act = 0.3047 ± 0.0075 and the pre-exponential factor, D_o = 1.621 × 10⁻³ m² s⁻¹. These values have been also used in the present case in order to estimate the temperature and pressure dependence of oxygen self-diffusion coefficients and various point defect parameters.

Fig. 2(a) shows the Arrhenius plot of oxygen self-diffusion coefficients D_cBΩ over the temperature range 300–1900 K and at elevated pressures (0–10 GPa), according to the cBΩ model (see eqn (2)). In Fig. 2(b), the diffusion coefficients are plotted as a function of pressure (0–10 GPa), at selected temperatures. The knowledge of the dependence of diffusion coefficients on pressure is associated with the related diffusion mechanisms in the host material, through the estimation of the corresponding activation volumes. By taking the natural logarithm of both sides in eqn (2) and differentiating afterwards with respect to pressure, the following equation for activation volume, υ^act (=(∂g^act/∂P)_T) is derived

Fig. 2 Oxygen self-diffusion coefficients in PuO₂ as a function of (a) inverse temperature (Arrhenius plot) in the temperature range 300–1900 K and (b) pressure (0–10 GPa), according to the cBΩ model. In (b), data are presented at temperatures above 700 K, in order to highlight the non-linear behavior of the diffusion coefficients versus pressure.

By ignoring the second term which is much smaller than the first one, eqn (7) allows the estimation of υ^act from the slope of the plot of ln D_cBΩ versus pressure.^20,45 Thus, the observed non-linear behavior of oxygen self-diffusion coefficients versus pressure in Fig. 2(b) implies that the activation volume υ^act does not remain constant over the entire pressure range at any calculated temperature. As a consequence, the contribution of more than one diffusion mechanism could be possible at elevated pressures. We shall return to this important defect parameter, later in this section.

The temperature and pressure dependence of the activation Gibbs free energy, g^act, (see eqn (1)) are illustrated in Fig. 3(a) and (b), respectively. Activation Gibbs energy decreases with increasing temperature while the effect of pressure is more pronounced as temperature increases. Additionally, g^act increases with increasing pressure while the slope of the curves decreases monotonically (Fig. 3(b)). The latter indicates that the activation volume υ^act also decreases with pressure (υ^act = (∂g^act/∂P)_T), in agreement with the pressure dependence of self-diffusion coefficients in Fig. 2(b).

Fig. 3 Activation Gibbs free energy g^act of oxygen self-diffusion in PuO₂ as a function of (a) temperature (300–1900 K) and (b) pressure (0–10 GPa).

In order to calculate the activation entropy and activation enthalpy from eqn (3) and (4), the volume thermal expansion coefficient β of the bulk material has been estimated from the well known relation β(T, P) = Ω⁻¹(∂Ω/∂T)_P and is depicted in Fig. 4 in a 3D surface plot, as a function of pressure (0–10 GPa) and temperature (300–1900 K). Our estimations of volume thermal expansion coefficient at zero pressure (β = (3.022–3.955) × 10⁻⁵ per °C) are in good agreement with the reported mean value (a = β/3 = 1.116 × 10⁻⁵ per °C) of the linear thermal expansion coefficient of stoichiometric PuO₂, measured from 25 °C to 1420 °C.⁴⁶

Fig. 4 3D surface plot of the volume thermal expansion coefficient β, (β = Ω⁻¹(∂Ω/∂T)_P) of PuO₂, as a function of temperature (300–1900 K) and pressure (0–10 GPa).

The temperature and pressure dependence of activation entropy s^act and activation enthalpy h^act are illustrated in Fig. 5 and 6, respectively. Activation entropy increases with temperature while large values are observed, ranging from 7 k_B to 14.5 k_B, depending on the pressure range. The corresponding high values of the term Ts^act, in conjunction with the activation Gibbs energy (refer to Fig. 3), results in high values of activation enthalpies, h^act (5.6–6.7 eV), over the whole temperature and pressure range (see Fig. 6). These estimated values of h^act are comparable with the activation energies reported recently for oxygen self-diffusion in UO₂ (5.6–6.5 eV) and ThO₂ (5.9–6.9 eV), over the same pressure range (0–10 GPa).^29,31

Fig. 5 Activation entropy s^act (in k_B units) of oxygen self-diffusion in PuO₂, as a function of (a) temperature (300–1900 K) and (b) pressure (0–10 GPa).

Fig. 6 Activation enthalpy h^act of oxygen self-diffusion in PuO₂, as a function of (a) temperature (300–1900 K) and (b) pressure (0–10 GPa).

Finally, the temperature and pressure dependence of activation volume are shown in Fig. 7(a) and (b). The values of activation volume range from less than one mean atomic volume (Ω_o) to about 2.5 Ω_o. A considerable increase in υ^act with increasing temperature appears in isobaric curves which exhibit an upward curvature i.e., at P = 0 and 1 GPa. However, the downward curvature observed at higher pressures restricts the values of υ^act close to Ω_o. In other words, the effect of pressure compensates the role of temperature to the variation of the activation volume.

Fig. 7 Activation volume υ^act of oxygen self-diffusion in PuO₂ as a function of (a) temperature (300–1900 K) and (b) pressure (0–10 GPa). The calculated values are expressed as a function of Ω_o, i.e., the mean atomic volume at zero pressure and temperature.

By taking into account the definition of activation volume υ^act (eqn (5)) in relation to eqn (1), the thermal expansion coefficient β^act and the compressibility κ^act of the activation volume are expressed as follows²⁰

These expressions enable the estimation of these two point defect parameters only from elastic and expansion properties of the bulk material. Their importance has been often underestimated in the literature by assuming that the compressibility of activation volume is equal to that of the bulk (κ^act ≈ κ = 1/B) while the same also holds for the volume thermal expansion coefficient (β^act ≈ β).^30,47,48 These disputed assumptions are of considerable importance in the field of geophysics where at the lower mantle of the earth, the pressures are of the order of the bulk moduli of the constituent minerals⁴⁹ and the dominant mechanism of deformation is diffusion-controlled creep.⁵⁰

In the present case of oxygen self-diffusion in plutonia, the quantities β^act and κ^act have been plotted as a function of temperature and pressure in Fig. 8 and 9, respectively. Two key points need to be emphasized. Firstly, the values of β^act and κ^act are more than one order of magnitude higher than the corresponding values of bulk plutonia. Secondly, a peculiar behavior of the above quantities is observed around 700 K. Specifically, β^act varies from 31.0 to 50.6 × 10⁻⁵ per °C at zero pressure and from 8.3 to 65.8 × 10⁻⁵ per °C at the maximum pressure, P = 10 GPa. At the same time, κ^act lies between (4.05–5.97) × 10⁻¹¹ Pa⁻¹ at room temperature and (5.02–7.90) × 10⁻¹¹ Pa⁻¹ at 1900 K, depending on the applied pressure. The ranges of the corresponding compressibility values, κ, of the bulk plutonia are (0.47–0.39) × 10⁻¹¹ Pa⁻¹ and (0.76–0.50) × 10⁻¹¹ Pa⁻¹ (refer to Fig. 1(b)). As concern the thermal expansion coefficient of the activation volume, at temperatures below 700 K, the effect of pressure causes the increase of β^act while at temperatures higher than 700 K, β^act decreases with pressure (refer to Fig. 8). The isotherm of 700 K corresponds to the minimum values of κ^act with respect to the other isotherms, over the entire pressure range. This “critical” temperature is low enough to be associated with any possible phase transition of the bulk but coincides roughly to the turning point in the plot of the specific heat of plutonia versus temperature, where the slope changes from increasing to decreasing.⁵¹ Since, specific heat (at constant pressure) is related directly with the thermal expansion coefficient, we expect these quantities to exhibit similar characteristics.

Fig. 8 3D surface plot of the thermal expansion coefficient of the activation volume β^act, as a function of temperature (300–1900 K) and pressure (0–10 GPa).

Fig. 9 3D surface plot of the compressibility κ^act of the activation volume as a function of temperature (300–1900 K) and pressure (0–10 GPa).

Large differences between β^act, κ^act and β, κ have also been reported by Sarlis and Skordas, in the case of oxygen self-diffusion in UO₂.³⁰ However, there is no indication of any reversible effect of temperature to β^act and κ^act, in their calculations. The compressibility of the activation volume for vacancy motion in the alkaline earth fluoride, CaF₂, which is isomorphous with actinide oxides, has been found to be forty times the compressibility of the host material. Additionally, a large (but negative) value of thermal expansion coefficient of activation volume has been reported at the same time.⁵² These reported values are consistent with our findings, suggesting that the fluorite structure with the low ionic packing factor should be responsible for the high compressibility of the activation volumes in each case. In many cases such as in metals or alkali halides, the compressibility of activation volume has been reported to be even 2–9 times greater than the corresponding value of the bulk, at zero pressure.^53–55 The fact that κ^act > κ arises from eqn (9), where (∂B/∂P)_T is always much greater than unity, while (∂²B/∂P²)_T takes negative values.^20,55

Recalling that the activation volume corresponds to the isothermal volume change of the crystal due to the formation and migration of a defect, we expect that, in the case of an oxygen vacancy formation that is accompanied by an outward relaxation, due to the electrostatic repulsive forces between the adjacent Pu cations, the increase of temperature (at ambient pressure) should result in an increment of the volume available to the vacancy, due to the lattice thermal expansion.⁴³ Taking into account that the volume change between the equilibrium state and the excited state at the saddle-point (migration volume) contributes further to the total change of activation volume, the thermal expansion coefficient of υ^act is expected to increase more than that of the crystal. Consequently, the compressibility of the activation volume which is related to the vacancies should also increase.

The validity of the cBΩ model has been recently verified in the cases of oxygen self-diffusion in actinide oxides (ThO₂, UO₂ and PuO₂) over a wide pressure and temperature range.^27–29,31 Nevertheless, in all cases, the estimations of the point defect parameters were restricted only to activation energies and activation volumes, without any extensive calculations over their reported pressure and temperature range. Our estimated values of the point defect parameters according to the cBΩ model, in conjunction with the previous reported values and some representative experimental and theoretical estimations of diffusivities in actinide oxides,^56–59 are summarized in Table 2. Our estimations of activation enthalpies are comparable with the theoretical predictions but deviate considerably from the experimental values that are scattered at a lower energy range. This discrepancy could be attributed to some extent, to the fact that the experimental measurements are referred to nominally stoichiometric oxides since an exact stoichiometry is impossible to be achieved.^56,59

Table 2 Calculated values of the parameter c^act, pre-exponential factor D_o, activation enthalpy h^act, activation entropy s^act, activation Gibbs energy and activation volume υ^act in actinide oxides (plutonia, thoria and urania) from the present and previous studies, according to the cBΩ model, in comparison with reported experimental and theoretical results

Oxide	Temperature/pressure	c^act	Do (m^2 s^−1)	h^act (eV)	s^act (kB)	g^act (eV)	υ^act (×10^−29 m^3)
a This study.b The mean atomic volume of PuO2 at T = 0 K, P = 0 GPa, is Ωo = 1.29738 × 10^−29 m^3.
PuO2	300–1900 K, 0–10 GPa^a	0.3047 ± 0.0075	1.621 × 10^−3	(5.64–6.65) ± 0.61	(7.27–14.48) ± 0.72	(3.45–6.18) ± 0.27	(1.06–3.20) ± 0.11^b
	1800–2400 K,^27 2400–3000 K (ref. 27)	0.3047 ± 0.0075, 0.0604	1.621 × 10^−3, 2.234 × 10^−8	—	—	—	—
	993–1293 K (ref. 56)	—	1.30 × 10^−6	2.09	—	—	—
	1000–1290 K (ref. 57)	—	1.19 × 10^−7	0.44	—	—	—
	773–1373 K (ref. 58)	—	1.87^+1−1.6 × 10^−6	1.94 ± 0.29	—	—	—
ThO2	300–1900 K, −10 to 10 GPa (ref. 31)	0.3388	3.442 × 10^−4	4.75–6.90	—	—	1.93–2.12
	1100–1473 K,^56 1473–1900 K (ref. 56)	—	1.00 × 10^−10, 5.73 × 10^−6	0.77, 2.16	—	—	—
UO2	300–1900 K, 0–10 GPa (ref. 29)	0.3052	1.277 × 10^−4	5.65–6.52	—	—	1.45–1.77
	∼2000–2400 K (ref. 60)	—	—	5.7	—	—	—
	1548–1923 K (ref. 61)	—	2.6 × 10^−5	2.6	—	—	—

---

To the best of our knowledge, the only reported values of activation volumes were also derived in the basis of the cBΩ model for UO₂ and ThO₂, over the temperature range, 700–1500 K (see Table 2), but without any estimation of their pressure dependence.^29,31 Our estimated values (2.07–2.81 × 10⁻²⁹ m³) over the same temperature range are comparable with these reported values. In a first thought, activation volumes close to Ω_o (∼1.3 × 10⁻²⁹ m³) could be consistent with a vacancy-mediated diffusion, since the formation of an oxygen vacancy corresponds to a formation volume close or greater than Ω_o, considering the contribution of the positive relaxation volume around this defect. At this point, we should mention that the observed curvature of the diffusion coefficients with pressure (see Fig. 2(b)) could be solely attributed to the pressure dependence of the activation volume only if a linear relation between ln D and υ^act is valid.^20,62 However, this is excluded in our case, since υ^act varies in a non-linear way with pressure (see Fig. 7(b)) and thus the plot of ln D vs. υ^act deviates from linearity, over the entire pressure range, at any given temperature. Thus, the diffusion process should be more complicated than a single vacancy mechanism and an additional process should be taken into consideration. This is compatible with the vacancy-mediated diffusion through the formation of short-lived Frenkel pair formation (pseudo-Frenkel pair) that has been considered in our MD calculations for actinide oxides.^36,37 Thus, in the case of PuO₂, the migration pathway includes the energy to create the pseudo-Frenkel pair (formation energy) and the energy to move an oxygen vacancy (migration energy). For example, in the case of UO₂ and ThO₂, the oxygen pseudo-Frenkel energies were estimated to be 4.42 eV and 4.70 eV, respectively.³⁷ As concern the corresponding activation volumes in PuO₂, at ambient pressure, the creation of the pseudo-Frenkel pairs corresponds to a formation volume close to zero,⁶³ while the outward relaxation of the lattice should be positive but small. So, our high estimated values of υ^act are attributed to the high values of its thermal expansion coefficient, rather than the diffusion mechanism itself. The thermal expansion of the lattice at elevated temperatures is accompanied by a drastic change of the thermal expansion of activation volume and thus, υ^act increases considerably. However, at higher pressures, the temperature-induced volume change is smaller and we expect the activation volume to increase to a lesser extent than at zero pressure.

4. Conclusions

In this work, a detailed investigation of oxygen self-diffusion in PuO₂ has been carried out in the frame of the cBΩ thermodynamic model which allows the estimation of important point defect parameters, such as activation enthalpy, activation volume, activation compressibility and thermal expansion coefficient of activation volume, as a function of pressure and temperature. All the necessary “input” data for the implementation of the model i.e., the bulk modulus and the mean atomic volume, were derived from classical MD simulations over the ranges 300–1900 K and 0–10 GPa, due to the lack of relevant experimental data. Our estimations of the activation enthalpies and the corresponding activation volumes are compatible with the dominant self-diffusion mechanism of oxygen vacancies migration through the formation of pseudo-Frenkel-pairs. It is noteworthy that, the compressibility and the thermal expansion coefficient of the activation volume are about one order of magnitude greater than the corresponding values of the bulk material, a finding that had been underestimated in many cases since now. Thus, the knowledge of υ^act alone is not always sufficient to describe correctly the related diffusion mechanism. The combination of MD calculations with the cBΩ thermodynamic model is proved to be a powerful tool for the investigation of diffusion kinetics in plutonia. A detailed analysis in the frame of the cBΩ model should be also carried out for self- or hetero-diffusion in other fluorite-structured or minor actinide oxides, in order to ascertain any possible similarities or differences among them.

Acknowledgements

A. C. and M. E. F. are grateful for funding from the Lloyd's Register Foundation, a charitable foundation helping to protect life and property by supporting engineering-related education, public engagement and the application of research. Work performed by M. W. D. C. was funded by the U.S. Department of Energy, Office of Nuclear Energy, Nuclear Energy Advanced Modeling Simulation (NEAMS) program. Los Alamos National Laboratory, an affirmative action/equal opportunity employer, is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. DE-AC52-06NA25396.

References

1. M. Lung and O. Gremm, Nucl. Eng. Des., 1998, 180, 133.
2. M. Kazimi, Am. Sci., 2003, 91, 408.
3. M. Stan, Mater. Today, 2009, 12, 20.
4. W. E. Lee, M. Gilbert, S. T. Murphy and R. W. Grimes, J. Am. Ceram. Soc., 2013, 96, 2005.
5. S. C. Middleburgh, G. R. Lumpkin and R. W. Grimes, Solid State Ionics, 2013, 253, 119.
6. O. S. Valu, O. Benes, R. J. M. Konings, J. C. Griveau and E. Colineau, J. Phys. Chem. Solids, 2015, 86, 194.
7. A. Witze, Nature, 2014, 515, 484.
8. C. A. Londos, N. Sarlis, L. G. Fytros and K. Papastergiou, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 6900.
9. N. Sarlis, C. A. Londos and L. Fytros, J. Appl. Phys., 1997, 81, 1645.
10. K. E. Sickafus, R. W. Grimes, J. A. Valdez, A. Cleave, M. Tang, M. Ishimaru, S. M. Corish, C. R. Stanek and B. P. Uberuaga, Nat. Mater., 2007, 6, 217.
11. T. Ichinomiya, B. P. Uberuaga, K. E. Sickafus, Y. Nishiura, M. Itakura, Y. Chen, Y. Kaneta and M. Kinoshita, J. Nucl. Mater., 2009, 384, 315.
12. M. Bertolus, M. Freyss, B. Dorado, G. Martin, K. Hoang and S. Maillard, et al., J. Nucl. Mater., 2015, 462, 475.
13. R. W. Grimes and C. R. A. Catlow, Proc. R. Soc. London, Ser. A, 1991, 335, 609.
14. K. Govers, S. Lemehov, M. Hou and M. Verwerft, J. Nucl. Mater., 2007, 366, 161.
15. A. Chroneos, M. J. D. Rushton, C. Jiang and L. H. Tsoukalas, J. Nucl. Mater., 2013, 441, 29.
16. D. C. Parfitt, M. W. D. Cooper, M. J. D. Rushton, S.-R. G. Christopoulos, M. E. Fitzpatrick and A. Chroneos, RSC Adv., 2016, 6, 74018.
17. P. Varotsos and K. Alexopoulos, Phys. Rev. B: Condens. Matter Mater. Phys., 1977, 15, 4111.
18. P. Varotsos and K. Alexopoulos, Phys. Rev. B: Condens. Matter Mater. Phys., 1977, 15, 2348.
19. P. Varotsos and K. Alexopoulos, Phys. Rev. B: Condens. Matter Mater. Phys., 1980, 22, 3130.
20. P. Varotsos and K. Alexopoulos, Thermodynamics of Point Defects and their Relation with the Bulk Properties, North-Holland, Amsterdam, 1986.
21. P. Varotsos, Solid State Ionics, 2008, 179, 438.
22. B. H. Zhang and X. P. Wu, Appl. Phys. Lett., 2012, 100, 051901.
23. V. Saltas and F. Vallianatos, Mater. Chem. Phys., 2016, 181, 204.
24. V. Saltas, A. Chroneos and F. Vallianatos, RSC Adv., 2016, 6, 53324.
25. E. S. Skordas, Solid State Ionics, 2014, 261, 26.
26. M. Lazaridou, C. Varotsos, K. Alexopoulos and P. Varotsos, J. Phys. C: Solid State Phys., 1985, 18, 3891.
27. A. Chroneos, M. E. Fitzpatrick and L. H. Tsoukalas, J. Mater. Sci.: Mater. Electron., 2015, 26, 3287.
28. A. Chroneos and R. V. Vovk, Solid State Ionics, 2015, 274, 1.
29. M. W. D. Cooper, R. W. Grimes, M. E. Fitzpatrick and A. Chroneos, Solid State Ionics, 2015, 282, 26.
30. N. V. Sarlis and E. S. Skordas, Solid State Ionics, 2016, 290, 121.
31. M. W. D. Cooper, M. E. Fitzpatrick, L. H. Tsoukalas and A. Chroneos, Mater. Res. Express, 2016, 3, 065501.
32. D. S. Aidhy, P. C. Millett, T. Desai, D. Wold and S. R. Phillpot, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 104107.
33. R. Devanathan, L. Van Brutzel, A. Chartier, C. Guéneau, A. E. Mattsson, V. Tikare, T. Bartel, T. Besmann, M. Stane and P. Van Uffelen, Energy Environ. Sci., 2010, 3, 1406.
34. D. Rupasov, A. Chroneos, D. Parfitt, J. A. Kilner, R. W. Grimes, S. Y. Istomin and E. V. Antipov, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 172102.
35. S. Plimpton, J. Comput. Phys., 1995, 117, 1.
36. M. W. D. Cooper, M. J. D. Rushton and R. W. Grimes, J. Phys.: Condens. Matter, 2014, 26, 105401.
37. M. W. D. Cooper, S. T. Murphy, P. C. M. Fossati, M. J. D. Rushton and R. W. Grimes, Proc. R. Soc. A, 2014, 470, 20140427.
38. M. S. Daw and M. I. Baskes, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 29, 6443.
39. F. Devynck, M. Iannuzzi and M. Krack, Phys. Rev. B: Condens. Matter Mater., 2012, 85, 184103.
40. D. C. Parfitt, C. L. Bishop, M. R. Wenman and R. W. Grimes, Phys. Rev. B: Condens. Matter Mater., 2010, 22, 175004.
41. P. Varotsos, W. Ludwig and K. Alexopoulos, Phys. Rev. B: Condens. Matter Mater. Phys., 1978, 18, 2683.
42. K. Alexopoulos and P. Varotsos, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 24, 3606.
43. M. J. Aziz, Defect Diffus. Forum, 1998, 153–155, 1.
44. V. Saltas and F. Vallianatos, Mater. Chem. Phys., 2015, 163, 507.
45. F. Vallianatos and V. Saltas, Phys. Chem. Miner., 2014, 41, 181.
46. M. Tokar, A. W. Nutt and T. K. Keenan, Nucl. Technol., 1973, 17, 147.
47. S. Karato, Phys. Earth Planet. Inter., 1981, 24, 1.
48. D. R. Mills, S. C. Parker and A. Wall, Philos. Mag. A, 1991, 64, 1133.
49. N. V. Sarlis and S. Skordas, J. Phys. Chem. A, 2016, 120, 1601.
50. J. P. Poirier and R. C. Liebermann, Phys. Earth Planet. Inter., 1984, 35, 283.
51. F. L. Oetting, J. Nucl. Mater., 1982, 105, 257.
52. J. J. Fontanella, M. C. Wintersgill and C. Andeen, Phys. Status Solidi B, 1980, 97, 303.
53. P. Varotsos, W. Ludwig and C. Falter, J. Phys. C: Solid State Phys., 1978, 11, L311.
54. P. A. Varotsos, J. Phys., Lett., 1977, 38, L455.
55. P. Varotsos and W. Ludwig, J. Phys. C: Solid State Phys., 1978, 11, L305.
56. K. Ando and Y. Oishi, J. Nucl. Sci. Technol., 1983, 20, 973.
57. R. L. Deaton and C. J. Wiedenheft, J. Inorg. Nucl. Chem., 1973, 35, 649.
58. A. S. Bayoglu, A. Giordano and R. Lorenzelli, J. Nucl. Mater., 1983, 113, 71.
59. G. E. Murch, C. Richard and A. Catlow, J. Chem. Soc., Faraday Trans. 2, 1987, 83, 1157.
60. T. Arima, K. Yoshida, K. Idemitsu, Y. Inagaki and I. Sato, IOP Conf. Ser.: Mater. Sci. Eng., 2010, 9, 012003.
61. J. F. Marin and P. Contamin, J. Nucl. Mater., 1969, 30, 16.
62. M. Kumari and N. Dass, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 844.
63. J. Oberschmidt and D. Lazarus, Phys. Rev. B: Condens. Matter Mater. Phys., 1980, 21, 5813.

---