Computational Modeling of Cation Diffusion in Isolated Nanocrystals of Mixed Uranium, Plutonium and Thorium Dioxides

1. Introduction

The most widely used fuel for nuclear reactors at present is uranium dioxide UO₂ enriched in the fissile isotope U-235. Unfortunately, the available reserves of uranium are limited, which forces us to consider alternative fuels. The use of mixed MOX fuel (U,Pu)O₂ containing the fissile isotope Pu-239 is considered promising [1,2]. In the future, it is possible to use both uranium and plutonium in combination with thorium, since thorium-232 is the raw isotope for the production of uranium-233 [3]. In the future, uranium, plutonium and thorium dioxides are of interest as fuel for high-temperature gas-cooled reactors, fast reactors with lead, sodium and helium coolants, and even space power reactors [3,4,5,6].

The main reason for the possible degradation of fuel characteristics under the influence of neutron irradiation are radiation-stimulated phenomena, including the destruction of the crystal lattice by ballistic collision cascades, the accumulation of radiogenic gases and other nuclear decay products, changes in the stoichiometry, structure and grain size of the fuel crystals [7,8,9,10,11,12,13,14,15]. Oxygen diffusion in oxide fuel is a relatively fast process. However, oxygen migration is not sufficient for mass transfer. The rates of processes associated with transport of mass (grain growth, sintering, creep, plastic deformation, recrystallization) are determined by the much slower diffusion of uranium, plutonium and thorium cations [13,14].

In mixed uranium-plutonium and uranium-plutonium-thorium fuel, in the presence of temperature gradients or radiation damage, the difference in the diffusion coefficients of uranium, plutonium and thorium cations may act as a mechanism for the spatial separation of these elements, creating undesirable heterogeneity of composition. Thus, a comparison of the diffusion coefficients of different types of cations in the mixed oxides (U,Pu,Th)O₂ is of interest.

Obtaining experimental data to characterize the disordering of irradiated nuclear fuel is difficult due to the complexity of the defect structure and extreme conditions of the reactor core. Thus, there are a large number of works devoted to computational modeling of defect formation in UO₂, PuO₂ and ThO₂ [10,12,16,17,18,19,20,21,22,23,24,25,26,27]. The studied features include thermo physical, mechanical and electrical properties [15,16,17], the electron subsystem [19,20,21,22], disordering mechanisms, collision cascades and thermal spikes [28,29,30,31,32,33], phase transitions [3,34,35], crack propagation [18], non-stoichiometry effects [36]. The research methods were first-principles calculations, as well as static and dynamic modeling using both empirical and a priori interaction potentials [10,22,23,24,37,38,39].

Until recently, molecular dynamics (MD) modeling of the diffusion of intrinsic ions in (U,Pu,Th)O₂ crystals was limited to the anion sublattice or highly disordered systems due to the low mobility of cations. Cation migration was studied near the surface, cavities, grain boundaries, or in the presence of artificially created defects [40,41,42]. In [43], MD modeling of uranium diffusion in a defect-free bulk of UO₂ nanocrystals with a free surface serving as a source of cation vacancies was performed. A recent study [44] applied a similar approach to studying the diffusion of thorium cations in the bulk of ThO₂ nanocrystals. In the present work, the technique of [44] is used to model cation diffusion in mixed oxides (U_x,Pu_y,Th_1-x-y)O₂.

2. Modeling Methodology

In the present work, the model systems for studying the diffusion of intrinsic cations of thorium, plutonium, and uranium were isolated nanocrystals of ThO₂, PuO₂ and crystallites of mixed compositions (U_0.25Pu_0.5Th_0.25)O₂, (Pu_0.25Th_0.75)O₂ (Figure 1). The nanocrystals contained from 5460 to 43848 particles. At the beginning of the computational experiments, the model crystallites had the shape of regular octahedrons, the faces of which were formed by crystallographic planes of the (111) type (Figure 1a). During the simulation, the crystallites evolved into a more optimal shape of octahedrons with truncated vertices (Figure 1b).

The interaction of particles was described by a set of pair potentials MOX-07 [10], which was supplemented by a compatible thorium-oxygen potential proposed later in [46]. The use of pair interaction potentials somewhat reduces the accuracy of the modeling compared to the multiparticle potentials Cooper-2014 (CRG, [47]), increasing the simulation performance, which was critical for the present work. In [30], it was shown that the MOX-07 and CRG potentials give close energies of intrinsic disordering of UO₂ and PuO₂ crystals. In addition, the MOX-07 and [46] potentials made it possible to quantitatively reproduce the superionic transition temperatures in UO₂ [12] and ThO₂ [44] crystals.

Potentials MOX-07 account for the Coulomb interaction of effective charges of the particles, the valence repulsion of overlapping electron shells, and the dispersion attraction in the form Equation (1). The parameters of these potentials are given in [10].

U_ij(R_ij) = K_E⋅q_iq_j / R_ij + A_ij⋅exp(−B_ij) − C_ij/R_ij⁶.

(1)

The presence of a free surface ensured the natural influx of cation vacancies into the crystal bulk, which provided equilibrium concentration of the thermal Schottky defects as implementation of the vacancy mechanism of cation diffusion.

As the temperature decreased, the rate of generation of the thermal vacancies became too low for direct molecular dynamics simulation, so the temperatures in the present work were limited to the superionic state of the crystals that is close to melting.

The objective of this work was to calculate the diffusion coefficient of cations within the volume of the model nanocrystals. Near the surface, the diffusion coefficients are increased compared to the bulk. To exclude the surface effect, the diffusion of cations in this work was studied inside a central spherical region, the boundary of which was located at a distance of 1.75⋅a (a denotes the lattice constant) from the nearest surface (Figure 2).

During the computational experiments, the mean squares of the displacements were calculated, which were averaged over 20 computational experiments at each temperature. The linear sections of the dependences <a²(t)> on time t were used to calculate the diffusion coefficient D using the formula

<a²(t)> = 6Dt.

(2)

The diffusion of cations was studied in the course of computational experiments lasting up to 450 ns with a step of integration of the equations of motion of ∆t = 3∙10^–15 s. The equations of motion were integrated using the “leapfrog” method in a form mathematically equivalent to the Verlet integration [48].

Accumulation of computational error could lead to the emergence of translational motion of the crystallite, as well as to rotation around the mass center. Such motion would create an artificial increase in the mean square of the particle displacement, which would overestimate the diffusion coefficients calculated from Equation (2). In order to exclude non-diffusion movement of cations, the mass center of the model systems was fixed at the origin of coordinates. At each step of the MD simulation, the solid-state rotation of the crystal was compensated. Despite this, at simulation times of hundreds of nanoseconds, a significant rotation of the cation sublattice accumulated. This rotation was excluded when processing the simulation results using a special original procedure. An example of the dependences <a²(t)> obtained before and after the rotation compensation is shown in Figure 3.

During a computational experiment, the temperatures were kept constant using the Berendsen thermostat [49]. This thermostat was chosen for its computational efficiency and because it had the lowest temperature fluctuation amplitude in our calculations. However, it is not ergodic.

The Berendsen thermostat adjusts the temperature T so that deviations ΔT from the required value T₀ decrease with time t according to the exponential law ΔT ~ exp(-t/τ), where τ is the thermostat parameter. With increasing parameter τ, the particle energy distribution approaches ergodicity. In [50], it was demonstrated that for τ > 200⋅∆t, where ∆t is the integration time step in molecular dynamics simulations, the Berendsen thermostat achieves a quality comparable to that of the Nosé-Hoover thermostat [51,52]. In the present work, τ = 3000⋅∆t = 9 ps.

To verify the use of the Berendsen thermostat, we obtained cation diffusion coefficients in the (Pu_0.25Th_0.75)O₂ crystal of 5460 particles while stabilizing temperature using the stochastic velocity rescaling method (SVR, [50]), which ensures ergodic thermostating. The temperature regulation intensity of the stochastic velocity rescaling method follows the same exponential law, ΔT ~ exp(-t/τ), as in the case of the Berendsen thermostat. We maintained the value of τ at 3000⋅∆t = 9 ps for the stochastic velocity rescaling.

Figure 4 presents the averaged mean square displacements of thorium cations and oxygen ions over 20 computational experiments at a temperature of T = 3350 K. The difference in the slopes of the two lines falls within the margin of error observed in this work.

We note that the diffusion coefficients of different cation types within the same crystallite showed nearly identical values in our simulations. Figure 5 presents the temperature dependence of the combined diffusion coefficient for thorium and plutonium cations in (Pu_0.25Th_0.75)O₂ crystals, calculated using both Berendsen and stochastic velocity rescaling thermostats.

As evident from the plots in Figure 5, the linear fits approximating these dependencies in ln(D) = f(1/kT) coordinates are virtually identical. Therefore, we conclude that the use of the Berendsen thermostat for subsequent calculations in this work is well justified.

The cation diffusion coefficients in this work depended on the melting temperature of the model crystallites, which was determined by their composition. Moreover, to properly compare diffusion coefficients, it was necessary to consider the melting temperatures (T_∞) of infinite-size systems. The T_∞ values were estimated by extrapolating the T_melt(N^−1/3) dependence to N^−1/3 = 0, where N represents the number of particles in the nanocrystal.

The T_melt(N^−1/3) values were determined as the temperatures at which melting of the model crystallites initiated within 3 ns after the start of computational experiments in 90% of cases. As demonstrated for UO₂ in our previous work [53], this approach enables accurate determination of melting temperatures in macroscopic crystals. The T_melt(N^−1/3) relationships in this work showed nonlinear behavior, similar to Ref. [53], though they became linear at large system sizes (Figure 6).

The T_melt(N^−1/3) dependencies in our study exhibited nonlinear behavior, consistent with the results reported in [53]. However, in the large-size limit (R→∞), the relationships converged to the linear asymptotic form

T_melt(R) = T_∞⋅(1−R₀/R) = T_∞⋅(1−(N₀/N)^1/3),

(3)

where T_∞ represents the bulk melting temperature and R is the characteristic crystal size, scaling as N^1/3. The size parameters R₀ and N₀ correspond to a nanocrystal that would melt at zero temperature. We note that analytical models describing nanocrystal melting (see, e.g., Ref. [54]) yield dependencies of the form given by Equation (3) in the large-size limit.

In order to estimate the T_∞ values, we extrapolated the T_melt(N^−1/3) dependences to the N^−1/3 → 0 limit. This was achieved by applying Equation (3) to the three highest T_melt data points corresponding to the largest crystal sizes in our study (systems comprising 15960, 27600, and 43848 particles), as illustrated in Figure 6.

The temperature ranges investigated in this work for each crystallite composition extended from T ≈ 0.8⋅T_∞ up to the melting points of the model crystallites. These ranges corresponded to the superionic state of these systems. Lower temperatures were inaccessible for direct simulation due to low cation mobility.

The required calculation performance in this work was achieved by using an in-house code with parallelization of calculations on CUDA architecture graphics processors. The same code was used earlier in [44,46].

3. Results and Discussion

3.1. Cation Diffusion Coefficients

In all the computational experiments, the diffusion coefficients of cations of all types (uranium, plutonium, thorium) practically coincided when these ions were part of the same crystallite. This result corresponds to the experimental data described in [14]. According to these data, the temperature dependences of the diffusion coefficients of uranium, plutonium and thorium ions in the ThO₂ matrix coincided. In the context of our work, this coincidence can be explained by the fact that the diffusion movement of cations consisted of their collective movement along the chain towards the vacancies entering the crystal volume from the surface. On the other hand, the diffusion coefficients of cations depended on the composition of (U_xPu_yTh_1−x−y)O₂. This change correlated with the melting temperature, which increased with the thorium content. At a given temperature, the model nanocrystals with a lower melting point were characterized by higher cation diffusion coefficients (Figure 7).

The correlation between melting temperature and cation diffusion coefficients in structural analogs of (U,Pu,Th)O₂ crystals is well established. As demonstrated in [14], plotting these coefficients against reduced temperature (T/T_∞, with T_∞ being the macroscopic melting point) reveals nearly identical curves for UO₂ and CaF₂ crystals.

In the present work, the existence of the discussed relationship was verified for mixed oxides (U_xPu_yTh_1−x−y)O₂. It is evident from Figure 7 that the size of the model crystallites did not exert a clear effect on the diffusion coefficients of cations. Therefore, to obtain the reduced temperatures T/T_∞, we used the T_∞ values calculated by extrapolating the dependence T_melt on model crystallite size N to an infinite size, as described in the Modeling Methodology section.

The obtained values T_∞ of the melting point of the (U_xPu_yTh_1−x−y)O₂ crystallites are given in Table 1. For comparison, the experimental melting points of pure thorium and plutonium dioxides are also shown. It is evident that the model melting temperatures in the present work are overestimated in comparison with the experiment. However, the tendency for the melting temperature to increase upon transition from pure ThO₂ to pure PuO₂ is reproduced correctly.

The T_∞ values from Table 1 are used to build the dependences of the diffusion coefficients on temperature in the ln(D(T_∞/T)) coordinates that are shown in Figure 8. It is evident that this choice of coordinates brings all the ln(D) values to one region near a line, which is located close to the linear extrapolations of low-temperature experimental data known for uranium and thorium dioxides [14]. The experimental Arrhenius plots ln(D) vs (T_∞/T) for UO₂ and ThO₂ demonstrate near-identical behavior. It can be concluded that the differences in the diffusion coefficients of cations in (U_xPu_yTh_1−x−y)O₂ crystals are indeed associated with melting temperatures.

The melting of nanocrystals in the present study began with the formation of a surface melt, which subsequently propagated into the bulk. In the previous work [53] it was shown that at temperatures close to melting, a molten layer was formed on the surface of model UO₂ crystals, the proportion of molecules in which increased as melting approached. A surface disordered in this way can serve as a source of vacancies diffusing into the bulk of the crystal. Accordingly, the relationship between the bulk diffusion coefficients and melting temperatures obtained in the present work can be a consequence of the existence of the surface melt.

To investigate the potential relationship between cation diffusion coefficients and crystal characteristics influencing the melting process, we calculated the enthalpy of fusion L and the average surface energy γ for (U_xPu_yTh_1−x−y)O₂ crystallites consisting of 43848 particles. The obtained values are listed in Table 2. For each composition, these values correspond to temperatures reduced by approximately 100 K relative to T_∞.

Within the framework of our modeling, there was a decrease in the heat of fusion L upon the transition from PuO₂ to ThO₂, and the surface energy γ decreased in sequence ThO₂ → (Th_0.75Pu_0.25)O₂ → (Th_0.25U_0.25Pu_0.5)O₂, as shown in Table 2. The calculated heat of fusion of PuO₂ was close to the experimental values known for plutonium dioxide (66.5 ± 5.4 kJ/mol [57]) and uranium dioxide (70 ± 4 kJ/mol [58]). The correlation between these trends and the changing cation diffusivity will be the subject of future investigations.

3.2. Activation Energies of Cation Diffusion

The mechanism of cation diffusion in this work was their movement towards cation vacancies that migrated from the surface through the bulk of the model crystallites. The obtained effective activation energies for cation diffusion in crystals of different compositions (Figure 7) coincided within the error margin determined by the root-mean-square displacement fluctuations observed in the simulations. The values of these energies (E_D) were in the range from (10.0 ± 0.6) eV to (10.6 ± 0.7) eV. These values are significantly overestimated compared to low-temperature experimental recommendations for thorium and uranium dioxides (6.5 and 5.6 eV, respectively [14]).

The energies obtained can be interpreted using the well-known thermodynamic theory of Lydiard and Matzke [14,59]. Below, we designate cation vacancies as V⁴⁻, and anion vacancies as V²⁺. In the ideal solution limit, the concentration of neutral vacancy complexes V⁴⁻⋅2V²⁺ relative to cation sublattice sites follows the relation

[V⁴⁻⋅2V²⁺] = [V⁴⁻⋅2V²⁺]₀⋅exp{−E(V⁴⁻⋅2V²⁺)/kT},

(4)

where E(V⁴⁻⋅2V²⁺) represents the formation energy of the complex, and [V⁴⁻⋅2V²⁺]₀ denotes a temperature-independent pre-exponential factor. Molecular dynamics simulations should yield an effective formation energy matching E(V⁴⁻⋅2V²⁺).

The normalized concentrations of single cation and anion vacancies are related as

[V⁴⁻] = const⋅[V²⁺]⁻²⋅exp{−E_Sh/kT}.

(5)

Here, E_Sh denotes the formation energy of a classical Schottky defect, which is a non-interacting pair of cation vacancy and two anion vacancies. In the superionic state, the anion vacancy concentration becomes temperature-independent as the anion sublattice attains maximal configurational entropy, reaching the limit of possible disorder. Consequently, the effective formation energy of isolated vacancies equals E_Sh. On the other hand, at temperatures below the superionic transition

[V²⁺]² = [V²⁺]₀²⋅exp{−E_AF/kT},

(6)

where E_AF is the energy of anti-Frenkel disordering of the anion sublattice. Thus, at low temperatures, which correspond to the experimental data, the effective energy of formation of single cation vacancies is reduced compared to E_Sh to the value

E(V⁴⁻) = E_Sh − E_AF.

(7)

Furthermore, we consider the possible formation of an electrically neutral combination of defects V⁴⁻⋅V²⁺ + V²⁺, consisting of a charged vacancy cluster V⁴⁻⋅V²⁺ paired with a distant compensating anion vacancy (V²⁺). In the superionic phase, the formation energy of this complex, E(V⁴⁻⋅V²⁺ + V²⁺), equals the sum of E(V⁴⁻⋅2V²⁺) and the work needed to separate one anion vacancy to infinity. In the normal crystalline state, the effective energy of formation of the complex V⁴⁻⋅V²⁺ + V²⁺ decreases to the value

E_eff(V⁴⁻⋅V²⁺ + V²⁺) = E(V⁴⁻⋅V²⁺ + V²⁺) − E_AF/2.

(8)

The coefficient of vacancy diffusion of cations is proportional to the concentration of cation vacancies [V⁴⁻] or clusters [V⁴⁻⋅2V²⁺], [V⁴⁻⋅V²⁺] involved in the migration process. Thus, the effective activation energy of diffusion is the sum of the formation energy of the indicated defects and the height of the potential barrier E_M separating the equilibrium positions of the cation.

Table 3 presents the values of the intrinsic disorder energies of UO₂, PuO₂, and ThO₂ crystals calculated by static methods. Results from prior studies employing the MOX-07 potentials [10] (identical to those used herein) are labeled as “MOX-07”. The values calculated using the Cooper-2014 potentials proposed by M.W.D. Cooper, M.J.D. Rushton, and R.W. Grimes [47] are designated with the abbreviation CRG. The energies calculated with the MOX-07 and CRG potentials are close to each other, being overestimated relative to the experimental values characterizing the cation sublattice.

The overestimation of the cation sublattice disordering energies is characteristic of the known classical interaction potentials for UO₂. The energies of Schottky trio formation calculated using such potentials are in the range from 7.7 to 14.7 eV [43], while the experimental data are in the range from 6 to 7 eV [14]. This systematic difference indicates changes in the electronic state of cations that are not taken into account by classical potentials. For comparison, Table 3 also shows the values of disordering energies obtained by first-principles calculations within the framework of the density functional theory (DFT) with the local density approximation (LDA) or the generalized gradient approximation (GGA) [59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,81]. Experimental values are marked with the abbreviation “exp“.

It is known that the Hubbard-U correction can be used to correctly calculate the band gap in crystals with valence f-electrons (UO₂, PuO₂). In [60,66], the use of DFT+U led to an increase in the calculated disordering energies of UO₂ compared to the LDA and GGA approximations without the Hubbard correction. The authors of [66] used the occupation matrix control to avoid obtaining metastable states. On the other hand, in [59] the use of a similar matrix-controlled technique in LDA+U and GGA+U calculations yielded relatively low values of the energies of anti-Frenkel disordering (3.3 eV) and the formation of a classical Schottky defect (6.0 eV).

The analysis of the literature data did not allow us to reveal a systematic difference in the disordering energies of (U,Pu,Th)O₂ crystals corresponding to different DFT approximations. However, there is a group of results in which the formation energies of identical defects in UO₂, PuO₂ and ThO₂ are close to each other. We believe that the similarity of the disordering energies of all three oxides is in agreement with the empirical data, since these oxides behave similarly in diffusion experiments [14].

As previously discussed, the temperature-dependent diffusion coefficients of oxygen and cations in UO₂, PuO₂, and ThO₂ crystals exhibit converging behavior when plotted in ln(D) versus T_∞/T coordinates. Exact coincidence of these curves would imply the following relationship between effective activation energies:

E_D/kT_∞(ThO₂) = E_D/kT_∞(UO₂) = E_D/kT_∞(PuO₂).

(9)

Using the experimentally established melting temperatures of 3120 K for UO₂ [75] and 3665 K for ThO₂ [55], along with the recommended uranium migration energy of 5.6 eV from [14], Equation (9) predicts a thorium migration energy in ThO₂ of 6.6 eV:

E_D(ThO₂) = T_∞(ThO₂)/T_∞(UO₂)⋅E_D(UO₂) = 6.6 eV.

(10)

This prediction is in quantitative agreement with the experimental value of 6.5 eV from the review [14].

Relation (9) may provide reasonable approximation for the intrinsic disorder energies. If valid, the values of E_AF and E_Sh should increase along the series PuO₂ → UO₂ → ThO₂, exhibiting a 20 % change from PuO₂ to ThO₂. Our static calculations satisfy this predicted trend.

Regarding the V⁴⁻⋅2V²⁺ cluster, we limited ourselves to considering the linear configuration <111>, which is a V²⁺−V⁴⁻−V²⁺ chain oriented along the [111] direction; the formation energies of other variations of this cluster are close to E(V⁴⁻⋅2V²⁺)₁₁₁ [30,64,72].

The cation migration energies E_M calculated both in the present and in other works are overestimated compared to the experimental recommendations (Table 3). There are calculations indicating a significant decrease in the migration energy of cation vacancies in the case of their movement as part of vacancy clusters. For example, in the work [65] the migration energy of a cation vacancy in uranium dioxide as part of the V⁴⁻⋅2V²⁺ cluster decreased compared to the migration energy of a single vacancy from 3.09 eV to 2.19 eV. On the other hand, the migration energy of thorium cations in ThO₂ (4.47 eV) considered by the authors of [65] in a similar manner, did not change.

In the present work, the lattice statics method was used to calculate the migration energies of a cation through a vacancy complex (V⁴⁻⋅2V²⁺)₁₁₁ in thorium and plutonium dioxides. The obtained values of E_M(V⁴⁻⋅2V²⁺) for ThO₂ and PuO₂ were 4.4 eV and 5.5 eV, respectively. These energies closely match the migration energies of single cation vacancies (4.5 eV for PuO₂ and 5.45 eV for ThO₂), which were also computed in the present work using the lattice statics approach.

According to Equations (5-8) and the calculation data for MOX-07 potentials from Table 3, the results of dynamic modeling obtained in this work can only correspond to the diffusion of cations via the bound cluster V⁴⁻⋅2V²⁺ with effective activation energies

E_D,Static(ThO₂) = E(V⁴⁻⋅2V²⁺)₁₁₁ + E_M(V⁴⁻⋅2V²⁺) = 6.9 + 5.5 = 12.4 eV,

(11)

E_D,Static(PuO₂) = E(V⁴⁻⋅2V²⁺)₁₁₁ + E_M(V⁴⁻⋅2V²⁺) = 5.0 + 4.4 = 9.4 eV.

(12)

The static calculation of the diffusion activation energy obtained for PuO₂ agrees well with molecular dynamics results. However, for thorium dioxide, the static prediction E_D,Static(ThO₂) = 12.4 eV significantly overestimates the dynamic activation energy E_D = 10.3 ± 0.3 eV (Figure 7). The discrepancy may originate from the collective cation motion involving temporary displacements of certain cations into interstitial sites. This complex cation diffusion mechanism in ThO₂ was previously identified in [44].

As temperature decreases, the superionic state of (U,Pu,Th)O₂ crystals transitions to a conventional crystalline phase containing either charged impurity-vacancy complexes V⁴⁻⋅V²⁺ or neutral clusters V⁴⁻⋅2V²⁺, as well as single cation vacancies. Notably, the effective diffusion activation energies derived from Table 3 data show no reduction in these cases, maintaining values comparable to those in the superionic regime. Assuming approximately equal migration energies for cations through both charged and neutral vacancy clusters, we obtain

E_D,Static(V⁴⁻⋅V²⁺ in ThO₂) = E(V⁴⁻⋅V²⁺ + V²⁺) + E_M(V⁴⁻⋅2V²⁺) − − E_AF/2 = 9.7 + 5.5 − 4.5/2= 12.9 eV;

(13)

E_D,Static(V⁴⁻⋅V²⁺ in PuO₂) = E(V⁴⁻⋅V²⁺ + V²⁺) + E_M(V⁴⁻⋅2V²⁺) − − E_AF/2 = 7.0 + 4.4 − 4.5/2 = 9.45 eV;

(14)

E_D,Static(single V⁴⁻ in ThO₂) = E_Sh − E_AF + E_M = 12.7 + 5.45 − 4.5 = 13.7 eV;

(15)

E_D,Static(single V⁴⁻ in PuO₂) = E_Sh − E_AF + E_M = 9.5 + 4.4 − 3.9 = 10.0 eV.

(16)

The static calculations presented above demonstrate that the superionic-to-normal crystalline phase transition is unlikely to reconcile the effective diffusion activation energies obtained from molecular dynamics simulations with low-temperature experimental data [14]. To improve agreement, it would be necessary to reduce both the calculated Schottky disorder energies and the cation vacancy migration energies.

4. Conclusions

According to the modeling performed in the present work, the diffusion coefficients of cations of all types in the bulk of the same (U_xPu_yTh_1-x-y)O₂ crystal should be very close. This result is consistent with the experimental data from [14] and can be explained by the vacancy migration mechanism, in which a group of different cations moves along a chain towards one vacancy (or vacancy cluster). At the same time, the diffusion coefficients in the present work changed with the melting temperature of the nanocrystals, which rose with increasing thorium content.

The cation diffusion coefficients exhibited an inverse correlation with the melting points of (U_xPu_yTh_1-x-y)O₂ crystallites at constant temperature. When plotted as ln(D) versus T_∞/T (where T_∞ represents the bulk melting temperature), the temperature dependencies for different compositions converged toward a common curve, demonstrating the universal scaling of diffusion behavior with melting point.

The values of the cation diffusion coefficients computed for superionic (UₓPuᵧTh₁₋ₓ₋ᵧ)O₂ nanocrystals agree with extrapolated low-temperature diffusion data from macroscopic UO₂ and ThO₂ crystals [14]. However, the model effective diffusion activation energies remain systematically overestimated relative to experimental values. Our analysis reveals no mechanism for reducing these activation energies during the transition from the superionic to normal crystalline phase. Improved agreement would require refining both the calculated Schottky defect formation energies and cation vacancy migration barriers.