Modeling Superionic Behavior of Plutonium Dioxide

Introduction

Plutonium dioxide PuO2, generally used in mixed compound (MOX), is an important material as a fuel and a stable storage instrument in nuclear reactors. Our knowledge of this material is limited due to toxicity and radioactivity, which make it difficult to study. The melting temperature of PuO₂is an important parameter in studying phase diagram of its mixed oxide, generally with UO₂, at high temperature. There have been a few experiments that aim to determine the melting temperature of PuO₂ [1–3], which was apprised to be around 2,700 K. However, during the traditional heating techniques, the measurements were affected either by the reaction of oxygen sample with the metallic containment or by the loss of oxygen in the atmosphere [4–6].

Recently, the melting behavior of stoichiometric PuO₂ has been studied for the first time by the fast laser heating and multi-wavelength spectro-pyrometry to reduce the side effect mentioned earlier [6]. Based on this measurement the melting point of PuO₂ has been determined as 3,017±28K. Moreover, the temperature evaluation of physicochemical properties of PuO₂ are limited up to 2,000 K which is well below the melting temperature [7,8]. Even though a diffuse transition called Bredig transition is a highly expected feature of all the actinide oxide yet it has not been proved experimentally for PuO₂. In order to observe whether PuO₂ display the Bredig transition to the superionic phase or not, one has to produce some physicochemical properties up to the melting point. However, during the enthalpy measurements of PuO₂ by Ogard, partial melting of sample above 2,370 K was observed that has been thought, from the private communication [7], as a result of the reaction between the container and PuO₂. Therefore, Fink did not include the phase transition in the non-linear least-square fit to the experimental data due to the large scatter data points between 2,160 and 2,370 K and lack of data above 2,370 K. Later on, Ralph [8] has applied the method of quasi-local linear regression (QLLR) to the enthalpy data by Ogard and no peak has been produced in both the real Cp and the artificial Cp. It is interpreted that the absence of the Cp peaks may be the result of the systematic error in the regression together with the random error in the data.

At this point, molecular dynamics simulation turns out to be a useful tool for studying the properties that may be essential for nuclear facilities to run safely. Molecular dynamics simulation of such systems requires developing a reliable potential for the interactions. Within the nuclear fuel materials, uranium dioxide UO₂ is one of the most studied systems. There have been a number of pair potentials developed to understand thermophysical and transport properties of UO₂ at solid and liquid phases. Govers et al. [9,10] have made a broad comparison to assess the range of applicability of these potentials. For PuO₂, previous studies have come up with some pair potentials, generally in the form of Buckingham, Morse, embedded atom model (EAM) and shell model to investigate PuO₂ [11–15]. Recently, Cooper et al. [15] have reported a many-body potential model to describe the thermomechanical properties of actinide oxides between 300 and 3,000 K. However, the Bredig transition has not been considered at any stage of this work.

PuO₂, like other actinide dioxides, is believed to be a type II superionic conductor, which displays a rapid but continuous increase in the ionic conductivity by heating at about 0.8Tm. The extent of anion disorder increases within the superionic phase. Ionic conductivity occurs almost entirely with anion diffusion. Although the phase transition is well above the operating temperature of nuclear facilities, anion Frenkel pair concentration begins to increase at this level. Oxygen diffusion rapidly increases around 2,160–2,370 K, while Pu ions show almost no diffusive behavior up to 3,000 K. This behavior is a well-known second-order phase transition that reflects as a λ peak in specific heat at constant pressure Cp at the transition temperature. After 3,000 K, Pu ions also show an increase in ionic conductivity comparable with oxygen ions that signal the system to undergo a first-order phase transition, namely melting transition. So the detailed knowledge of the Bredig transition and the melting of PuO₂ at atomic level is important in order to enhance the nuclear fuel performance and the accident conditions.

The latest value of the melting temperature of PuO₂ determined by laser-heating technique gave us the opportunity to investigate the thermally induced Bredig transition and the melting. For this purpose, with the present article, a new semi-empirical rigid ion potential is introduced, and thermophysical properties calculated via classical molecular dynamics simulation in constant pressure–temperature ensemble are presented for the wide range of temperature (300—3,600 K).

Potential model

There exist a few semi-empirical potential model in literature developed for PuO₂ [11–15]. Yamada, Kurosaki and Arima [11–13] have used Morse- and Buckingham-type pair interaction potential function. Buckingham type with shell model pair interaction potential has been introduced by Chu [14]. Cooper added the term many-body EAM to Buckingham- and Morse-type pair interaction potentials.

In general, the parameters for these potential functions O–O interaction have been generally taken from the potential developed earlier forUO₂. On the other side, Pu–Pu and Pu–O parameters have been determined by trial-and-error method, aiming to reproduce the experimental lattice constant at low temperatures. In obtaining the parameters of the pair potential in the present study, we also take the O–O parameters from the potential developed for UO₂ in our previous work [16]. The other parameters were adjusted in order to reproduce low-temperature experimental values of lattice parameter, bulk modulus, elastic constants and cohesive energy. Moreover, we have also aimed to model the Bredig transition to the superionic phase, which is highly expected but has not been yet observed experimentally. The potential function has been originally proposed by Vashishta-Rahman [17] and is given by

The first term stands for the Coulomb interactions, second term contains core repulsions, third term is the effective monopol-induced dipole interaction and the last term is Van der Waals interaction. The potential parameters are given in Table 1.

Molecular dynamics simulation

Molecular dynamics technique was performed for 5 × 5 × 5 cell constructed with 500 cations (Pu⁴⁺) and 1,000 anions (O^2–) arranged as CaF₂-type structure. The MD code called as MOLDY [18] was adopted to carry out the calculations. The Ewald’s sum technique is used to account for the long-range Coulomb interactions. The positions and velocities of the ions are calculated by integrating the Newton’s equation of motion using Beeman’s algorithm, which is predictor–corrector type, with the time step Δt=1.0×10−15s. Nose-Hoover thermostat and Parrinello–Rahman constant pressure methods are applied to control the temperature and pressure in the constant pressure–temperature (NPT) ensemble. The simulations were run for 50 ps: the first 10 ps was used for the equilibration process and the properties of the system were calculated by averaging over the following 40 ps. The temperature is varied from 300 to 3,600 K at 100 K intervals. Additional calculations have been carried out with different number of atoms each for 50 K temperature intervals: 864 (Pu⁴⁺)–1,728 (O^2–), 1,372 (Pu⁴⁺)–2,744 (O^2–) and 2,048 (Pu⁴⁺)–4,096 (O^2–); 6×6×6, 7×7×7 and 8×8×8, respectively, to make sure that number of atoms are enough to calculate thermodynamic properties and also to be able to determine the type of the transitions. In order to calculate properties, these simulations have been performed for 100 ps, where the equilibration run is 30 ps and data is accumulated for 70 ps.

Results and discussion

The values of lattice parameter, bulk modulus, elastic constants and cohesive energy are given in Table 2 and compared with the available experimental data and ab initio results. The bulk modulus and elastic constants were calculated from the utilization of known Birch–Murnaghan equation of state [19, 20], which is in good agreement with experiment and ab initio results [21–25]. Contrary to the previous MD results [11–15], in the present paper the lattice parameter is underestimated about 1.3 % and 3 % compared to the ab initio calculations and experimental data, respectively [22, 23]. However, bulk modulus obtained with the previous potentials has been overestimated and varied between 200 and 239 GPa [11–15] when compared with the experimental result 178 GPa [21]. In the present study, the amount of deviation of bulk modulus from the experimental study is about 5 %, whereas other studies are between 12 % and 34 %. The main reason for these results is the weight of observables while finding the potential parameters. Groups place great emphasis on the lattice constant while developing potential parameters. In this study, we distribute the weight constant unevenly on lattice constant, bulk modulus, elastic constants and energy as we are performing fitting procedure. Temperature dependence of energy is also taken into account in order to mimic the fast-ion phase transition. Discrepancy between the calculated value and experimental data of C₁₂ are presented in Table 2. This surely can be modified by changing the potential parameters but in that the indication of the phase transition can be lost.

L/L0 is calculated, where L0 is the lattice constant at 300 K, as a function of temperature together with the other simulation results [11–15] and the experimental data [22], which are available up to 2,000 K are presented in Figure 1. Figure 2 shows the linear thermal expansion ΔL/L0. All potentials produce a similar thermal expansion of the lattice between 300 and 2,000 K; discrepancy between the simulation and the experiment increases as the temperature increased. At about 2,100 K, there is a sign of a change in trend of the evolution of the data with temperature. As recommended by Fink [7], transition to the superionic phase may be expected at temperature between 2,160 and 2,370 K.

Figure 1:

Lattice parameter evolution with temperature. Experimental data are from TPRC data series [22].

Figure 2:

Evolution of relative linear thermal expansion of present results, MD data and experimental data with temperature.

The temperature dependence of the pair correlation functions gijr given in Figure 3 was also considered in our calculations aiming to have quantitative measure of the spatial local order–disorder of the atomic structure of PuO₂.

Figure 3:

Pair correlation functions of PuO2 at three different temperatures.

The radial distribution function is defined as follows:

Δn(r) is the average number of atoms at a distance between r and r+Δr. ρ is the atomic density and Δr is the width of the shell. The peaks of gijr of all three pairs of ions become lower and broader due to the larger thermal vibrations of ions in their lattice sites as the temperature increased up to 2,000 K. Right after 2,000 K, goor has very small oscillations and overlap of the principal peaks after the first peak, indicating a considerable degree of disorder of the oxygen sublattice while the plutonium ions are still located in the lattice sites with larger vibrations. This liquid-like behavior of the oxygen sublattice at solid phase at about 2,100 K is interpreted as the onset of the transition to the superionic phase of PuO₂.

For a species of N particle, the mean square displacement (MSD) is calculated as

where rn(t) is the position of particle n at time t. The MSD of oxygen ions at different temperatures is shown in Figure 4. The linear change of MSD with time indicates that the oxygen ions become more diffusive at about 2,100 K, while the MSD of Pu ions remains constant. Figure 5 illustrates the diffusion coefficient of both Pu and O ions in Arrhenius diagram (log D versus 1/T) that are calculated from the slope of the MSD data using Einstein relation

Figure 4:

Mean square displacement versus time at different temperatures.

Figure 5:

Evolution of diffusion coefficient with temperature for both plutonium and oxygen ions.

It is clearly evident from the anomalous increase of diffusion coefficient of oxygen ions at ~2,100 K that there is an onset of Bredig transition to the superionic phase. While the oxygen diffusion continuously increases beyond 2,100 K, the plutonium ions show almost no diffusive behavior. As experimental data for oxygen diffusion coefficient for PuO₂ is only available up to 1,400 K [26], we are not able to compare self-diffusion coefficient. As the temperature increases, Pu ions also show diffusive behavior after 3,000 K and D_Pu is less than D_O about an order of 10.

The enthalpy change ΔH=HT−H298 is presented in Figure 6 and compared with the experimental data taken from Fink [7] and MD simulation with shell model potential by Chu. Evidently, there is a discontinuity in calculated ΔH at about 2,100 K and very good agreement with experimental data before and after this temperature. This justifies the potential parameters and validates the constructed model of PuO₂.

Figure 6:

Enthalpy change of PuO2 with temperature.

The existence of the thermally activated transition into superionic phase can be confirmed by a λ-peak in the heat capacity at constant pressure CpT. Many fluorite-type ionic crystals exhibit such a transition at about 0.8Tm [27]. The heat capacity was evaluated from the variation of the internal energy ET with temperature at constant pressure and presented in Figure 7:

Figure 7:

Temperature dependence of heat capacity for PuO2 for the present study at 100 K intervals and comparison with experimental data. The inset shows the results for the different box sizes at 50 K intervals.

A very weak increase in C_p up to 1,750 K is interpreted as increase in the anharmonicity of the lattice vibrations. Further increase in temperature results with the additional increase in calculated C_p and a λ shape peak is clearly produced at the critical temperature T_c = 2,055 K, which is close to the expected value of the phase transition temperature in fluorite type of ionic crystals [7]. In order to make sure that it is a λ-peak transition in specific heat, we have increased the number of data by varying the temperature at 50 K intervals in NPT ensemble for MD simulation calculations. These calculations have been carried out for three different supercell boxes: 6 × 6 × 6 (864 Pu-1728 O), 7 × 7 × 7 (1372 Pu-2744 O) and 8 × 8 × 8 (2048 Pu-4096 O) to show the consistency of the calculations in terms of number of atoms. Results for heat capacity at constant pressure around the transition region are shown in Figure 7 as inset. It is clearly evident that the system undergoes a thermally activated λ-type Bredig transition to the superionic phase. Moreover, the linear increase of MSD of oxygen ions with time at about 2,100 K while that of plutonium remains constant presented in Figure 4 and anomalous increase of diffusion coefficient of oxygen ions at around the same temperature presented in Figure 5, that are due to the premelting of oxygen sublattice, also support the indication of lambda transition in heat capacity.

Conclusions

A model of PuO₂ as a superionic conductor has been constructed via parameterization of a rigid ion potential with partially ionized atoms. The thermophysical properties of the model have been investigated up to 3,600 K within the classical molecular dynamics in the NPT ensemble. A sign of an anomalous behavior in the lattice parameter and oxygen–oxygen pair correlation function was observed at about 2,100 K. This was reflected as a sharp increase in the self-oxygen diffusion coefficient to the value of a typical liquid system, a discontinuity in enthalpy change and as a λ-type peak in the constant pressure heat capacity at about the same temperature. These features indicate that the system undergoes a thermally activated Bredig transition to the superionic phase and validates the constructed superionic model of PuO₂. As it is also pointed out by Fink [7], more accurate high-temperature enthalpy measurements are needed to put this issue beyond the doubts.