Nuclear fuels are almost exclusively actinide-based, with thorium, uranium and plutonium being the most common elements used to provide the fissile material. Of the three, the most abundant substance is uranium, a relatively common mineral in the Earth’s crust, which occurs in low concentrations in most rocks, soils and sea water, while it is primarily found in uraninite or pitchblende deposits. Uranium exists in U⁴⁺, U⁵⁺, and U⁶⁺ oxidation states when in combination with oxygen, thus UO₂ forms the lowest available stoichiometry and UO₃ the highest. However, as shown in Figure 1, the uranium-oxygen system is complex, with a large number of intermediate oxides between these two compositions and a number of different polymorphs at each stoichiometry. These include U₄O₉ and U₃O₇, which are closely related to the cubic fluorite structure of UO₂, as well as U₃O₈ and UO₃ which adopt distinct layered-type structures. Bridging the gap between these is U₂O₅ reported to have a mixture of fluorite-type and layered-type polymorphs. Over the years a considerable body of research has amassed that deals with fluorite structured UO₂, due to its role as the primary fuel material in nuclear reactors. However, the literature is lacking on fundamental research into the higher uranium oxides which also play crucial roles in the fuel cycle and, in the case of both U₃O₈ and UO₃, could soon be used as fuels themselves (Hopper et al., 2002). Although more work is starting to emerge on these two materials, there is still a significant paucity of data on the less stable U₄O₉, U₃O₇ and U₂O₅ oxides, which could provide valuable insight into the process of oxidation in UO₂ and fluorite-based structures. Thus, there is broad scope for more detailed investigation of the different polymorphs at each stoichiometry to improve our understanding of their structures and properties.

A systematic study of the phase diagrams of light actinide oxides was initiated at Los Alamos National Laboratory in 2017, and considerable progress has been made since then. Specifically, the melting curves of all the six stable stoichiometric uranium oxides, UO₂, U₄O₉, U₃O₇, U₂O₅, U₃O₈ and UO₃, were calculated using ab initio quantum molecular dynamics (QMD) simulations implemented with VASP (Vienna Ab initio Simulation Package), as well as one complete phase diagram in one case (UO₂). Some information on the location of solid-solid phase boundaries has been generated in a few other cases. This study is in progress, and its results will be reported elsewhere as they become available. In this work we present the systematics of the ambient melting temperatures of the stoichiometric uranium oxides. This systematics is aimed at bridging a gap between different versions of the ambient phase diagram of uranium-oxygen system. The determination of the phase-transition boundaries for the U-O system is of great importance in the nuclear industry. Knowledge of the melting transition in the nuclear fuel is particularly important in the analysis of hypothetical meltdown accidents, as it defines the structural limit of a combustible element. Moreover, due to a possible failure of the cladding during an accident, the fuel could come into contact with the coolant. If the latter is water, as in most reactors, strongly oxidizing conditions can be produced, under which the oxygen content of the fuel can be significantly increased. Hence, a precise knowledge of the fuel melting point dependence on the oxygen content is also of primary importance for the analysis of hypothetical mishaps.

Although a number of stable stoichiometric U-O compounds can be formed, only three of them are prevalent: UO₂, U₃O₈ and UO₃. U₃O₈ is the first isolatable layered oxide; all of the fluorite-type ones, U₄O₉, U₃O₇ and U₂O₅ oxidize and convert to U₃O₈ over time, or disproportionate to UO₂ + U₃O₈ depending on the conditions. U₃O₈ is commonly produced from UO₂ oxidation as a kinetically controlled product (Rousseau et al., 2006), via the reaction 3 UO₂ + O₂ → U₃O₈, but can also be formed from reduction of UO₃ at high T (Wen et al., 2013), via 6 UO₃ → 2 U₃O₈ + O₂. UO₃, the third prevalent uranium oxide, as a thermodynamically controlled product (Loopstra et al., 1977). Assuming that it is dominated by these three stable prevalent oxides, the high-T portion of the phase diagram of the U-O system was modeled by Babelot et al. (1986) and is shown in Figure 1. Specifically, the values of the ambient congruent melting points of U₃O₈ and UO₃ obtained in 5) are 2450 K and 1800 K, respectively. The only available corresponding experimental datum is that for U₃O₈ (as UO_2.66), 2010 ± 100 K (Manara et al., 2005), ∼ 20 % lower. All the other stoichiometric oxides are assumed to form a continuous UO_2+x solid solution which exhibits incongruent melting and a miscibility gap between the corresponding solidus and liquidus. This phase diagram was subsequently modified by (Manara et al., 2005) in view of the corresponding experimental results; see Figure 1.

In this work we present the results of ab initio quantum molecular dynamics (QMD) simulations on the ambient melting points of all the six stable stoichiometric uranium oxides. As already mentioned above, at high T the three fluorite-type oxides oxidize and decompose, which makes the corresponding high-T melting experiments very challenging. Another factor contributing to difficulty with high-T experiments is the sensitivity of the corresponding stoichiometry to oxygen content in the atmosphere and (most importantly) the lack of inert (unreactive) container material. For this reason, experimental melting results of 6) may not be accurate enough. In contrast, computer simulations are free from these difficulties: any stoichiometric (as well as non-stoichiometric) oxide can be simulated and studied as regards its thermal stability, melting curve (including the ambient melting point as the starting point of its melting curve), relative solid-solid phase stability, etc.

The ambient melting points the systematics of which is presented in this work should be understood as the starting (zero pressure) points of the corresponding melting curves. Each melting curve has been calculated using the Z method implemented with VASP. The Z method was developed in Refs. (Belonoshko et al., 2006; Belonoshko et al., 2008), and its practical use is described in detail in (10). In each case, the electronic structure of U was represented by [Xe 4f¹⁴ 5d¹⁰] 5f³ 6s² 6p⁶ 6d¹ 7s² (14 valence electrons), and that of O by [He] 2s² 2p⁴ (6 valence electrons). Localization of the f-electrons of U was achieved by using the DFT + U methodology (Dorado et al., 2009) following the Dudarev scheme (Dudarev et al., 1998a) for which only the difference in the values of Hubbard coefficients U and J matters rather than each of them individually. The use of DFT + U methodology to achieve f-electron localization prevents the appearance of delocalized f-electrons which would have resulted in non-integer uranium oxidation states, a nonphysical result considering recent experimental observations (Kvashnina et al., 2013). We used the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional (Perdew et al., 1996). Of all the GGA formulations that we tried for UO₂, PBE_sol + U predicts lattice parameters closest to experiment, while PBE + U gives the formation energy (E _form) closest to experiment. PBE (with no Hubbard coefficients) overestimates E _form by almost 1.5 eV, PBE + U underestimates by ∼ 0.2 eV, and PBE_sol + U underestimates by ∼ 0.6 eV. We have chosen to use PBE + U for our ab initio study, mainly because its E _form is the best match of the experimental one, which is the clear manifestation of its ability to accurately capture the UO₂ energetics. The values of U and J known in the literature are (in eV) 4.5 and 0.54, from the experimental x-ray photoemission spectra of actinide dioxides (Yamazaki and Kotani, 1991; Kotani and Yamazaki, 1992), or 4.5 and 0.51, from a DFT + U based theoretical analysis (Dorado et al., 2009). In addition to the theoretical study of UO₂ of Ref. (Dorado et al., 2009), the values of U = 4.5 and J = 0.51 were used in earlier theoretical studies of UO₂ in (Dudarev et al., 1998b; Vathonne et al., 2014) as well as those on U_1−y Pu _y O₂ mixed oxide (MOX) fuel (Dorado and Garcia, 2013; Njifon et al., 2018; Njifor and Torres, 2020). The above values of U and J imply that in the Dudarev scheme U–J ≈ 4.0 eV.

The proper choice of U is also important for the system under consideration to have the correct physical properties: the FM, AFM, or NM ground state, the values of the magnetic moment per atom and the band gap in agreement with experiment, the correct lattice constant (or density), etc. We used the occupation matrix control (OCM) method which allows one to study the physical properties of systems of either of the three magnetic states as a function of U. OCM for VASP was suggested by (Allen and Watson, 2014) and has become available in VASP starting with the version vasp5.3.5. (Our present study was carried out using vasp5.4.4.) Its use for both UO₂ and PuO₂ was very recently demonstrated by Chen and Kaltsoyannis (2022). Employing their approach, we studies the properties of UO₂ as a function of U using the Dudarev scheme (which is equivalent to U–J in the standard DFT + U formulation). In this case, imposing (Chen and Kaltsoyannis, 2022) i) the AF ground state, ii) magnetic moment of 1.8 μ _B , iii) the ground state lattice constant of 5.470 Å, and iv) the band gap of 2.25 ± 0.25 eV leads to (in eV) i) U ∼ > 3 , ii) 2.5 ∼ < U ∼ < 3 (U = 4 gives ∼ 1.9 μ B ), iii) 3 ∼ < U ∼ < 4 , and iv) 4 ∼ < U ∼ < 5 , so that, optimizing between the four criteria gives 3.5 and 4 as equally best overall choice of U for UO₂. In our study we take U = 4 eV taking into account the arguments of the previous paragraph. Since no (U, J) parameter sets are known for other uranium oxides, we use U − J = 4.0 eV for all the stoichiometric uranium oxides, to ensure both the transferability of the “effective” U-U and U-O inter-atomic interaction potentials and the ease of the reproducibility of our results by other researchers.

In each of the six cases of the stoichiometric uranium oxides considered here, we first prepared the suprecell that was used for melting simulations with the Z method. The supercell was prepared by relaxing the solid structure that the corresponding uranium oxide melts from, with the convergence of total energy to 10^–6 eV/atom and that of Hellmann-Feynman forces on each atom to 0.01 eV/Å. Size of the supercell varied from ∼ 300 atoms to ∼ 800 atoms; see below for more detail. In each of the six cases, the corresponding supercell was subject to a set of initial temperatures (T ₀) separated by an increment of 250 K and run with QMD in the NVE ensemble, for a total of up to 20,000 time steps of 1.5 fs each, i.e., up to 30 ps, to determine the melting temperature (T _m ). During this running time, in each of the six cases the system fully equilibrates. That is, if a QMD run continues past 30 ps, both T and P do not move away from the corresponding equilibrium values (which we take to be the outcome of this run) but fluctuate around them, as Figures 2–7 clearly demonstrate. Note that a larger UO₂ system of 960 atoms discussed in the end of this paper requires a little longer simulation time of 25,000 time steps (37.5 ps) for full equilibration. But the smaller and larger UO₂ systems are entirely consistent with each other as regards the corresponding equilibrium p = 0 values of T _m . Uncertainty of the value of T _m intrinsic to the Z method is therefore 125 K, half of the increment of T ₀ (Burakovsky et al., 2015), which constitutes 4%–7% of T _m (the largest uncertainty of 7% is for UO₃, and the lowest one of 4% for UO₂). Each of our six ambient melting points corresponds to a pressure in the interval (−1 GPa, 1 GPa), i.e., 0 ± 1 GPa. Assuming that the initial slope of the melting curve (dT _m /dP at p = 0) is ∼ < 100 K/GPa [the only known value of this slope for uranium oxides is the one for UO₂: (92.9 ± 17.0) K/GPa (Manara et al., 2010)], a P uncertainty of 1 GPa translates into a T _m uncertainty of δ T m ∼ d T m d P δ P ∼ < 100 K which is within the 125 K uncertainty of the method itself. This seems to be the case for both UO₂ and UO₃ for which the average values of P are, respectively, ∼–1 GPa and ∼ 1 GPa; see Figures 3, 7.

Thus, the results of our QMD simulations of the ambient T _m are expected to be quite accurate overall. Now we discuss these results in more detail for each of the six stoichiometric uranium oxides.

As Figures 2–7 clearly demonstrate, although the Z-method simulations of uranium oxides exhibit some basic features of those of monatomic substances, e.g., during melting T goes down while P goes up, unlike the case of a monatomic substance, the plots of the time evolution of T and P are not exact mirror images of each other (for a monatomic substance, they are exact mirror images of each other because they oscillate in anti-phase and their oscillations are of similar amplitudes since the total energy, ∼ P V + k _B  T, is conserved). Moreover, comparison of the corresponding pairs of the plots of the time evolution of T and P show that over some intervals of running time both T and P increase (or decrease). For example, in Figures 2, 3 for UO₂ both T and P increase over the 15,500–17,500 interval of running time. Let us dwell upon this issue in some more detail.

Upon the comletion of the superionic transition, the system consists of two subsystems: the fast one (“f”) composed of light (oxygen) ions, and the slow one (“s”) composed of heavy (uranium) ions. As such, the total energy of the system if E = E _f + E _s; E _f ∼ P _f V _f + k _B T _f, E _s ∼ P _s V _s + k _B T _s. In a general case, the temperature and pressure of such a complex system is T = t _f T _f + t _s T _s and P = p _f P _f + p _s P _s, where the weights t _f, t _s, p _f and p _s depend on the corresponding equations of state, the ionic mass ratio, etc. During the equilibration stage T _f → T _s → T and P _f → P _s → P. However, the corresponding time evolution of T and P is directly related to the time evolution of, respectively, T _f and T _s, and P _f and P _s which, in turn, depend on the corresponding weights t _f and t _s, and p _f and p _s. Thus, it is in principle possible that during the equilibration stage T and P go up or down in anti-phase, as well as that T and P increase or decrease simultaneously. But in our QMD simulations, the effect of a simultaneous T-P increase/decrease may be an artifact related to the size of the computational cell. Based on general grounds, we expect this effect to depend on the P-T balance for the fast subsystem which may be shifted towards P for a smaller system. Indeed, for a smaller supercell, the VASP periodic boundary conditions introduce additional “walls” inside the bulk of a system, which enhances P _f, in the assumption that the fast subsystem is approximated by a liquid or gas taking the entire system volume. With increasing size of the supercell, these additional “walls” will gradually disappear; correspondingly, the effect itself will diminish and eventually disappear altogether.

To test this hypothesis, we carried out another suite of QMD simulations of UO₂ in cubic fluorite structure choosing a larger 4 × 4 × 5 supercell of 960 atoms (320 U plus 640 O). Choosing a larger supercell is unfeasible, since for the 960-atom system the corresponding total number of valence electrons is N _e = 8320, on the verge of the current VASP capability of 10⁴ atoms; VASP simulations of systems with larger N _e become virtually endless. Figures 8, 9 demonstrate, respectively, the time evolution of T and P of this system. As clearly seen in these figures, i) the effect of a simultaneous T-P increase/decrease is now absent, and ii) the plots of T and P are virtually the mirror images of each other. In other words, increasing the size of the computational cell virtually eliminates the effects related to the appearance of a fast subsystem. Finally, we note that the numerical value of the ambient T _m for a larger system is virtually identical to that for a smaller system (i.e., ≈ 3100 K). That is, smaller systems considered in this work are fully converged as regards the corresponding values of T _m . In this respect, our present results on the ambient T _m s of uranium oxides should be considered as accurate, which was already pointed out in Section 3.

Another point, worth of dwelling upon, is the discrepancy between our theoretical results and some experimental results, e.g., those of (Manara et al., 2005), as well as the discrepancy between different experimental results themselves, e.g., those of (Manara et al., 2005; Latta and Fryxell, 1970), as seen in Figure 1. The possible explanation for this discrepancy that readily suggests itself is the oxidation of the original substance, in terms of the accommodation of additional oxygen atoms (interstitial defects) in the crystal lattice of the original structure during experiments which may influence their outcome. To characterize this effect quantitatively, a through theoretical study of the oxydation process is required which goes far beyond the scope of our work. Specifically, at low concentrations these defects take the form of single oxygen atoms, however, as more oxygen is added to the original lattice, the defects aggregate to form clusters, such that the knowledge of the clustering scheme is required for the quantitative characteristics of the discrepancy between the theoretical and experimental melting points.

Thorough combined experimental and theoretical studies of the factors that may be responsible for the observed discrepancies in the values of T _m in different experiments were recently conducted by two research groups (Efipano et al., 2020; Fouquet-Métivier et al., 2023). In (35) the melting behavior of americium-uranium MOX under different atmospheres was investigated, and it was concluded that uranium-rich samples melt at temperatures significantly lower (around 2700 K) when they are laser-heated in a strongly oxidizing atmosphere (compressed air) compared to the melting points (beyond 3000 K) registered for the same compositions in an inert environment (pressurised Ar). This behavior was interpreted on the basis of the strong oxidation of such samples in air, leading to lower T _m s in the corresponding experiments. The results of 6) shown in Figure 1 are obtained using laser heating in pressurized (up to 300 MPa) atmosphere aimed at suppressing as much as possible the evaporation from the heated surface. In (36) the phase diagram of U-Pu-O system was studied. For uranium-rich samples, strong interaction with tungsten crucibles, and some interaction with rhenium crucibles was detected which explains why lower T _m s are systematically obtained using thermal arrest technique used in older experiments. The results of (7) shown in Figure 1 are obtained using this technique. The influence of the atmosphere in the laser heating tests was also highlighted. Specifically, uranium-rich samples tend to oxidize in air, which leads to a progressive decrease of the corresponding T _m s, and those that result from melting in argon atmosphere were recommended as more reliable.

Here we briefly summarize the main results of this work.

We have calculated the ambient melting points of the six stable stoichiometric uranium oxides via ab initio QMD simulations using the Z method implemented with VASP. The six melting points are listed in Table 1, along with the corresponding ambient melting densities. Because experimental values of the ambient melting densities are not available in the literature yet, our theoretical results should serve as a guideline for future experiments.

Both the melting transition and the superionic transition preceding melting can be characterized quantitatively in each of the six cases of QMD simulations.

The ab initio ambient melting points of the three prevalent uranium oxides, UO₂, U₃O₈ and UO₃, obtained in our work are in good agreement with the corresponding experimental and/or theoretical data. Those of the remaining three, U₄O₉, U₃O₇ and U₂O₅, as well as that of U₃O₈, can be assigned a common value of 2465 ± 65 K. This is consistent with the 3-phase (horizontal) equilibrium line on the U-O phase diagram at 2450 K shown in Figure 1. This figure clearly demonstrates that a family of the six ambient melting points map out a solid-liquid transition boundary consistent with the high-temperature portion of the phase diagram of uranium-oxygen system suggested by Babelot et al. (1986) rather than that of Manara et al. (2005).

In view of our findings of the values of the congruent melting points of the three less prevalent stoichiometric uranium oxides, the uranium-oxygen phase digram shown in Figure 1 should be reconsidered appropriately. This will be discussed in more detail in one of our future publications.