Of all the physical properties of a material, melting behavior is a fundamental property closely related to its structure and thermodynamic stability; it has thus always been a crucial subject of research. 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 in determining fuel performance and safety. Moreover, PuO 2 and UO 2 are two endpoints of the phase diagram of MOX, so their ambient T m s are fundamental reference points.

The current knowledge of the T m of MOX is limited, and the literature does not converge on the unique behavior of T m as a function of x . Carbajo et al. (2001) and Guéneau et al. (2008) 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. However, the studies of Kato et al. (2008a), Kato et al. (2008b), De Bruyckner et al. (2010), and Böhler et al. (2014) 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. Although the melting behavior of uranium–plutonium MOX has been experimentally addressed in many recent studies, it still lacks a sound theoretical foundation. In particular, the melting behavior of MOX has not yet been adequately modeled based on general thermodynamics principles or using an equation of state (EOS) approach.

This uncertainty in the melting behavior of MOX is directly 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 research groups have reported measurements of T m of PuO 2 using various experimental techniques. Most are summarized in a review by Carbajo et al. (2001), who recommended a value of 2701 ± 35 K based on measurements circa 1960s using the thermal arrest technique on tungsten-encapsulated samples. More recently, Guéneau et al. (2008) recommended a value of 2,660 K based on published data and the thermodynamic modeling of the Pu–O system. However, in the same year Kato and colleagues presented their experimental findings (Kato et al., 2008a; Kato et al., 2008b) which called into question the then commonly-accepted value of T m for PuO 2 and proposed a considerably higher value of ∼ 3000 K. While using essentially the same thermal arrest technique as in the 1960s, Kato and colleagues paid particular attention to not only maintaining the exact O/M ratio as previous research had done but also to the effect of sample–crucible interactions. This way, they could attribute lower values of T m in previous studies to extensive interactions between PuO 2 samples and tungsten—typical crucible material in this range of T . The shortcomings thus indicated have been properly taken into account in the most recent experimental studies on the melting of PuO 2 . These studies were carried out using a containerless laser-heating technique to produce values of T m for PuO 2 : 3017 ± 28 K (De Bruyckner et al., 2010) and 3050 ± 59 K (Böhler et al., 2014). The most recent theoretical value of T m for PuO 2 is 3046 ± 135 K (Burakovsky et al., 2023).

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 by Fouquet-Métivier et al. (2023) as follows. (i) The interaction of the U 1 − x Pu x O 2 samples of high Pu content with tungsten crucibles, and in some x > 0.5 cases with rhenium crucibles, drives 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 non-ideal behavior of the UO 2 and PuO 2 partitions of MOX, resulting 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 O/M close to 2 of 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. Here, we present a theoretical model for the melting curve (liquidus) of a mixture and apply it to MOX which is considered a mixture of pure UO 2 and PuO 2 .

We consider the case of ideal mixing, where the constituents of the mixture do not effectively interact with each other, which would otherwise result in the volume of the mixture being different from the sum of the volumes of its constituents. We assume that the mixture is of the form A 1 − x B x and 0 ≤ x ≤ 1 and that no stoichiometric A n B m (both m and n are integers ≥ 1 ) compound exists; otherwise, additional arguments should be invoked regarding the enthalpy of formation of such a compound, which will result in a modification of the analytic form of the melting curve which we now derive. Thus, any eutectic is neglected. Our derivation is based on the assumption that both of the melting curves of pure constituents are known in the analytic form of the pressure ( P ) dependence of the melting point: T m = T m ( P ) . This form can be either a simple polynomial fit, such as a quadratic or cubic, or a more sophisticated Simon–Glatzel form T m ( P ) = T m ( 0 ) ( 1 + P / a ) b , where T m ( 0 ) is the ambient T m and a , b = c o n s t , etc. However, nothing else is known about the mixture, except the two melting curves; for example, their EOSs are not available. Thus, the model will predict the melting curve of a mixture based on the melting curves of its components only. Finally, we assume that both constituents of the mixture are molten and are in a state of thermal equilibrium with each other at a common temperature T m so that the two partial pressures are such that T m 1 P 1 = T m 2 P 2 = T m . (1) In Equation 1 T m 1 ( P ) and T m 2 ( P ) are the melting curves of constituents 1 and 2, respectively, and P 1 and P 2 are the corresponding partial pressures at T = T m . We thus model the system’s liquidus, which is the true melting line along which all the constituents are molten; our approach generally applies to systems, the constituents of which are not completely miscible so that the system’s phase diagram may have a miscibility gap. Thus, in what follows, we use the subscript ℓ instead of m in order to associate the melting point with liquidus and not confuse it with solidus.

First, we develop the appropriate mixing rules. We begin with the EOS. If the cold ( T = 0 ) EOSs of constituents 1 and 2 are, respectively, P 1 ( ρ ) and P 2 ( ρ ) , their finite- T counterparts can be written as P 1 ρ , T = P 1 ρ , 0 + α 1 B T , 1 T , P 2 ρ , T = P 2 ρ , 0 + α 2 B T , 2 T , (2) where α and B T are, respectively, the thermal expansion coefficient and isothermal bulk modulus at temperature T . This form of the thermal EOS does not explicitly take into account the T -dependence of the bulk modulus, and/or the P - (or T -) dependence of the thermal expansivity, and is therefore approximate. Indeed, since P T , V = P T 0 , V + ∫ T 0 T α B T d T ,

Equation 2 results from the above relation (with T 0 = 0 ) provided that α ⋅ B T ≈ const—a potential weak T -dependence of B —is complemented by a similar weak dependence of α to keep their product (roughly) T -independent. As our previous theoretical studies reveal, thermal EOS of the form of Equation 2 holds for many substances, such as copper (Baty et al., 2021a), silver (Baty et al., 2021b), palladium (Baty et al., 2024), and body-centered cubic bismuth (Burakovsky et al., 2024a), among others. In fact, along the corresponding melting curves, thermal EOS of the form Equation 2 is virtually exact.

We use the pressure mixing rule (law of additive volumes, or the Amagat–Leduc model) which requires that the pressures of the components be equal at a chosen mixture composition, total volume, and temperature. This pressure mixing rule thus describes the pressure equilibrium of the mixture. The set of the equations that describe the mixture at pressure equilibrium are P 1 ρ 1 , T = P 2 ρ 2 , T (3) and 1 − x ρ 1 + x ρ 2 = 1 ρ . (4) where ρ 1 and ρ 2 are, respectively, the densities of constituents 1 and 2, ρ is the density of the mixture, and x is the (fractional) mass percentage of constituent 2 (without any loss of generality we consider constituent 1 as a host and constituent 2 as a dopant): x = M 2 / ( M 1 + M 2 ) , 1 − x = M 1 / ( M 1 + M 2 ) . Equation 4 is equivalent to M 1 / ρ 1 + M 2 / ρ 2 = ( M 1 + M 2 ) / ρ , which is the total volume of the mixture being the sum of the volumes of its constituents. It follows from the above relation that ρ = ρ 1 ρ 2 x ρ 1 + 1 − x ρ 2 = 1 − x ρ 1 + x ρ 2 1 + 1 − x x ρ 1 − ρ 2 2 ρ 1 ρ 2 . Hence, provided that | ρ 1 − ρ 2 | ≪ ρ 1 , ρ 2 , the density of the mixture is quasi-additive: ρ ≈ ( 1 − x ) ρ 1 + x ρ 2 . This is the case for MOX since the values of the densities of pure UO 2 and P u O 2 , ≈ 10.97 and 11.46 g/ c m 3 , respectively (Carbajo et al., 2001), are indeed very close, and the density of U 1 − x Pu x O 2 ± y is described by ρ ( x , y ) = ( 10.97 + 0.49 x ) ( 1 + 0.0039 y ) (Carbajo et al., 2001). Here x can be either mass or atomic (or volume) percentage discussed below because, for MOX, the two are essentially identical.

Note that in Equation 3 the value of T cannot be arbitrarily high—that is, the range of T for the applicability of the pressure mixing rule is limited. Indeed, at the liquidus temperature T ell all the constituents of the mixture are molten; hence, each individual P is determined by the corresponding melting curve T m = T m ( P ) such that, for constituent i , T ell = T m i ( P i ( T ell ) ) . Since different constituents have different melting curves T m i = T m i ( P ) , the corresponding values of P i ( T ell ) are different as well. Then, at T = T ell the value of P of the mixture is generally not equal to any of the P i ( T ell ) s and is given by a mixing rule different from Equations 3 and 4. This new mixing rule is discussed in more detail in the next section, taking into account that at the liquidus point, the mixture is at temperature rather than pressure equilibrium.

By switching from mass percentage x to (atomic) volume percentage X , via X = x x + 1 − x m 2 m 1 , 1 − X = 1 − x 1 − x + x m 1 m 2 . (5) In Equation 5 m 1 and m 2 are the atomic masses, respectively (hence x = X X + ( 1 − X ) m 1 m 2 and 1 − x = 1 − X 1 − X + X m 2 m 1 ) , Equation 4 converts to v = 1 − X v 1 + X v 2 , (6) where v , v 1 , and v 2 are the atomic volumes of the mixture and its constituent. Here, we define the atomic mass of the mixture as m = ( 1 − X ) m 1 + X m 2 . Equation 6 represents Zen’s law of the additivity of the atomic volume of a mixture (Zen, 1956) (analogous to the quasi-additivity of the density of the mixture discussed above). We note that in the case of real mixing, formulas for the total volume V = V 1 + V 2 (resulting from Equation 4), or its atomic-volume analog Equation 6 must be replaced with, respectively, V = V ̃ 1 + V ̃ 2 and v = ( 1 − X ) v ̃ 1 + X v ̃ 2 . Here, “ ∼ ” indicates the partial (effective) volume in the mixture, which may be larger or smaller than that in the ideal (non-interactive) case depending on whether the mixture constituents effectively attract or repel each other. It can be shown that in this case v = ( 1 − X ) v 1 + X v 2 + ( 1 − X ) X v def , where v def is the defect of the volume additivity of the ideal-mixing case; here v def itself may be a function of X . Our model can in principle be formulated in the case of non-ideal mixing, but such a formulation would go well beyond the scope of our present work.

Since at fixed T , P = − d E / d V , the analog of the pressure mixing rule Equations 3, 4 for energy is E ρ , T = 1 − X E 1 ρ 1 , T + X E 2 ρ 2 , T . (7) In Equation 7 E , E 1 , and E 2 are the energies of the mixture and its constituents. Here, E 1 ( ρ , T ) and E 2 ( ρ , T ) are directly related to, respectively, P 1 ( ρ , T ) and P 2 ( ρ , T ) .

The “specific” Gibbs function—that is, the Gibbs function per unit mass which is consistent with the pressure mixing rule—is G = 1 − x G 1 + x G 2 . (8)

Indeed, the specific volume of the mixture ( ∼ 1 / ρ ) is V = ( ∂ G / ∂ P ) T ≡ G P = ( 1 − x ) V 1 + x V 2 , which is equivalent to Equation 4. We note that, although this formulation is not conventional, it can be found in the literature, such as in Duvall and Taylor (1971), where it was used for the description of the shock compression of a two-component mixture. It then follows from Equation 8 that the specific entropy of the mixture is S = − ( ∂ G / ∂ T ) P ≡ − G T = ( 1 − x ) S 1 + x S 2 . The isothermal compressibility ( 1 / B T ) is β = − ( 1 / V ) ( ∂ V / ∂ P ) T ≡ − G P P / V ; hence, β V = ( 1 − x ) β 1 V 1 + x β 2 V 2 . That is, 1 ρ B T = 1 − x ρ 1 B T , 1 + x ρ 2 B T , 2 . (9)

Note that Equation 9 follows directly from Equation 4 written as the total volume of the system being the sum of the volumes of its constituents: V = V 1 + V 2 . Then, under a small pressure change of Δ P , the total volume change is Δ V = Δ V 1 + Δ V 2 ; therefore, since B T ≡ − V ( ∂ P / ∂ V ) T , B T , 1 ≡ − V 1 ( ∂ P / ∂ V 1 ) T and B T , 2 ≡ − V 2 ( ∂ P / ∂ V 2 ) T , V B T = − Δ V Δ P = − Δ V 1 Δ P + Δ V 2 Δ P = V 1 B T , 1 + V 2 B T , 2 . (10) Equation 10 is equivalent to Equation 9 in view of Equation 4. We will also need the thermal expansion coefficient, α = ( 1 / V ) ( ∂ V / ∂ T ) P ≡ G P T / V or α V = ( 1 − x ) α 1 V 1 + x α 2 V 2 , which is equivalent to α ρ = 1 − x α 1 ρ 1 + x α 2 ρ 2 . (11)

Now, dividing Equation 11 by Equation 9, upon some algebra, we arrive at α B T = 1 − x α 1 B T , 1 + x α 2 B T , 2 + 1 − x x α 1 B T , 1 − α 2 B T , 2 ρ B T ρ 1 B T , 1 − ρ B T ρ 2 B T , 2 .

Note that Equation 9 allows the following parametrization: ρ 1 B T , 1 = ρ B T 1 ± C x , ρ 2 B T , 2 = ρ B T 1 ∓ C 1 − x , (12) where C is assumed to be a constant such that, regardless of the value of C , Equation 9 holds true because of the identity ( 1 − x ) ( 1 ± C x ) + x ( 1 ∓ C ( 1 − x ) ) = 1 . In the following, we keep the upper sign in front of C in the two denominators of Equation 12 so that C itself can be of either sign (or zero). We consider C as the only free parameter of this formulation, the value of which must be determined independently—based on the available experimental information on the liquidus of the system. We assume that no other experimental information is available; in particular, on the values of ρ B T , which would have otherwise allowed the calculation of the value of C directly from Equation 12. Then, via Equation 12, ρ B T / ( ρ 1 B T , 1 ) − ρ B T / ( ρ 2 B T , 2 ) = C , which we use in the above relation for α B T to finally obtain α B T = 1 − x α 1 B T , 1 + x α 2 B T , 2 + 1 − x x C α 1 B T , 1 − α 2 B T , 2 . (13)

We now consider a two-component mixture at temperature equilibrium at the liquidus point ( P , T ℓ ) and assume that each of the components is described by thermal EOS of the form Equation 2. Then, for the pressures of the mixture and its components at T = T ℓ , P ρ , T ℓ = P ρ , 0 + α B T T ℓ , P 1 ρ 1 , T ℓ = P 1 ρ 1 , 0 + α 1 B T , 1 T ℓ , P 2 ρ 2 , T ℓ = P 2 ρ 2 , 0 + α 2 B T , 2 T ℓ .

Since P 1 ( ρ 1 , 0 ) = P 2 ( ρ 2 , 0 ) = P ( ρ , 0 ) ≡ P (the pressure mixing rule), it then follows that ( P ℓ , 1 ≡ P 1 ( ρ 1 , T ℓ ) , P ℓ , 2 ≡ P 2 ( ρ 2 , T ℓ ) , and P ℓ ≡ P ( ρ , T ℓ ) ) α B T T ℓ = P ℓ − P , α 1 B T , 1 T ℓ = P ℓ , 1 − P , α 2 B T , 2 T ℓ = P ℓ , 2 − P , the use of which in Equation 13 multiplied by T ℓ , leading to P ℓ = 1 − x P ℓ , 1 + x P ℓ , 2 + 1 − x x C P ℓ , 1 − P ℓ , 2 . (14)

Equation 14 is our formula for the liquidus of a two-component mixture. It is in fact the new mixing rule mentioned above. Since both P ℓ , 1 = P ℓ , 1 ( T ) and P ℓ , 2 = P ℓ , 2 ( T ) are assumed to be available, provided that the value of C is known, this formula gives the value of the melting P of a mixture at any given melting temperature T ℓ * via P ℓ , 1 * = P ℓ , 1 ( T ℓ * ) and P ℓ , 2 * = P ℓ , 2 ( T ℓ * ) .

It is important to note the general features of the above formulas for the product of α B T and the liquidus of a mixture. Considering Equation 14 as a representative example, its general features are that it (i) is symmetrical under the simultaneous permutations 1 ↔ 2 , x ↔ ( 1 − x ) , C ↔ − C , (ii) satisfies the boundary conditions P ℓ ( x = 0 ) = P ℓ , 1 and P ℓ ( x = 1 ) = P ℓ , 2 , and (iii) satisfies the self-mixture condition where any pure substance can be considered a mixture with itself; hence, the choice of P ℓ , 1 = P ℓ , 2 = P should lead to P ℓ = P at any x . As clearly seen, this is achieved by the presence of ( P ℓ , 1 − P ℓ , 2 ) on the right-hand side of Equation 14.

Both T ℓ ( 0 ) in Equation 15 and T ℓ ′ ( 0 ) in Equation 16 as functions of x have a local extremum (either minimum or maximum) at X = C − 1 2 C . (17)

If this extremum does occur for the mixture, then 0 < X < 1 , which puts a constraint on C : | C | > 1 , such that C > 1 ⇔ 0 < X < 0.5 , C < − 1 ⇔ 0.5 < X < 1 . (18) Equation 18 is consistent with the invariance under the simultaneous permutations x ↔ ( 1 − x ) and C ↔ − C . In this case, the values of T m ( 0 ) and T m ′ ( 0 ) at x = X are T ℓ X 0 = C − 1 2 T ℓ , 1 ′ 0 T ℓ , 2 0 − C + 1 2 T ℓ , 2 ′ 0 T ℓ , 1 0 C − 1 2 T ℓ , 1 ′ 0 − C + 1 2 T ℓ , 2 ′ 0 (19) and T ℓ ′ X 0 = 4 C T ℓ , 1 ′ 0 T ℓ , 2 ′ 0 C − 1 2 T ℓ , 1 ′ 0 − C + 1 2 T ℓ , 2 ′ 0 . (20)

If, however, | C | ≤ 1 , the extremum occurs either outside the physical region 0 ≤ x ≤ 1 of the mixture or at one of its endpoints (if | C | = 1 ), both the local maximum and minimum are at the two endpoints, and T ℓ ( 0 ) is a monotonically decreasing or increasing function of x . In this case, Equations 19, 20 do not apply.

Based on general considerations, it is expected that C = O ( 1 ) . In the following example of the application of our model to a real system, C is identically 1. The system that we are considering here is a mixture of silicon and germanium which, according to Olesinski and Abbaschian (1984), form a continuous solid solution without any eutectics. According to the literature (Jayaraman et al., 1963; Deb et al., 2001; Yang et al., 2004; Yang and Jiang, 2004; Kubo et al., 2008; Pasternak et al., 2008), the low- P melting curves of Si and Ge are, respectively, T m Si ( P ) = 1690 − 63 P and T m Ge ( P ) = 1210 − 50 P . Hence, this is an example of a mixture of “anomalous” melters, which both exhibit T m decreasing with P (for the vast majority of substances, T m increases with P , which is the normal case). Hence, in Equation 15 T ℓ , 1 ( 0 ) = 1690 , T ℓ , 1 ′ ( 0 ) = − 63 , T ℓ , 2 ( 0 ) = 1210 , and T ℓ , 2 ′ ( 0 ) = − 50 . Then, the best fit of the form Equation 15 with the above parameters to the experimental data of Olesinski and Abbaschian (1984) brings up C = 1 . In Figure 2, the resulting T ℓ ( x ) is compared to the liquidus of Si 1 − x Ge x from the experiment.

The computational details of our QMD simulations can be found in our previous work on this subject (Burakovsky et al., 2023; Burakovsky et al., 2024b). The mixed U O 2 - PuO 2 system is modeled as a substitution alloy in which some U atoms (at randomly chosen lattice sites) are replaced by Pu atoms according to the corresponding PuO 2 content. Specifically, for our simulations of a 768-atom U 0.25 Pu 0.75 O 2 system, 576 U atoms chosen at random are replaced with 576 Pu atoms (alternatively, 192 Pu atoms of pure 768-atom PuO 2 systems are replaced with 192 U atoms chosen at random).

Our melting simulations are carried out using the so-called Z method (Burakovsky et al., 2015). In these simulations, the U 0.25 Pu 0.75 O 2 supercell is subject to a set of initial temperatures ( T 0 ) separated by an increment of 250 K and run with QMD in the N V E ensemble, for a total of 5,000–6,000 time steps of 2.5 fs each—up to a total of 15 ps of running time—to determine T m at the corresponding melting pressure ( P m ) . 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 (shown as solid lines in Figures 4 and 5). In this case, P m = 0 corresponds to the unit cell of a = b = c = 5.58 Å or ρ = 10.50 g/ c m 3 .

As seen in Figures 4 and 5, since the beginning of the run, after ∼ 600 time steps (1.5 ps), T decreases and P increases (since in the N V E ensemble the total energy E ∼ k B T + P V is conserved). This is a signature of a superionic transition. Due to a 15-fold difference in the atomic masses of U/Pu and O, the O sublattice becomes less stable than for U/Pu, and when sufficiently high T it disorders first, such that the anions ( O − ) start flowing through the ordered structure of the cations ( U + / P u + ). Such a (superionic) phase transition accompanied by a rapid increase in ionic conductivity has been observed in many diatomic systems. It was observed in both pure UO 2 (Dworkin and Bredig, 1968) and PuO 2 (Chroneos et al., 2015; Günay et al., 2016). Figures 4 and 5 demonstrate its occurrence in U 0.25 Pu 0.75 O 2 .

Thus, the first drop in T (increase in P ) corresponds to the activation of the O flow. This process takes ∼ 1000 time steps ( ∼ 2.5 ps). The system of quasi-static cations and mobile anions then equilibrates, and the second drop in T (increase in P ) occurs after a total of ∼ 1400 time steps ( ∼ 3.5 ps). This second drop in T is associated with the disordering of the U/Pu sublattice—a true melting transition—and the corresponding ( P , T ) point lies on the system’s liquidus. The melting process takes ∼ 600 time steps (1.5 ps). The emerging liquid equilibrates at ( P , T ) ≈ ( 0 , 2900 K ) , which is the liquidus point of U 0.25 Pu 0.75 O 2 .

Hence, our ab initio ambient melting point of U 0.25 Pu 0.75 O 2 appears to be 2900 ± 135 K. Uncertainty of the value of T m intrinsic to the Z method is 125 K, half of the increment of T 0 (Burakovsky et al., 2015), which turns out to constitute ∼ 4 % 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 / d P at P = 0) is ∼ 90 K/GPa, as predicted by the model (Equation 16), a P m uncertainty of 0.5 GPa translates into a T m uncertainty of δ T m ∼ d T m d P δ P ∼ < 45 K. Therefore, the combined uncertainty of T m in our QMD simulations is ≈ 135 K, which is within 5% of T m . Thus, our results on the ambient T m of U 0.25 Pu 0.75 O 2 are expected to be quite accurate overall.

Our value of T m ( 0 ) = 2900 K for U 0.25 Pu 0.75 O 2 is ∼ 200 K above the value of ∼ 2700 K suggested by our previous QMD simulations of smaller systems (Burakovsky et al., 2023). Hence in our simulations, some size effects are indeed present and should therefore be taken into account. Assuming that the results on smaller supercells at x = 0.7 and 0.9 should be corrected by adding ∼ 200 K to the corresponding values of T m s, the two corrected values would be fully consistent with the theoretical liquidus curve.