Ambient Melting Behavior of Stoichiometric Uranium-Plutonium Mixed Oxide Fuel ^†

Abstract

Mixed oxides of uranium and plutonium (MOX) are currently considered as a reference fuel for the new generation of fast breeder reactors such as ASTRID. The key factor determining the performance and safety of a fuel such as MOX is its operational limits in the application environment which are closely related to the material’s structure and thermodynamic stability. They are in turn closely related to the ambient (zero pressure) melting point $\left(T_m\right)$; thus, $T_m$ is an important engineering parameter. Furthermore, PuO${}_2$ and UO${}_2$ are two endpoints of the phase diagram of MOX; therefore, their ambient $T_m$s are fundamental reference points. However, the current knowledge of the $T_m$ of MOX is limited and controversial as several studies available in the literature do not converge on the unique behavior of $T_m$ as a function of $x.$ Specifically, some studies produced $T_m$ as a monotonically decreasing function of x such that, with $T_m$ of UO${}_2$$(x=0)$ of 3150 K, $T_m$ of PuO${}_2$$(x=1)$ is ∼2650 K, while other studies resulted in $T_m$ having a local minimum at $0.5<x<1$ such that $T_m$ of PuO${}_2$ is ∼3000 K, so that the difference between the two values of $T_m$ is as high as 350 K. In this study, using the ab initio Z method implemented with the Vienna Ab Initio Simulation Package (VASP), we carry out a suite of quantum molecular dynamics simulations to obtain the ambient $T_m$ of MOX at several values of $x,$ $0<x<1,$ including the two end points $(x=0,$ $x=1).$ Our results agree with the behavior of $T_m$ of MOX as a function of x having a local minimum at $x=0.7$ and $T_m$ of PuO${}_2$ of 3050 K. Our study suggests potential ambient density–melting point systematics of MOX which may be useful in subsequent research on MOX such as its thermoelasticity modeling.

1. Introduction

1.1. Mixed Oxide Fuel (MOX)

Of all its alternatives, nuclear power seems to be one of the most reliable energy sources to achieve our goal of a world free of human-made CO${}_2$ emissions [1]. However, nuclear power has always been under scrutiny due to its potential impact on the environment, and the safety issues related to its use (this concern comes from accidents such as Chernobyl or Fukushima). The choice between the two possible nuclear fuels, uranium and mixed oxide (MOX), as the most adequate for the future of worldwide power has been under debate. MOX has the advantage of being less disposable as it produces less intermediate- and high-level waste. Furthermore, MOX is planned to be the fuel for the new generation of fast breeder reactors such as the Advanced Sodium Technological Reactor for Industrial Demonstration (ASTRID).

Uranium–plutonium MOX is the system U${}_{1-x}$Pu${}_x$O${}_{2-y}$, where x is the plutonium fraction; $y=0$ and $y>0$ correspond to stoichiometric and substoichiometric MOX, respectively. With low plutonium content ($x<0.1$), it is used as a nuclear fuel in several thermal neutron reactors around the world. MOX with a higher plutonium content is expected to be a favorable fuel for fast neutron reactors (FNRs). The constraint of an oxygen-to-metal (O/M, M = U + Pu) ratio less than or equal to 2 is chosen as a safety precaution to protect the steel cladding from corrosion during irradiation in FNRs, even though this substoichiometry also has negative consequences such as a reduced thermal conductivity or the two-phase composition of the MOX. In this work, we will consider stoichiometric MOX only.

1.2. Ambient Melting Behavior

Of all the physical properties of a material, melting behavior is a fundamental property closely related to its structure and thermodynamic stability, and therefore has always been a crucial research subject. The ambient (zero pressure) melting point ($T_m$) is also an important engineering parameter as it defines the operational limits of a material in its application environment. It becomes critical in nuclear engineering where the thermo-mechanical stability of a nuclear fuel element is a key factor determining fuel performance and safety. Furthermore, PuO${}_2$ and UO${}_2$ are two endpoints of the phase diagram of MOX; therefore, their ambient $T_m$s are fundamental reference points.

The current knowledge of $T_m$ of MOX is limited and controversial as the literature studies do not converge on the unique behavior of $T_m$ as a function of $x.$ Specifically, studies [2,3] produced $T_m$ as a monotonically decreasing function of x such that, with $T_m$ of UO${}_2$ $(x=0)$ of 3150 K, $T_m$ of PuO${}_2$ $(x=1)$ is ∼2650–2700 K, while studies [4,5,6,7] resulted in $T_m$ having a local minimum at $0.5<x<1$ such that $T_m$ of PuO${}_2$ is ∼3000–3050 K, so that the difference between the two values of $T_m$ is as high as 350 K.

This controversy in the melting behavior of MOX is directy related to the ambiguity in the value of the ambient melting point of PuO${}_2$ as well as the melting behavior of MOX with high Pu content.

Since the 1960s, several groups have reported measurements of $T_m$ of PuO${}_2$ using various experimental techniques. Most of them are summarized in a review by Carbajo et al. [2] who recommended a value of $2701\pm 35$ K based on measurements circa 1960s using the thermal arrest technique on tungsten-encapsulated samples. More recently, Guéneau et al. [3] recommended a value of 2660 K based on published data and the thermodynamic modeling of the Pu-O system. However, in the same year, Kato et al. presented their experimental study [4,5] in which the then-commonly-accepted value of $T_m$ for PuO${}_2$ was called into question, and a considerably higher value of ∼3000 K was proposed. While using essentially the same thermal arrest technique as in the 1960s, Kato et al. paid particular attention not only to maintaining the exact O/M ratio, as the previous investigators did, but also to the effect of sample–containment interactions. In this way, they could attribute lower values of $T_m$ in the previous studies to extensive interactions between PuO${}_2$ samples and tungsten, typical crucible material in this range of $T.$ The most recent experimental values of $T_m$ for PuO${}_2$ are $3017\pm 28$ K [6] and $3050\pm 59$ K [7].

The difficulties with the experimental determinations of the melting behavior of the stoichiometric MOX as well as the shortcomings of the experimental techniques used for these determinations are summarized in the recent review by Fouquet-Métivier et al. [8]. Specifically, (i) the interaction of the U${}_{1-x}$Pu${}_x$O${}_2$ samples of high Pu content with tugsten crucibles, and in some $x>0.5$ cases with rhenium crucibles, drive the corresponding $T_m$ systematically lower; (ii) Pu content affects the $T_m$ measurements such that laser heating which is currently the most widely used technique causes a non-ideal behavior of the UO${}_2$ and PuO${}_2$ partitions of MOX which results in a distinct minimum of $T_m$ as a function of x at $0.4<x<0.8;$ (iii) the measured values of $T_m$ are affected by the O/M ratio which varies in the atmosphere during the experiments, so that it is not clear whether the measured $T_m$ corresponds to an O/M of close to 2 in the stoichiometric case or if it is much less than 2.

Hence, the clarification of the behavior of $T_m$ as a function of x requires further study. Such study, based on first principles, is the subject of our work.

The reliability of the first-principle-based methods is widely accepted, specifically regarding the calculation of the melting curves and even complete phase diagrams including both melting curves (solid–liquid transition lines) and solid–solid phase transition boundaries. For example, comparison of the calculated melting curves of zinc [9], platinum [10] and osmium [11] to the corresponding experimental data of Refs. [9,10,12] reveals excellent agreement in each of the three cases. Similarly, complete first-principle-based phase diagrams of copper [13] and silver [14] are in good agreement with the experimental data of Refs. [15] and [16], respectively.

As regards the actinides and their compounds, a first-principle-based study of their melting curves or phase diagrams has not been undertaken yet. This is mainly related with the computational intensity of such studies and the availability of the adequate resource. Specifically, a 768-atom PuO${}_2$ system studied in our work, with the corresponding core–valence electron structure of both Pu and O (16 valence electrons of Pu and 6 of O), has a total of 7168 valence electrons (rather close to the current “unofficial” limit of 10,000 beyond which it is virtually impossible to run computer simulations to completion) which makes the simulations of such a system extremely time consuming. Thanks to the adequate computational resource at Los Alamos, we could complete these simulations within a reasonable amount of time of several months.

We now switch to the detailed discussion of the results of our first-principle-based study and the methods that we used.

2. Methods

In this paper, we calculate the ambient melting points of MOX at several values of $x,$ $0\le x\le 1,$ including the two end points, using the ab initio Z method based on quantum molecular dynamics (QMD) simulations implemented with the Vienna Ab Initio Simulation Package (VASP). The Z method is described in detail in [11]. Our approach is based on density functional theory (DFT) with local density approximation (LDA). We note that the standard implementation of DFT with LDA in the case of PuO${}_2$ predicts the ground state with either ferromagnetic (FM) or antiferromagnetic (AFM) order which is inconsistent with the nonmagnetic (NM) ground state observed in experiment. Although the disagreement may be due to strong correlations of f-electrons of Pu, more advanced computational schemes such as standard LDA + U or hybrid DFT which take into account the strong correlation effects of f-electrons are known to fail to generate nonmagnetic ground states [17]. Only the schemes modified by the inclusion of spin–orbit coupling in addition to the correlation effects have been shown to generate the correct ground state [17]. Hence, here we use such an LDA + U computational scheme with spin–orbit coupling.

2.1. LDA + U Methodology

The use of the LDA + U methodology is not specific to UO${}_2$, but is a standard approach to dealing with oxides in general. In particular, for UO${}_2,$ it is used to achieve f-electron localization which prevents the appearance of delocalized f-electrons which would have resulted in non-integer uranium oxidation states, a nonphysical result considering recent experimental observations [18].

For our theoretical study the electronic structure of U was represented by [Xe 4f${}^{14}$ 5d${}^{10}$] 5f${}^3$ 6s${}^2$ 6p${}^6$ 6d${}^1$ 7s${}^2$ (14 valence electrons), that of Pu by [Xe 4f${}^{14}$ 5d${}^{10}$] 5f${}^6$ 6s${}^2$ 6p${}^6$ 7s${}^2$ (16 valence electrons), and that of O by [He] 2s${}^2$ 2p${}^4$ (6 valence electrons). The LDA + U approach was implemented with the Dudarev scheme [19] for which only the difference in the values of Hubbard coefficients U and J matters rather than each of them individually. We used (in eV) $U=4.0$ and $J=0$ for both U and Pu which seems to be the widely adopted choice. See, e.g., [17] or [20], in which the following sets of $(U,J)$ are chosen for U and Pu, respectively: $(4.5,0.5)$ and $(4.7,0.7),$ so that in either case $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 [21] and has become available in VASP starting with the version vasp5.3.5. Its use for both UO${}_2$ and PuO${}_2$ was very recently demonstrated by Chen and Kaltsoyannis [22]. Using their approach, we studied the properties of both UO${}_2$ and PuO${}_2$ as a function of $U.$ In the case of UO${}_2,$ imposing [22] (i) the AF ground state, (ii) magnetic moment of 1.8 ${\mu }_B,$ (iii) the ground state lattice constant of 5.470 Å, and (iv) the band gap of $2.25\pm 0.25$ eV leads to (in eV) (i) $U\stackrel{>}{\sim }3,$ (ii) $2.5\stackrel{<}{\sim }U\stackrel{<}{\sim }3$ ($U=4$ gives $\sim 1.9 {\mu }_B$), (iii) $3\stackrel{<}{\sim }U\stackrel{<}{\sim }4,$ and (iv) $4\stackrel{<}{\sim }U\stackrel{<}{\sim }5,$ so that, optimizing between the four criteria gives 3.5 and 4 as equally best overall choice of U for UO${}_2.$ Similarly, in the case of PuO${}_2,$ imposing [22] (i) the NM ground state, (ii) the ground state lattice constant of 5.395 Å, and (iii) the band gap of 1.8 eV leads to (in eV) (i) $3.5\stackrel{<}{\sim }U\stackrel{<}{\sim }4,$ (ii) $4\stackrel{<}{\sim }U\stackrel{<}{\sim }5,$ and (iii) $4.5\stackrel{<}{\sim }U\stackrel{<}{\sim }5,$ so that, optimizing between the three criteria gives 4 and 4.5 as equally best overall choice of U for PuO${}_2.$ In our study, we take $U=4$ eV for both UO${}_2$ and PuO${}_2$; besides being one of the best overall choices, in view of the above considerations, it will ensure both the transferability of the “effective” inter-atomic interaction potentials within the UO${}_2$ and PuO${}_2$ components of MOX of a given composition and the ease of the reproducibility of our results by other researchers.

2.2. Melting Simulations

Both UO${}_2$ and PuO${}_2$ maintain their ambient cubic fluorite structure (Fm$\overline{3}$m) in the entire range of T up to the corresponding melting points. Specifically, for UO${}_2$ it is $\rho = 10.97$ g/cm${}^3,$ $a=b=c=5.47$  Å [23] while for PuO${}_2$ $\rho =11.7$ g/cm${}^3,$ $a=b=c=5.39$ Å [24]. It is described in terms of a 12-atom unit cell containing four uranium atoms in face-centered cubic positions and eight oxygen atoms filling the tetrahedral sites; see Figure 1.

For our melting simulations of MOX, we used a supercell of 324 atoms $(3\times 3\times 3).$ This supercell size ensures the convergence of total energy of each QMD simulation to ${10}^{-6}$ eV/atom and that of Hellmann–Feynman forces on each atom to 0.01 eV/Å. The mixed UO${}_2$-PuO${}_2$ system is modeled as a substitution alloy in which some U atoms (at randomly chosen lattice sites) are replaced by Pu ones according to the corresponding PuO${}_2$ content. For instance, in the case of U${}_{0.75}$Pu${}_{0.25}$O${}_2$ 81 uranium atoms chosen at random are replaced with 81 plutonium atoms, etc.

In the melting simulations, the corresponding supercell is subject to a set of initial temperatures $\left(T_0\right)$ separated by an increment of 250 K and run with QMD in the $NVE$ ensemble, for a total of up to 25,000 time steps of 1 fs each for $x=0$ and $x=1$ (both pure UO${}_2$ and PuO${}_2$) and up to 15,000 time steps of 1.5 fs each for $x\ne 1$, i.e., up to 25 ps in either case, to determine $T_m$ at the corresponding melting pressure $\left(P_m\right).$ As the system equilibrates upon the completion of the melting process, the values of $T_m$ and $P_m$ are determined from the corresponding running averages which are shown as solid lines in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8. Note that in any molecular dynamics simulations, either classical or ab initio quantum ones, similar to those discussed here, both T and P oscillate during a run. Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 detail such oscillations. Their magnitude depends on both the complexity of the substance under investigation and the system size: the more complex the substance and the smaller the system size the larger the amplitude of the oscillations, and vice versa. In our case the energy transfer between the lighter oxygen and heavier actinide subsystems during the QMD run results in larger-than-usual amplitude of the oscillations, up to ∼500 K for T and ∼5 GPa for $P.$ However, we emphasize that it is the running averages only that matter not the oscillations around them which are just a molecular dynamics artifact. Over the running time as long as 20,000 time steps, uncertainty related to these oscillations is only a fraction of a percent, so that the corresponding running averages exhibit virtually the exact values of T and $P.$ Hence, the running averages during the latest portion of the run yield the values of $T_m$ and $P_m.$ Uncertainty of the value of $T_m$ intrinsic to the Z method is 125 K, half of the increment of $T_0$ [11], which turns out to constitute ∼2% of $T_m.$ Uncertainty of the value of $P_m$ in our simulations is ∼0.5 GPa. Assuming that the initial slope of the melting curve $(d{T}_m/dP$ at P = 0) is $\stackrel{<}{\sim }100$ K/GPa (see below), a $P_m$ uncertainty of 0.5 GPa translates into a $T_m$ uncertainty of $\delta T_m\sim \frac{d{T}_m}{dP}\delta P\stackrel{<}{\sim }50$ K. Therefore, the combined uncertainty of $T_m$ in our QMD simulations is $\approx 135$ K, which is within 5% of $T_m$ in each case. Thus, our results on the ambient $T_m$ are expected to be quite accurate overall.

3. Results

Our results of the QMD simulations of the ambient points of MOX are summarized in

Table 1. The ab initio melting points of U${}_{1-x}$Pu${}_x$O${}_2$, $(P_m,T_m\pm \mathsf{\Delta }{T}_m),$ obtained from the Z method implemented with VASP, along with the values of lattice constants and densities at the corresponding $T_m$.

x	${\mathit{T}}_{\mathit{m}}$ (K)	$\mathsf{\Delta }{\mathit{T}}_{\mathit{m}}$ (K)	Lattice Constant (Å)	Density (g/cm${}^3$)
0	3150	135	5.700	9.685
0.25	2910	135	5.653	9.984
0.5	2870	135	5.624	10.20
0.7	2670	135	5.585	10.46
0.9	2750	135	5.572	10.58
1	3050	135	5.590	10.50

.

Figure 2 and Figure 3 show the (ambient) melting run [11] of Fm$\overline{3}$m-UO${}_2$ at a density of 9.685 g/cc (lattice constant 5.700 Å) which corresponds to the first line in Table 1. As seen in Figure 2 and Figure 3, at the beginning of the run, the system equilibrates as a UO${}_2$ solid (i.e., a quasi-static cation/anyon system). After ∼6000 time steps (6 ps) T goes down and P goes up (since in the $NVE$ ensemble the total energy, $E\sim k_{B}T+PV,$ is conserved). This is a signature of a superionic transition. Due to a 15-fold difference in the atomic masses of U and O, the oxygen sublattice becomes less stable than the plutonium one, and at sufficiently high T it disorders first, such that the anions (O${}^{-}$) start flowing through the ordered structure of the cations (U${}^{+}$). Such a (superionic) phase transition accompanied by a rapid increases in ionic conductivity has been observed in many diatomic systems. We note that for UO${}_2$, this transition was discovered by Dworkin and Bredig [25] from a curvature change of the T dependence of the enthalpy.

Thus, the first drop in T (increase in P) corresponds to the activation of the oxygen flow. This process takes ∼3000 time steps (∼3 ps). Then the system of quasi-static cations and mobile anyons equilibrates, and the second drop in T (increase in P) occurs after a total of ∼13,000 time steps (∼13 ps). This second drop in T is associated with the disordering of the uranium sublattice, i.e., a true melting transition. The melting process takes ∼5000 time steps (5 ps). The emerging liquid equilibrates at $(P,T)\approx (0,3150\mathrm{K})$.

Figure 2. Time evolution of temperature in the Fm$\overline{3}$m-UO${}_2$ melting run at $\rho =9.685$ g/cc. The initial temperature is 8500 K.

Figure 2. Time evolution of temperature in the Fm$\overline{3}$m-UO${}_2$ melting run at $\rho =9.685$ g/cc. The initial temperature is 8500 K.

Figure 3. The same as in Figure 4 for the time evolution of pressure (in kbar; 10 kbar = 1 GPa).

Figure 3. The same as in Figure 4 for the time evolution of pressure (in kbar; 10 kbar = 1 GPa).

Figure 4 and its P counterpart shown in Figure 9 demonstrate the (ambient) melting run of Fm$\overline{3}$m-PuO${}_2$ at a density of 10.50 g/cc (lattice constant 5.590 Å). We have obtained its three melting points close to $P=0$ based on which its ambient melting point was calculated (which corresponds to the last line in Table 1); see the next section. Just as for UO${}_2,$ the occurrence of the superionic phase transition is very clearly seen. For PuO${}_2$ it was discovered by Chroneous et al. [26] from the analysis of oxygen self-diffusion coefficients. Superionic behavior of PuO${}_2$ was modeled in [27] using classical molecular dynamics.

Figure 4. Time evolution of temperature in the Fm$\overline{3}$m-PuO${}_2$ melting run at $\rho =10.49$ g/cc. The initial temperature is 8000 K.

Figure 4. Time evolution of temperature in the Fm$\overline{3}$m-PuO${}_2$ melting run at $\rho =10.49$ g/cc. The initial temperature is 8000 K.

Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 show the (ambient) melting runs of Fm$\overline{3}$m-U${}_{1-x}$Pu${}_x$O${}_2,$ $x=0.25,$ 0.5, 0.7 and 0.9. The superionic phase transition occurs in each case, although it is not seen as clearly as for both pure UO${}_2$ and PuO${}_2$.

Figure 5. The same as in Figure 4 for U${}_{0.75}$Pu${}_{0.25}$O${}_2$.

Figure 5. The same as in Figure 4 for U${}_{0.75}$Pu${}_{0.25}$O${}_2$.

Figure 6. The same as in Figure 4 for U${}_{0.5}$Pu${}_{0.5}$O${}_2$.

Figure 6. The same as in Figure 4 for U${}_{0.5}$Pu${}_{0.5}$O${}_2$.

Figure 7. The same as in Figure 4 for U${}_{0.3}$Pu${}_{0.7}$O${}_2$.

Figure 7. The same as in Figure 4 for U${}_{0.3}$Pu${}_{0.7}$O${}_2$.

Figure 8. The same as in Figure 4 for U${}_{0.1}$Pu${}_{0.9}$O${}_2$.

Figure 8. The same as in Figure 4 for U${}_{0.1}$Pu${}_{0.9}$O${}_2$.

Figure 9. The same as in Figure 3 for U${}_{1-x}$Pu${}_x$O${}_2,$ $0<x<1.$ Only the running averages are shown in either case.

Figure 9. The same as in Figure 3 for U${}_{1-x}$Pu${}_x$O${}_2,$ $0<x<1.$ Only the running averages are shown in either case.

To conclude this section, we note that, since in the $NVE$ ensemble the total energy, $E\sim k_{B}T+PV,$ is conserved, one should expect the plots of the time evolution of T and P during a QMD run to be exact mirror images of each other (i.e., when T drops, P increases, and vice versa). This is indeed the case for each of the elemental solids studied using the Z method, e.g., Refs. [11,13,14], which all offer the corresponding plots. In our case of stoichiometric MOX, for smaller 324-atom systems such P-T, mirror reflection is not clearly observable, although their time evolution during the QMD run follows the right tendency—when T drops, P increases. This mirror-reflection behavior becomes almost exact in the case of a larger 768-atom system, as Figure 10, Figure 11, Figure 12 and Figure 13 clearly demonstrate. The observed violations of the P-T “mirror reflectivity” for a smaller system may be related to the 15-fold difference in the atomic masses of M (=U + Pu) and O and the way VASP calculates T and $P;$ however, the detailed analysis of the reasons for these violations goes far beyond the scope of this work. We note that, in any event, the results on the 324-atom systems are reliable since these systems are fully converged in regard to the system size, as discussed in more detail below.

The following section offers the detailed discussion of our results.

4. Discussion

We have obtained six ab initio melting points of MOX, including the two endpoints, which are listed in Table 1. The best polynomial fit to these six $T_m$s, Equation (2) below, results in the melting line of MOX which is one of the fundamental features of its phase diagram. Our results appear to agree with Refs. [4,5,6,7] on the behavior of $T_m$ of MOX as a function of x having a local minimum at $x=0.7$ and $T_m$ of PuO${}_2$ of 3050 K. In what follows, we dwell upon several key issues related to our study: (i) the validity of our results in regard to the convergence with respect to the system size; (ii) the initial slope of the melting curve of PuO${}_2,$ one of the end members of MOX, in comparison to that of UO${}_2,$ the other end member for which it is known experimentally; (iii) the ambient phase diagram of stoichiometric MOX which is determined by its melting line based on the size $T_m$s determined in this work; and (iv) potential ambient density–melting point systematics of MOX which may be useful for its thermoelasticity modeling.

4.1. Energy Convergence of Our Results with Respect to the System Size

Although full energy convergence is achieved in each QMD simulation, the accuracy of our results in regard to the system size must also be ensured. We checked that the full energy convergence of our results with respect to the system size is achieved for both UO${}_2$ and PuO${}_2.$ Since all the MOX systems were studied using supercells of exactly the same size as that of both UO${}_2$ and PuO${}_2,$ namely $(1-x)$ UO${}_2$ + x PuO${}_2,$ full energy convergence must have been achieved in each case. Specifically, in the case of UO${}_2,$ a supercell of 324 atoms $(3\times 3\times 3)$ is enough to achieve the energy convergence of up to 1 meV/atom, or uncertainty of ∼10 K as regards the value of $T_m,$ which is an order of magnitude less than the uncertainty of the Z method itself. Indeed, direct comparison of total energies as the outcomes of the 324- and 768-atom $(4\times 4\times 4)$ runs, for both UO${}_2$ and PuO${}_2,$ clearly shows that the 1 meV/atom energy convergence was already achieved. Figure 10 and Figure 11 demonstrate the (ambient) melting run of the 768-atom UO${}_2$ system. In this case $T_m\approx 3200$ K, which is consistent with its 324-atom counterpart of ≈3150 K within the uncertainty of the method itself. Similarly, Figure 12 and Figure 13 demonstrate the (ambient) melting run of the 768-atom PuO${}_2$ system. In this case $T_m\approx 3100$ K, which is also consistent with its 324-atom counterpart of ≈3050 K within the uncertainty of the method itself. It is clearly seen in these figures that $P_m$ has also converged, just as $T_m$.

Figure 10. The same as in Figure 2 for the 768-atom UO${}_2$ system.

Figure 10. The same as in Figure 2 for the 768-atom UO${}_2$ system.

Figure 11. The same as in Figure 3 for the 768-atom UO${}_2$ system.

Figure 11. The same as in Figure 3 for the 768-atom UO${}_2$ system.

Figure 12. The same as in Figure 4 for the 768-atom PuO${}_2$ system.

Figure 12. The same as in Figure 4 for the 768-atom PuO${}_2$ system.

Figure 13. The same as in Figure 3 for the 768-atom PuO${}_2$ system. Only the running average is shown.

Figure 13. The same as in Figure 3 for the 768-atom PuO${}_2$ system. Only the running average is shown.

4.2. Initial Slope of the Melting Curve of PuO₂

The initial melting slope is an important thermodynamic characteristic of a substance because (i) it relates the melting point with the latent heat of fusion, $L_m,$ and volume change at melt, $\mathsf{\Delta }{V}_m,$ via the Clausius–Clapeyron formula $d{T}_m/dP=T_m\mathsf{\Delta }{V}_m/L_m,$ hence it is important for the validation of the thermodynamic modeling using, e.g., Calphad’s approach, and (ii) it is related to the parameters a and b of the Simon–Glatzel (SG) melting equation $T_m\left(P\right)=T_m\left(0\right){(1+P/a)}^b,$ via $d{T}_m{/dP|}_{P=0}=T_m\left(0\right)b/a,$ hence it is important for the validation of the melting curve as a whole and of the corresponding SG equation.

Theoretical calculations and thermal expansion estimates reveal that at T∼3000 K, the lattice constant of the fluorite structure of PuO${}_2$ is ∼5.6 Å [24]. We therefore ran melting simulations for a set of three supercells having a lattice constant of ∼5.6 Å to determine the low-P segment of the melting point of PuO${}_2$ in order to obtain its initial slope to estimate the uncertainty in T associated with a potential uncertainty in $P.$ Our results are summarized in

Table 2. The three ab initio melting points of PuO${}_2$, $(P_m,T_m\pm \mathsf{\Delta }{T}_m),$ obtained from the Z method implemented with VASP.

Lattice Constant (Å)	Density (g/cm${}^3$)	${\mathit{P}}_{\mathit{m}}$ (GPa)	${\mathit{T}}_{\mathit{m}}$ (K)	$\mathsf{\Delta }{\mathit{T}}_{\mathit{m}}$ (K)
5.597	10.46	−1.8	2780	125
5.590	10.50	−0.6	3010	125
5.587	10.51	0.8	3090	125

.

The three melting points of PuO${}_2$ are shown in Figure 14. Taking into account the corresponding uncertainties, the resulting analytic form that describes them best is

$$T_m\left(P\right)=3045.9+92.5P-25.9P^2,$$

(1)

where P stands for pressure, which can be considered as the low-P ($-2\mathrm{GPa}\le P\le 1\mathrm{GPa}$) melting curve of PuO${}_2.$ Thus, the ab initio ambient melting point of plutonim oxide that results from our QMD simulations is $T_m\approx 3046$ K, which is in excellent agreement with the most recent experimental studies: $3017\pm 28$ K [6] and $3050\pm 59$ K [7]. The initial slope of the PuO${}_2$ melting curve given by Equation (1), 92.5 K/GPa, is consistent with the known value of this slope for UO${}_2,$ the other endpoint of MOX: $92.9\pm 17.0$ K/GPa [28]. We note that the ambient melting points of PuO${}_2$ and UO${}_2$ appear to be close to each other as well: ∼3050 K and ∼3150 K (e.g., 3147 K [28]), respectively.

Interpolating the corresponding values in Table 2 gives the following values of the lattice constant and density at the ambient melting point of PuO${}_2$ of 3050 K which we adopt based on our QMD simulations: 5.590 Å and 10.50 g/cm${}^3$. All these values correspond to the last line in Table 1.

4.3. The Ambient Phase Diagram of Stoichiometric MOX

Figure 15 compares our theoretical results on the ambient melting points of U${}_{1-x}$Pu${}_x$O${}_2$ at $x=0,$ 0.25, 0.5, 0.7, 0.9 and 1, summarized in Table 1, to the experimental data from the literature. In this respect, Figure 15 represents the current knowledge of the ambient phase diagram of stoichiometric MOX.

We note excellent agreement betwen our melting line, as the best fit to our melting data from Table 1,

$$T_m\left(x\right)=3147.4-142.06x-2334.8x^2+2347.8x^3,$$

(2)

(shown as black line in Figure 15) and the experimental data of Kato et al. [4] (shown as green bullets), as well as good agreement with the results of Lyon and Baily [29] (shown as empty squares) up to x∼0.8. Other experimental data shown in Figure 15 come from [7,30,31,32]. We also note that for MOX with high PuO${}_2$ content $(x\stackrel{<}{\sim }1)$, disagreement between our QMD results and those of Refs. [7,30] is beyond three standard deviations; hence, additional research is needed towards the resolution of this discrepancy.

4.4. Potential Ambient Density–Melting Point Systematics of MOX

Our QMD results suggest potential ambient density–melting point systematics of MOX which may be useful in subsequent research on MOX such as its thermoelasticity modeling. Figure 16 shows both ${\rho }_m,$ the density of MOX at its ambient $T_m,$ and its ambient $T_m$ as a function of the PuO${}_2$ content.

It is clearly seen that ${\rho }_m$ and $T_m$ are sort of “mirror images” of each other, in the sense that when one increases, the other decreases, and vice versa. Taking into account uncertainties in $T_m,$ the following relation is found:

$${\rho }_m·T_m^{7/6}\approx \mathrm{const}=(1.12\pm 0.08){10}^5.$$

(3)

Thus, ${\rho }_m·T_m^{7/6}$ is a constant to within ∼93% accuracy. This relation can be used for, e.g., thermoelasticity modeling. The dislocation mediated melting model [33] suggests that $G_m/\left({\rho }_m{T}_m\right)\approx \mathrm{const}$ along the melting curve [34], where $G_m$ is the shear modulus at $T_m.$ Hence, in view of (3), one should expect that for MOX $G_m·T_m^{1/6}\approx \mathrm{const},$ which can be used as a constraint for modeling $G_m=G_m\left(\rho \right)$ as well as $G=G(\rho ,T)$ in the Preston–Wallace thermoelasticity model [35] $G(\rho ,T)=G(\rho ,T$=$0)\left(1-\alpha \frac{T}{T_m\left(\rho \right)}\right),$$\alpha \approx \mathrm{const}.$

5. Conclusions

In this study, we obtained the ambient melting points of MOX at several values of $x,$ $0<x<1,$ including the two end points $(x=0,$ $x=1),$ using the ab initio Z method based on QMD implemented with VASP, and compared them to experimental and independent theoretical results available in the literature. This comparison is shown in Figure 15 which demonstrates very good agreement with the literature results and ours. Specifically, our results appear to agree with Refs. [4,5,6,7] on the behavior of $T_m$ of MOX as a function of x having a local minimum at $x=0.7$ and $T_m$ of PuO${}_2$ of 3050 K. Figure 15 represents the current knowledge of the ambient phase diagram of stoichiometric MOX; this knowledge will be advanced in subsequent studies on the subject.

As Figure 16 clearly shows, the distinct maximum of ${\rho }_m={\rho }_m\left(x\right)$ corresponds to $x=0.9,$ while the distinct minimum of $T_m=T_m\left(x\right)$ to $x=0.7.$ We note that a very recent study of the MOX phase diagram using Calphad’s methodology [8] suggests the phase diagram of MOX (Figure 5 of [8]) which is very similar to ours, Figure 15, for which the melting line $T_m=T_m\left(x\right)$ also has a distinct minimum at $x=0.7$ but which is in better agreement with the experimental results by Böhler et al. [7] rather than those of both Kato et al. [4] and Lyon and Baily [29] that our MOX phase diagram is. This is why their value of $T_m$ at $x=0.7$ is ≈2975 K, about 300 K above ours. One additional source of uncertainty of experimental measurements may be the variation of the O/M ratio during the melting of MOX in the atmosphere which may contribute as much as $\sim 100$ K to the error in the exeprimental $T_m$ [36]. In any event, the above two values of $T_m$ for MOX with $x=0.7$ cannot be reconciled with each other within the uncertainties of our method itself; thus, additional studies are necessary to elucidate the nature of this discrepancy and to shed more light on the phase diagram of MOX.

Our study suggests potential ambient density–melting point systematics of MOX which may be useful in subsequent research on MOX such as its thermoelasticity modeling.