Thermodynamic non-ideality and disorder heterogeneity in actinide silicate solid solutions

Non-ideal thermodynamics of solid solutions can greatly impact materials degradation behavior. We have investigated an actinide silicate solid solution system (USiO4–ThSiO4), demonstrating that thermodynamic non-ideality follows a distinctive, atomic-scale disordering process, which is usually considered as a random distribution. Neutron total scattering implemented by pair distribution function analysis confirmed a random distribution model for U and Th in first three coordination shells; however, a machine-learning algorithm suggested heterogeneous U and Th clusters at nanoscale (~2 nm). The local disorder and nanosized heterogeneous is an example of the non-ideality of mixing that has an electronic origin. Partial covalency from the U/Th 5f–O 2p hybridization promotes electron transfer during mixing and leads to local polyhedral distortions. The electronic origin accounts for the strong non-ideality in thermodynamic parameters that extends the stability field of the actinide silicates in nature and under typical nuclear waste repository conditions.


INTRODUCTION
Uranothorite (U x Th 1−x SiO 4 ) is the only known natural silicate system where actinides, U and Th, form a continuous solid solution between thorite (ThSiO 4 ) and coffinite (USiO 4 ) 1-12 . Uranothorite is of geological interest because of its natural occurrence in uranium deposits 4,5,7,[13][14][15][16] , suggesting an enhanced stability field for the coffinite endmember. Furthermore, because zircon (ZrSiO 4 ) has been proposed as a possible phase for the incorporation and immobilization of Pu [17][18][19][20][21][22] , uranothorite, as an actinide solid solution system based on the zircon structure (I4 1 / amd), can provide critical insights into coexisting actinides in this crystalline matrix. How the mixing of various actinides impacts the structural and thermodynamic stabilities is essential to the understanding of behaviors of Pu in the zircon structure type. On the other hand, coffinite is known to be one of the key alternation phases of uraninite (UO 2 in nature) or UO 2 in the spent nuclear fuel (SNF) in a geologic repository [23][24][25][26][27][28][29][30] due to the contact with groundwater under anoxic and silica-enriched conditions 31,32 . However, whether uranothorite or other actinide silicate solid solutions can be more important sinks for U, Pu, and other minor actinides that control the dissolution and distribution of radionuclides is unknown and requires a fundamental thermodynamic assessment. Furthermore, the degradation or formation processes may also be influenced by material morphology. Particularly, the presence of nanophases and surface defects may be catalytic and initialize intragranular cracking or intragranular corrosion, expediting the degradation process 33 . Therefore, to demonstrate the feasibility of using zircon as a ceramic waste form for the immobilization of multiple actinides (e.g., U and Pu) under geological repository conditions and predict the fate of this actinide silicate solid solution under relevant mineralization and weathering processes, a better understanding of the thermodynamics and structures of uranothorite are essential.
These important differences between the two isostructural endmembers raise the key question of whether their intermediate compositions (uranothorite) carry transitional structural and energetic characteristics and whether those differences control the thermodynamically and kinetically driven degradation or transformation processes. Due to the similar chemical characters of actinide cations and close ionic radii of 8-coordinated Th 4+ , U 4+ , and Pu 4+ (1.19, 1.14, and 0.96 Å, respectively) 51 , an ideal mixing behavior in the zircon structure was initially expected. However, based on recent experimental work, as well as the present study, U x Th 1−x SiO 4 compounds have both structures and thermodynamic parameters that deviate negatively from the linear Vegard's law, suggesting non-ideality behaviors 12 . Specifically, a regular solution model was proposed for uranothorite given the negative enthalpy of mixing, which favors the formation of uranothorites by lowering their formation energies that facilitates accommodating U in the zircon structure 12 . By contrast, Ferriss et al. proposed different solid solution models using the density functional theory (DFT) and Monte Carlo methods by taking into account the four nearest homo-or hetero-adjacent atom interaction energies, which suggests positive mixing enthalpies and phase exsolution 52 . The disagreement between experimental and computational results motivated this comprehensive structural study on uranothorite, coffinite, and thorite to investigate (i) the origins of the thermodynamic non-ideality by investigating the atomic-scale arrangements and the electronic features in the local atomic-scale, (ii) the entropic term which has contributed significantly to the thermodynamic stability at elevated temperature (80-300°C generated from the decay heat) 53 , as is described by the mixing model of U and Th, and (iii) geochemical implication if a multi-element solid solution system, such as uranothorite, can extend the stability field of the immobilized actinides (e.g., U and Pu) in the waste matrix when stored in a geologic repository.
With the structural insights, we further revisited the previously published calorimetric data 12 and refined them with the updated subregular solid solution model associated with a negative enthalpy of mixing. The subregular model describes the excess Gibbs free energy or the enthalpy of mixing by Ωx(1−x), where Ω is concentration-dependent and controlled by two interaction parameters, Ω = w 12 (1−x) + w 21 x. Though high temperature calorimetry has been demonstrated to be sensitive to bond strength, it does not reveal the origin of the enthalpic stabilization for intermediate compositions. Generally speaking, such a negative deviation from ideal mixing can be attributed to (i) short-range ordering, (ii) local distortions of coordination polyhedra, (iii) electronic features, or any combination of these three possibilities 54 . All three aspects were examined in this work by neutron total scattering data of thorite, coffinite, and three uranothorite samples (U 0.11 Th 0.89 SiO 4 , U 0.35 Th 0.65 SiO 4 , and U 0.84 Th 0.16 SiO 4 ).
The diffraction data were analyzed by the Rietveld method for long-range structural information, and the diffuse scattering was subject to pair distribution function (PDF) analysis using real-space fitting that yields short-range structural features, such as the extent of cation ordering (or clustering) and polyhedral distortion. PDF analysis provides a direct determination of the atomic correlations. Fitting a PDF of a crystalline material may be performed using real-space Rietveld refinement or reverse Monte Carlo (RMC) modeling methods, both of which remove symmetry operations (unless introduced as constraints to the fit). In this work, both fitting techniques were used to explore the atomicscale arrangement (ordering, clustering, or disordering) of U and Th in the zircon structure, which was further assessed by unsupervised machine learning (ML) (K-means clustering) to reveal the mid-range (~2 nm) arrangement of cations in their sublattices as belonging to various clusters. Additionally, DFT calculations were performed on coffinite and thorite to reveal their electronic structures, particularly actinide covalency within the metal polyhedra, the difference of which may cause the thermodynamic non-ideality when U and Th are mixed.
Lastly, we performed multi-component system thermodynamic calculations using the refined subregular solid solution model from this work and the reported or extrapolated thermodynamic properties (enthalpy, Gibbs free energy, heat capacity and entropy) to evaluate the phase equilibria of uranothorite under various oxygen fugacities and moderately elevated temperatures [55][56][57] . The outcome of the thermodynamic calculations allows the prediction of whether uranothorite will be stabilized or be subject to degradation during possible mineralization and weathering processes in nature or in a geological repository.

RESULTS AND DISCUSSION Neutron diffraction and scattering data analyses
The long-range crystalline structures of coffinite, uranothorite, and thorite were refined by performing Rietveld analysis of the neutron diffraction (ND) data. Compared to X-ray diffraction, ND is more sensitive to oxygen and more suited to distinguish neighboring elements (e.g., U and Th) which have similar electron densities. The experimental diffraction data and fitted patterns can be found in Supplementary Fig. 1. From the final convergence of the refinements, the fitted unit-cell parameters, volumes and R wp values (ranging from 3.4% to 7.6%) are shown in Table 1. The obtained atomic coordinates of oxygen and atomic displacement parameters are listed in Supplementary Table 1. The real-space Rietveld refinement of the neutron total scattering data was then performed using PDFgui with the experimental G(r) data and fits shown in Supplementary Fig. 2, and fitted values listed in Table 1. Because there is no superstructure in the model defined in PDFgui, the fitted result represents a random atom configuration, in which U and Th atoms are considered to be fully disordered in the crystallographic sites. The fitted unit-cell parameters, atomic coordinates, displacement parameters, and R wp values (ranging from 17.4% to 28.9%) are summarized in Table 1 and Supplementary Table 2. These two obtained sets of lattice parameters are general in agreement with previous values determined from X-ray diffraction 12,50 . The unit-cell volume exhibits a slight increase with increasing Th addition (Fig. 2). This trend is consistent with the previous findings 12 and is supported by the substitution of U 4+ (1.00 Å, Shannon ionic radius in eightfold coordinated) by Th 4+ (1.05 Å) 51 .
To explore the mixing geometries of U and Th atoms beyond short-range distances (i.e., first three nearest-neighbor metal cation) and the extent of their ordering or disordering in the matrix, RMC modeling was performed on the scattering data of the three uranothorite samples (U11, U35, and U84). Compared to the real-space Rietveld method (PDFgui) that fits the PDF from a starting crystallographic model assuming infinitely large crystallographic planes, the RMC fitting is performed on the normalized reciprocal space scattering pattern, F(Q), and the PDF by varying the positions of atoms in a predefined atom "box" following a set of constraints, and then periodically comparing the model to data, which allows for the investigation of the arrangements of U and Th atoms in a 10 × 10 × 10 supercell (see the "Methods" section and Supplementary Information for details). Figure 3 provides the fitting results of RMC modeling. The starting mixing configurations for the three uranothorites are locally ordered structures (Fig. 4). The divergence of U and Th from the initial ordered models was achieved by atom swapping and translation per refinement step, which was constrained by distance window (Supplementary Table  3) and bond valence sum (BVS) constraints (Supplementary Table  4). Here we examine the atomic positions relaxed after RMC modeling as the "true" configuration to evaluate the mixing behavior of Th and U, in comparison with their "ordered" ones in the supercells (Fig. 4). Previously, the variation of the a lattice parameter through X-ray diffraction suggested that chemical ordering may occur in a direction parallel to the a-axis 12 . However, in our current assessment from the neutron scattering data, there seems to be no strong indication of ordering occurred in either the a-or the c-axis, which was also reflected in Fig. 4. There is also no preference of U-Th clustering along the (010) or (001) planes. The representations of the final atom models perpendicular to the Table 1. Unit-cell parameters obtained from ND Rietveld refinement and PDFgui refinement.

Sample
Composition  c-axis may be slightly misleading to visual inspection with respect to the extent of clustering in the systems analyzed. This is because U and Th atoms are closer-packed along the (001) plane; however, the (010) plane is actually closer-packed with 2 U/Th atoms per plane as compared to 1 U/Th atom per plane for the (001) plane. In order to understand how the arrangements of U and Th atoms propagate to nanoscale, a further detailed categorization of the distribution of all U and Th atoms in the 10 × 10 × 10 supercell was quantified by an unsupervised ML algorithm using the atomic arrangement from the RMC-relaxed superlattice "box" before and after refinement (Fig. 5).
Further, the U-O, Th-O, and U-Th partial PDFs, g(r), for the three uranothorites obtained from RMC modeling and PDFgui fitting are shown in Supplementary Fig. 3. The integration of g(r) over distance r yields the coordination number distribution function N (r) which is shown for each pair, U-O, Th-O, and U-Th in Supplementary Figs. 4 and 6. Three coordination shells are highlighted in red, yellow, and blue for the U-Th pairs in their nearest neighbors in the a-b planes, the second nearest, and the third in the a-b planes, respectively.

Solid solution model of uranothorite
By using the refined ND data, the unit-cell volume (V cell ) in terms of USiO 4 mole fraction x is fitted by the following polynomial function: The excess volume (ΔV mix ) can then be modeled by a subregular volumetric model based on Eq. (1), where ω G1 = 11.6 Å 3 and ω G2 = −41.7 Å 3 are the asymmetric volumetric interaction parameters. Such a relation follows a subregular model, which is similar to the previous finding that ΔV mix has a negative deviation from Vegard's law (ideal mixing) behavior. This updated model suggests more complicated mixing behavior in uranothorite which is asymmetric with the minimum volume of mixing located around 70 mol% of U (Fig. 2b). This updated model also explains our recent work on the hightemperature behavior of uranothorite, where U 0.90 Th 0.10 SiO 4 was found to decompose above 800°C to phases firstly with a decreasing trend of unit-cell volume, presumably having U content move from x = 0.9-0.8, and then a reversed increasing trend of unit-cell volume at higher temperatures, corresponding to the formation of U-depleted uranothorite (x < 0.8) 50 . The continuing shifts in chemical composition are driven by thermodynamics due to the mixing-induced enthalpic stabilization, evidenced by a previous calorimetric study where negative enthalpies of mixing (ΔH mix ) were reported 12 . Inspired by the subregular volumetric model better describing the uranothorite solid solution, we revisited the published ΔH mix data 12 , and refined them with a subregular model, where w 12 = −88.6 kJ/mol and w 21 = −155.2 kJ/mol, with adjusted R 2 = 0.9888 (Fig. 2c). The negative interaction parameters w 12 and w 21 indicate that the uranothorite solid solution is expected to be stabilized by the enthalpic effect. Such negative deviations from Vegard's law in both formation enthalpy and unit-cell volume data could originate from (i) a favorable mixing between Th and U in the local range, (ii) the distortion of U or Th occupied polyhedra, and/or (iii) the difference in bonding characters (5f 0 configuration for Th 4+ and 5f 2 configuration for U 4+ ) 58,59 . First, based on the results from PDF analysis, we compared the coordination number distribution  function N(r) of U-O, Th-O, and U-Th pairs determined from the RMC-derived models with the same quantities from a hypothetical random distribution model and show the comparison across three coordination shells in Fig. 6. No noticeable difference is observed, or there is no particular feature in the data that would lead to the locally ordered model. This is consistent with the PDFgui refinements which use the disordered model that produces good fits with acceptable R wp values.
In order to further quantitatively evaluate the degree of cationic ordering, we calculated the Warren-Cowley (WC) parameters 60,61 of U and Th based on the atomistic models from RMC fitting, and compared with the ordered and random configurations. Here the WC parameter, α i,j (r), is defined as follows: where P i;j r ð Þ is the probability of finding an atom of type j at r given an atom of type i at the origin, and c j is the average concentration of j-type atoms. Thus, the U-centered WC parameter α U,Th has a value ranging from −N Th /N U to 1 (and a similar definition for Th-centered WC parameter α Th,U ). For an assumed ideal random distribution of U and Th atoms in the sublattice sites, α U,Th (or α Th,U ) is zero and depicted as the dashed lines in Fig. 6. A positive value indicates a tendency for like pairing (i.e., clustering or formation of domains) and the positiveness with the upper limit of 1 represents the "strength" of such a tendency. A negative value indicates the tendency to form unlike pairs (U-Th pairs). Thus, on each U/Th ratio, we compared their WC parameters calculated from RMC-determined atomistic models and those from ordered and random configurations. In Fig. 6, the first three polyhedral coordination shells were color coded, with the coordinates in red, yellow, and blue corresponding to the U-Th pairs in their nearest neighbors in the a-b planes, the second nearest, and the third in the a-b planes, respectively, which are defined by the coordination number distribution function N Th-U (r). The calculated α, superimposed as open circles, are in almost-zero or very weakly negative values in the region for all three solid-solution compositions within the first three coordination shells, again suggesting the disordering of U and Th atoms. Despite the fact that a random mixing model is more probable, the ordered model cannot be completely rejected. Potential ordering may have also been underestimated due to the randomization of atoms which is intrinsically built into the RMC algorithm 62 . A low level of shortrange ordering may still be possible, as suggested by the slight negativity of WC parameters within a certain coordination shell (e.g., the first one in U11 or the third one in U35), because P i;j r ð Þ is obtained by statistically averaging all same-type atoms in the supercell and thus less sensitive to some finer local information.
Thus, ND in this work, previous XRD analysis 12 , and EXAFS analysis (performed on similar uranothorite compositions in a study by Labs et al.) 11 are all consistent with the above RMCderived interpretation and WC parameters. Particularly, Labs et al. 11 found no evidence for the formation of secondary phases and no significant degrees of atom-clustering through EXAFS and powder diffraction of hydrothermally-grown uranothorites. Lastly, a previous DFT calculation on uranothorite predicted a tendency to form local clusters with a positive enthalpy of mixing 52 , which contradicts the present experimental evidence. This is not too surprising when the first contribution to non-ideality is mostly considered that is usually important only when mixing cations have large size differences, thus less probable in uranothorite due to the close cation sizes of U 4+ (1.00 Å) and Th 4+ (1.05 Å).
To further examine the atomistic models and how the local disordering propagates to the outer coordination shells and further to macroscopic-scale structures, the unsupervised ML algorithm (K-means clustering, details provided in the "Methods" section) was implemented on the large dataset of coordinates of U and Th atoms in uranothorite supercells. Figure 5 provides the inertia-cluster number (used to determine the appropriate number of clusters per cell) and cluster sizes of U or Th atoms in uranothorite supercells (end member compositions yielded non-physical values and have been omitted). Noted that here the "cluster" does not represent non-repeating local U-U or Th-Th clusters in the conventional structural sense that is opposite of short-range ordering. Instead, in K-means clustering algorithm, it is a group of U or Th atoms that are found by ML to be of similar types. Thus, the distribution of such "clusters" can shed light on the local patterns of U or Th from the impact of cation mixing. Gaussian functions were used to fit the distributions (Supplementary Table 5). The K-means clustering results on RMC fits are compared with those of the "ordered configuration" (used as the "initial" model, exhibited in Fig. 4, for RMC modeling). The RMC modeling leads to a slightly broader distribution than the "ordered" configurations (Supplementary Table 5). For instance, the ordered model of U11 has σ = 1.9 compared to the RMC model which has σ = 3.3, where σ is the variance in the number of atoms per cluster. Therefore, the cluster sizes (i.e., the number of atoms per cluster) of RMC models overall resemble those of the ordered models in the med-range scale.
This is an interesting finding, suggesting that although U and Th are mixed randomly in the cationic sublattice, such disordering could locally differ, forming nanosized domains similar to what the ordered arrangement may introduce, despite the drastic change in the atom positions (Fig. 5). Specifically, the statistical results of Kmeans clustering on RMC-derived models determined the size of U "clusters" in U11 or Th "clusters" in U84 to be approximately 5-20 atoms (mean value 9.8 atoms), suggesting that U or Th cations form similar local environments in the size of the second cation coordination shells (color coded in yellow in Fig. 6). Such a mixing-induced feature is more pronounced near the equimolar composition. U35 has approximately 20-40 atoms for U "clusters" (mean value 27 atoms), generally equivalent to the size of the third cation coordination shell. Such locally contained disordering could cause the formation of nanodomains and hence nanosized heterogeneity. Similar disorder heterogeneity was recently highlighted in pyrochlore solid solution systems 63 . Nanoscale disorder heterogeneities may change the interpretation of the thermodynamic properties of uranothorite, possibly affecting chemical processes such as the selective dissolution of U and Th.
Structural and electronic origins for the thermodynamic nonideality With a cationic disordering model, instead of a short-range ordering model, suggested by RMC for uranothorites, the local distortion in U-and Th-occupied polyhedra becomes a more probable origin for the thermodynamic non-ideality. Note that the MO 8 polyhedra (M = U, Th) are triangular dodecahedra that can be considered as two intersecting tetrahedra with two sets of O 4 atoms (Fig. 1b), M(O1) 4 and M(O2) 4 6,64 . The first set consist of four oxygen atoms (O1) forming longer bonds with U (2.421 Å) 11 or Th (2.467 Å) 2 , described as the extended tetrahedra, sharing the edges with SiO 4 tetrahedra aligning along the c-axis. Such <U-O1> or <Th-O1> bonds are not very sensitive to U and Th composition, which is the reason for having a near linear trend of lattice parameter c when varying the U content 12,65 . On the other hand, four other oxygen atoms (O2) have shorter bonds with U (2.292 Å) 11 or Th (2.384 Å) 65 , described as the compressed tetrahedra, sharing corners with other compressed tetrahedra along the a-b planes, which are more sensitive to substitution 12,65 . Corner-sharing also make the disphenoidal M(O2) 4 easier to be relaxed to distorted geometry under chemical strain (e.g., cationic substitution) or physical strain (e.g., compression). The different susceptibility of the two types of tetrahedra to substitution is confirmed by RMC's results that <U-O2> expands by 2.4% while <U-O1> only expands by 1.4% when U content decreases to 11 mol% in uranothorite. Similarly, in Th sublattices, <Th-O2> shrinks by 3.2% compared to 2.2% for <Th-O1> (Table 2). Such local features have impacts on the long range which leads to a greater magnitude in unit-cell variance with an anisotropic behavior. Mixing introduces greater changes along the a (and b)-axis due to the flexibility of the second-type oxygen (O2), while encountering more resistance along the c-axis because of the insensitivity of the first-type oxygen (O1) and the more rigid edgesharing M(O1) 4 -Si(O1) 4 chains.
Furthermore, coffinite and thorite have different bond-length ratios (<M-O1>/<M-O2>) of 1.057 and 1.035, respectively, which suggests that the corresponding dodecahedra have slightly different local environments. This could be due to the distinct covalent tendency of U coupled with crystal field splitting similar to the Jahn-Teller distortion 54 , so that U 5f can strongly overlap with O 2p orbitals to form hybridization 66,67 . In order to investigate the bonding structure and its impacts on local geometry, we performed DFT + U calculations on coffinite and thorite. The U covalency is more directional along a-b planes in the compressed tetrahedra resulted in a shorter bond length of <U-O2>. The structure optimized by DFT + U method is in an antiferromagnetic (AFM) ordering, in agreement with the experimental results and previous calculations 68 (Supplementary Tables 6 and 7). As shown in the projected density of states (pDOS, Fig. 7a, d), the overlap between the U 5f orbital and O 2p orbital density of states from −8 to 0 eV demonstrates a strong hybridization of these two orbitals, which provides a covalency signature for the <U-O> bonding. On the contrary, a limited overlap of Th 6d orbital and O 2p orbital density of states from −8 to 0 eV proved a marginal hybridization and poor covalency of <Th-O> bonding. Furthermore, more hybridization between Th and O emerges above the Fermi level, which can enhance anti-bonding and decrease the <Th-O> covalent bonding strength.
To further understand the local bonding structure, charge density and electron localization function (ELF) were calculated and depicted in Fig. 7 and Supplementary Fig. 4 69 . Figure 7c, f shows that although most charge density is localized near U or Th nuclei, indicating a strong ionic bond signature, a small proportion of electron density is distributed in the interstices of <U-O> in coffinite, which is consistent with a previous DFT result that U 5f was found to contribute to <U-O> bonds and exhibit a mixture of covalent/ionic characteristics 68 . More importantly, contour plotting reveals that electron density of U overlaps more with O2 than O1, hence the <U-O2> bonding has more covalent character. The excess covalency in <U-O2> was observed to be more than that in <Th-O2> due to the excess density overlapping in <U-O2> bond (Fig. 7c, f). Specifically, in the USiO 4 matrix, oxygen attractors are in the shape of a meniscus heading to nearby U attractors, which indicates strong polarization between U and O attractors ( Supplementary Fig. 5). While in ThSiO 4 , this polarization level is smaller and hence the local electron density is nearly spherical ( Supplementary Fig. 5). Lastly, ThSiO 4 was found to have higher ELF values, again indicating stronger ionic bonding of <Th-O> bonds. In general, ionic bonding is stronger than covalent bonding, and the difference in bonding nature may also explain why ThSiO 4 is more thermodynamically stable than USiO 4 .
Such a difference in <U/Th-O> bonding also makes the M(O2) 4 tetrahedra more "compressed" in coffinite than in thorite ( Table 2 and Fig. 8). Interestingly, when U and Th mix in uranothorite with comparable amounts, such local polyhedral differences become "smeared" (Fig. 8), with UO 8 dodecahedra having more ThO 8 -like features and vice versa. Additionally, another shorter metal-oxygen bonds <M-O3> emerge upon mixing, further distorting the local disphenoidal geometries and forming an irregular triangular dodecahedron ( Table 2). The UO 8 polyhedral distortion could be a result of an initial distortion introduced by substituted Th that breaks the local symmetry and then attenuates the crystal field splitting, leading to the weakening of the covalency of <U-O2> and its elongation. Thus, the distortion in UO 8 is more pronounced in the U(O2) 4 tetrahedra, which makes it resembles its counterpart in ThO 8 . While in the U-or Th-rich uranothorite, the introduction of different type cation starts to distort the two disphenoidal M(O) 4 differently. We proposed an electronic interaction mechanism between UO 8 and ThO 8 dodecahedra to explain the variance in <M-O1>/<M-O2> bond ratio. Considering the excess covalent bonding property in <U-O2> bond could result in additional electron density accumulation in O2 atom, when U starts to enter the Th-rich matrix or vice versa, more electron density can transfer from <U-O2> bond to <Th-O1> bond through the bridging O2 atom (the same oxygen labeled as O1 atom for Th, Fig. 8b). Electron depletion from <U-O2> results in a less stable and hence longer <U-O2> bond while leaving <U-O1> untouched since it bridges to Th as O2 atom and remains mostly ionic. This results in a reduction of <U-O1>/<U-O2> bond ratio (Fig. 8a). On the other hand, electrons density flows into the <Th-O1> bond and could fill in the anti-bonding orbitals according to the pDOS and destabilize the <Th-O1> bond with an extended bond length. This also explains the increase of <Th-O1>/<Th-O2>. Such electron transfer interaction may be the main origin of the nonideal thermodynamic behavior of uranothorite.
Thermodynamic modeling of uranothorite and implications for material degradation As a conclusion of RMC analyses and DFT calculations, the above proposed geometrical distortion coupled with the alteration in electronic structures is responsible for stabilizing intermediate U/ Th composition in uranothorites with negative enthalpies of mixing and volume of mixing 54,70 . Another important implication for thermodynamics of uranothorite is that the configurational entropy may be determined by using the Boltzmann entropy formula (ΔS Boltzmann ¼ À4R 71,72 , where R is the ideal gas constant and x is the U concentration in the dodecahedral site. The configurational entropy due to mixing for uranium concentration of 0.11, 0.35, and 0.84 were calculated to be 2.88, 5.38, and 3.65 J·mol −1 ·K −1 (Supplementary Table 8). Such determination is usually not an easy case for other non-ideal solid- solution systems, because only random distribution of mixed ions will make the configurational entropy possible to be evaluated by the Boltzmann entropy ΔS Boltzmann , and certain favorable atomic arrangements in the short range will lead to the reduction from the theoretical maximum value 70,[73][74][75][76] . For a full thermodynamic assessment, vibration entropy and heat capacity (C p ) are also needed. However, due to the lack of experimentally determined C p and entropy data, we estimated them for USiO 4 Table 9). The vibrational entropy was also estimated from the integration of C p . Although there is even no experimental data for endmembers, USiO 4 and ThSiO 4 , we applied different theoretical methods to obtain heat capacities, including density functional perturbation theory-based calculations by Phonopy 78 , Neumann-Kopp's rule 79 , and n-dimensional Debye model 80 . A comparison of C p derived from the three methods and from previous work is shown in Supplementary Figs. 6-8, with their integrated entropy values summarized in Supplementary Table  9. All values are in general consistent within the error range, except for the n-dimensional Debye model, due to its limited application on multi-atom ceramic materials. By accepting the volume-entropy empirical method, the vibrational mixing entropies of uranothorite were further determined to be 0.20, −0.26, and −0.87 J·mol −1 ·K −1 , respectively for U11, U35, and U84 samples (Supplementary Table 8). Thus, configurational entropy dominates the entropy contribution based on these two   Table 8) can contribute positively to the thermodynamic stability of uranothorite at 300°C (peak temperature for the nuclear repository backfill barrier materials) 53 by −6 to −12 kJ/mol, and at 1000°C (ceramic waste processing temperature) 22,47,71 by −14 to −27 kJ/mol, respectively.
The availability of C p , ΔH, ΔS, and ΔG of this solid solution system from this work and from literature 12,40,49,81 now enables evaluation of changes in mineral's stability fields triggered by the presence of uranothorite in complex systems, which in this work were simplified by including the following solid phases with their thermodynamic parameters: coesite 82 2 85 , coffinite 40,49,81 , and thorite 12,49 . This evaluation was performed through thermodynamic modeling using the Hch/Unitherm code 87,88 . The data for uranothorite, coffinite, and thorite used in these calculations are summarized in Table 3. Calculations of equilibrium compositions were performed on three model systems by analyzing limited numbers of factors (temperature and redox condition) that can affect the formation of (i) coffinite in the first, (ii) uranothorite described by the ideal mixing model in the second, and (iii) uranothorite described by the revised subregular solid solution model (from section "Solid solution model of uranothorite") in the third calculation. Here the redox parameter was set in the calculations through fixed fugacity of O 2 , gas, the thermodynamic properties of which were taken from refs. 83,89 . The range of log(fO 2 ) was chosen from −40 to −8 to simulate a wide range of natural redox environments. Temperature effects were evaluated from 25 to 60°C, which, although lower than the average nuclear waste repository temperature, is relevant to the formation, weathering, or material degradation processes in nature. An additional reason for selecting such a narrow temperature interval is also in limited applicability of the mixing model developed for uranothorite: this model was developed for ambient conditions and its applicability to higher temperatures is not evident. However, we believe that the effects observed even for such a narrow temperature interval can serve as a good illustration of how accounting for the presence of solid solutions can change the predicted phase stability fields. These effects can potentially be manifested to much higher degree at higher temperatures.
These calculations can reveal the impacts of the uranothorite solid solution system and mixing models describing the system on the stability of U-bearing zircon-type phases under various log (fO 2 )-temperature conditions. The starting composition of the model system was selected to be 1 mol of thorite and 1 mole of coffinite for the solid solution system for simplification. Significant differences were identified for the three models ( Fig. 9 and Supplementary Table 9). Without the presence of uranothorite, coffinite only starts to be stabilized at 45°C but under very reducing conditions, log(fO 2 ) < −36 (Fig. 9a). This result is consistent with the natural observation that coffinite precipitation was found no greater than 50°C in the uranium diagenetic unconformity-related deposits in the Athabasca basin, where Th is present in only trace concentrations and uranothorite is not formed 90,91 . When uranothorite is possible to form (model 2), the formation temperature is lowered down to 40°C under very reducing conditions (Fig. 9b). At 55°C and above, due to the very high stability of thorite, the formation of uranothorite also becomes possible even under very oxidizing conditions. This comparison illustrated the positive impact of the solid solution system on the stabilization of coffinite endmember. In model 3, we improved the mixing behavior of U and Th in uranothorite by using the subregular solid solution model (Eq. 3). The formation of uranothorite was found to move down in temperature even at 25°C and log(fO 2 ) < −30 (Fig. 9c). The stability field of uranothorite with relatively high U content expands as temperature increases (Fig. 9c, f), which eventually covers the whole oxygen fugacity range at 50°C, suggesting the favorable stabilization of U-bearing uranothorites. These calculations illustrate that thorite overall can stabilize coffinite by absorbing U into the matrix to form uranothorite solid solutions under the conditions where pure coffinite is unstable (Fig. 9d, e). Furthermore, this work also demonstrates the significant contribution of the non-ideal mixing effect to the formation and stabilization of high U content uranothorites in much wider environmental conditions, which suggest the weathering or material degradation processes of uranothorite to be thermodynamically unfavorable.
Although the thermodynamics favors the stabilization of uranothorite over a wide range of oxygen fugacity and temperature ranges, uranothorite may still undergo material degradation due to kinetic or chemical processes. According to the discussion in section "Solid solution model of uranothorite", uranothorite exhibits nanoscale heterogeneity, which controls chemical reactions and processes occurring at the interfaces or surfaces that may promote selective dissolution of U and Th. For instance, U can be involved in redox reactions leading to the more mobile uranyl-based species, but Th remains tetravalent, so the dissolution rate of uranothorite could be enhanced due to the nanoscale heterogeneity in the material. Plus, these small sized heterogeneities could also yield a significant amount of highly active sites in the interfaces or surfaces, which can act as catalysts for reactions with the aqueous solution to promote dissolution [92][93][94] . As a consequence, the heterogeneity can potentially introduce the incongruent corrosion process and induce intragranular stress in the bulk structure due to dislocation activity. This may also partially explain that why in nature,~35 mol% U were mostly documented for uranothorite minerals that are also finegrained, which could be in fact the degraded product of more U-rich uranothorites after selective dissolution of U upon the contact of hydrothermal fluids and/or the sequential recrystallization 95 . Such conclusion may also be more general for other actinide solid Estimated from this work and c p model is discussed in Supplementary Note 6. b The Gibbs free energies of formation were obtained by using the enthalpies of formation in this table, the standard entropies, and auxiliary thermodynamic parameters 12,40,49,85,120 . solutions, particularly for those when Pu is mixed with U or with multi-elements in the zircon structure. Although Pu, due to its even stronger covalent character [96][97][98] , may be more strongly stabilized in zircon due to potential thermodynamic non-ideality when mixed with U, Th, and/or minor actinides, the possibility of selective dissolution of Pu due to nanoscale disorder heterogeneous should not be overlooked. In summary, we provide an updated insight of the atomic-scale arrangement of U and Th cations in uranothorite. Despite being a non-ideal mixture, mixing of U and Th in uranothorite is suggested, based on nPDF and ML algorithm, to follow a disordering model over the short range while being heterogeneous over the nanosized range. Local geometric distortions are driven by the electron transfer due to actinide covalency explain the thermodynamic non-ideality, specifically the negative deviations in volume and enthalpy from the Vegard's law. This discovery emphasizes the role of 5f electrons in local geometric distortion and the resultant thermodynamic anomalies when actinides (Th, U, Np, and Pu) mix with other elements (transition metals or lanthanides) in natural minerals or actinide nuclear waste forms. Lastly, the solid solution model describing uranothorite was revised according to the structural findings to a subregular solid solution model, which was further used in the multi-component system thermodynamic modeling to investigate the phase equilibria of uranothorite under various redox-temperature conditions. Due to the non-ideal mixing effect, the solid solution system was found to have an enhanced and extended stability field for forming uranothorite with higher U content where pure coffinite is unstable, suggesting its weathering or material degradation to be thermodynamically unfavorable. Thus, covalency-driven nanosized disorder heterogeneity and thermodynamic non-ideality were demonstrated to have important implications for actinidecontaining solid solution systems in general and their potential applications for immobilizing Pu and minor actinides.

Sample preparation
All the uranothorite samples were prepared through hydrothermal treatment, according to the procedure already described in literature 8,10,99 .
Mixtures of U(IV) and Th(IV) concentrated solutions prepared in hydrochloric acid were added dropwise to a slightly off-stoichiometric amount of Na 2 SiO 3 (≈5%) previously dissolved in water. The pH of each solution was adjusted to 11.0-11.5 using 8 M NaOH, then NaHCO 3 was finally added to buffer the solutions at pH = 8.7. The final mixtures were transferred into 23 mL Teflon containers placed in autoclaves (Parr-type instruments). Hydrothermal treatments were performed at 250°C for 1 week. The obtained precipitates were separated from the solutions by centrifugation at 4000 r.p.m., washed with water then ethanol, and finally dried overnight in air at room temperature. The final uranothorite samples contained amorphous silica and thorium-uranium dioxides as impurities. In order to proceed to their purification, a procedure involving successive washing steps in static and dynamic conditions was successfully developed. The complete elimination of impurities was reached in a reasonable time with respect to the chemical durability of the uranothorite solid solutions 10 . Chemical compositions of the samples were checked by scanning electron microscopy-wavelength dispersive spectroscopy (SEM-WDS), and those results were published previously 12 . Table 1 provides the naming convention for the sample compositions which were analyzed as part of this work.

ND and scattering
Samples were measured at the Nanoscale Ordered Materials Diffractometer (NOMAD), Spallation Neutron Source (SNS), and Oak Ridge National laboratory (ORNL) 100 . Prior to the measurements, the samples were annealed at 500°C for 12 h according to the thermogravimetric studies on these samples 12 . Samples were then contained in boron free silica NMR tubes (Wilmad) with an inner diameter of 4.2065 ± 0.0065 mm and a wallthickness of 0.38 mm. A diamond standard was used to calibrate the conversion of time-of-flight to Q for each pixel 100 . A vanadium rod with diameter 5.85 mm was measured for data normalization, taking into account effects such as the incident neutron spectrum and the spectral efficiency of the detectors. An identical empty NMR tube is measured as a background for the sample measurements. All sample measurements were collected for an accelerator proton charge of 4 C, which at the time of the experiment corresponded to roughly 75 min. The PDFs have been calculated with Q max = 31 Å −1 leading to a real-space resolution of 0.02 Å.

GSAS refinement
The obtained ND data were analyzed by the Rietveld method using General Structure Analysis System (GSAS) software 101 . The background was fitted by a Chebyshev function, and the peak profiles were fitted to pseudo-Voigt convolution functions with a peak asymmetry correction 102 . The starting structure parameters for the uranothorite refinement were  Table 9. taken from work done by Fuchs et al. 103 . Data from all detector banks were simultaneously refined for the lattice parameters, atomic coordinates, and isotropic atomic displacement factors of U, Th, Si, and O. The detailed Rietveld refinement procedures were also stated previously elsewhere [104][105][106][107] .

Real-space fitting
A PDF is the probability that an atomic correlation exists at a corresponding radial distance. The neutron PDF (nPDF) was fitted using a combination of PDFgui and RMCProfile software tools 108,109 . The first silicon-oxygen bond distance was truncated in order avoid fitting the contribution of the silica NMR tube due to the limited amounts of samples. Real-space Rietveld refinement was performed on neutron total scattering data by implementing PDFgui 109 . G PDF (r) data from the Fourier transformation is fitted in PDFgui, and its detailed mathematical form can be found in Supplementary Information. The higher weighted R in Table 1 when compared with real-space Rietveld refinements to GSAS Rietveld refinements is likely a function of the background intensity which can significantly influence the weighted R 110 .

RMC modeling
RMC modeling was performed on the same neutron total scattering data. The software implemented for RMC modeling was RMCProfile 108 . G RMC (r) and F(Q) data were used for RMC modeling. For RMC modeling through RMCProfile, the nPDF was normalized according to the single atom scattering definition and will be denoted as G RMC (r) in this work, which differs from G PDF (r) used in PDFgui by a mathematical relationship showed in Supplementary Information. The starting atomic configurations for RMC fitting were obtained by creating a unit cell comprised of U, Si, and O using the atomic positions and lattice parameters obtained through PDFgui fitting (thereby avoiding the need to estimate a mass density), arbitrarily substituting Th for U to match the desired composition, then extending a unit cell comprised of U, Th, Si, and O to a 10 × 10 × 10 supercell. The swapping of U and Th during RMC fitting was implemented to relax the structure from the ordering introduced into the initial configuration. RMC modeling of samples U11 and U35 was carried out using distance window constraints and BVS constraints applied in order to confine silicon, uranium, and thorium atoms within oxygen polyhedra. This helps avoid the formation of interstitial atoms of Si, U, or Th and oxygen atoms without charge compensation. During the RMC modeling of samples U35 and U11, 12,600,000 attempted moves with 10,100,000 attempted swaps of U and Th (swap rate of 0.8) were carried out. For sample U84, a reasonable fit could only be obtained by significantly reducing the duration of the fit to 8,800,000 attempted moves with 7,080,000 attempted swaps. It is worth noting for reference that a swap rate of 0.5 indicates that the number of attempted atom swaps will equal the number of accepted atom displacements. The r-grid spacing was 0.02 Å and the maximum atomic moves were restricted 0.02 Å. The weighing ratio [F(Q): G RMC (r)] assigned to data during RMC fitting was 100:7.

Evaluation of cluster sizes through unsupervised ML
A K-means cluster algorithm was used to group the atoms in the atom boxes. An empirical relationship commonly referred to as the "inertia" parameter was used to ascertain the appropriate number of clusters of varying atom populations 111,112 . K-means clustering is one of the unsupervised learning algorithms for cluster analysis, and the distribution from Sci-kit was used for this work 111 , where more information can be found in the source literature. In K-means clustering, the algorithm partitions a set of available data (in this case the Euclidean atomic coordinates) into a predefined number of non-overlapping groups (where K is the predefined number). The clusters defined by the algorithm are selected such that they reduce an "inertia" parameter which is ∑(x i −c j ) 2 , where x i is the position of the ith atom and c j is the position of the jth centroid. Empirically, the number of K clusters that best describes the dataset can be approximated by the "shoulder method" where the eponymous "shoulder" is the point at which the slope of the curve approaches zero on the graph of inertia vs. number of clusters. Selecting too few clusters results in high inertia values, suggesting the dataset is under-fit while selecting too many clusters does not significantly affect the inertia, indicating the dataset is over-fit. For this study, the value corresponding to the shoulder was K = 50 clusters. The clustering algorithm was performed using U atoms from the RMCProfile atom models of samples U11 and U35, and Th atoms for sample U84. Clustering calculations were performed for each RMC atom model, where the starting partitions were selected from 40 configurations prior to refinement. Histograms were generated based on the averages of five clustering trials.

Computational methods
DFT + U calculations were performed by using Vienna Ab initio Simulation Package (VASP) 5.4.4 (ref. 113 ). Generalized gradient approximation (GGA) approximation was used for the exchange-correlation function, and projector-augment wave (PAW) was used to orient the effect of core electrons 114,115 . DFT-D3 correction with Becke-Johnson damping was added to interpret van de Waals interaction in the system 116 . Effective U parameters for 5f electrons in uranium and thorium was set as 4.0 eV to correctly describe the highly localized f electrons 117,118 . The valence electrons of uranium, thorium, silicon, and oxygen were respectively set as 5f 3 6s 2 6p 6 6d 1 7s 2 , 6s 2 6p 6 6d 2 7s 2 , 3s 2 3p 2 , and 2s 2 2p 4 . In ionic and electronic optimization steps, total energy was converged to 1 × 10 −5 eV and force between ions was converged to 1 × 10 −4 eV/Å. Based on the convergence study, integration of the Brillion zones of coffinite and thorite were done in a Gamma grid of 5 × 5 × 5 and cutoff energy was set as 550 eV. Considering the magnetic structure of USiO 4 , both ferromagnetic (FM) and colinear 1−k AFM phases have been tested in this study.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon request or can be accessed at https://doi.org/10.17632/wp87pcpx6h.1.