## 2 Methods
### Melting Simulations  
Both UO\({}_{2}\) and PuO\({}_{2}\) maintain their ambient cubic fluorite structure (Fm3m) 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\) A [23] while for PuO\({}_{2}\)\(\rho\ =\ 11.7\) g/cm\({}^{3}\), \(a=b=c=5.39\) A [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  
Figure 1: The ambient fluorite structure of both UO\({}_{2}\) and PuO\({}_{2}\). Uranium/plutonium atoms are shown in green, and oxygen atoms in orange.  
\(10^{-6}\) eV/atom and that of Hellmann-Feynman forces on each atom to 0.01 eV/A. 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 \((T_{0})\) 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\
eq 1\), i.e., up to 25 ps in either case, 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 which are shown as solid lines in Figures 2-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. Figures 2-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 actinides subsystems during the QMD run results in larger-than-usual amplitude of the oscillations, up to \(\sim\)500 K for \(T\) and \(\sim\)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 \(\sim\)2% of \(T_{m}\). Uncertainty of the value of \(P_{m}\) in our simulations is \(\sim\)0.5 GPa. Assuming that the initial slope of the melting curve \((dT_{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{dT_{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.