Molecular Dynamics Simulation of Oxygen Diffusion in (Pu_xTh_1-x)O₂ Crystals

1. Introduction

Plutonium dioxide is used as a nuclear fuel in a mixture with uranium dioxide. The presence of the plutonium-239 isotope increases the average fission neutron yield, which facilitates the breeding of fissile isotopes. A similar application of PuO₂ in a mixture with thorium dioxide is also feasible; this would enable the production of the fissile uranium-233 isotope through neutron capture by thorium-232 [1]. PuO₂ and ThO₂ crystals share the same fluorite structure (Fm-3m) and have similar lattice parameters, thus allowing for the formation of a mixed oxide (Pu_xTh_1−x)O₂ [1].

Within a nuclear reactor core, fuel crystals are subjected to intense radiation damage. To model phenomena related to mass and charge transport under neutron irradiation, as well as to improve techniques for the electrochemical recovery of plutonium and thorium from irradiated fuel [2,3], it is essential to study the oxygen transport mechanisms in (Pu_xTh_1−x)O₂ oxides.

Oxygen diffusion in pure PuO₂ and ThO₂ oxides has been experimentally investigated [4,5,6,7]. However, interpreting the underlying diffusion mechanisms is challenging due to discrepancies in the reported data, which vary in both the values of the diffusion coefficients and the activation energies (in the case of ThO₂). In particular, the relative contribution of oxygen transport via interstitial sites versus anion vacancies remains unclear.

Near the melting point, PuO₂ and ThO₂ crystals may transition to a superionic state, characterized by an extreme disordering of the oxygen sublattice. A transitional region exists between this state and the "normal" crystalline phase [8], where the mechanisms of oxygen transport may change with temperature. The available experimental data pertain to lower temperatures and thus do not characterize this specific region. Therefore, computational modeling of oxygen diffusion in PuO₂, ThO₂, and (Pu_xTh_1−x)O₂ crystals is relevant.

Computational methods have been used to calculate the energies of Frenkel and Schottky disorder in PuO₂ and ThO₂ crystals, as well as the activation energies for oxygen transport under such disordered conditions [9,10,11,12,13]. These static calculations employed both first-principles approaches [9,10,11] and classical interatomic potentials [12,13]. Both first-principles and classical models yielded quantitative agreement between the calculated energies and experimental estimates. However, these static calculations were limited to relatively simple single-particle transport models, which may not accurately represent the actual processes occurring at high temperatures, particularly near the superionic transition.

Studies on the dynamic modeling of oxygen diffusion in crystalline nuclear oxide fuels have primarily focused on uranium dioxide, which is a structural analog of PuO₂ and ThO₂ [14,15,16,17]. In most of these works, periodic boundary conditions were used. The periodic boundary conditions eliminate the presence of surface. In contrast, the work [14] modeled isolated UO₂ nanocrystals with free surfaces, which led to an improved agreement between the calculated oxygen diffusion coefficient and experimental data. Further investigation into the influence of surfaces as a source of point defects involved in oxygen transport is therefore of significant interest. This has become one of the objectives of the present work.

The model systems in this work were isolated (Pu_xTh_1−x)O₂ nanocrystals. To the best of our knowledge, oxygen diffusion in such systems has not been previously modeled using the molecular dynamics method. A methodology analogous to the one employed here was used in earlier studies to investigate oxygen diffusion in hypostoichiometric UO_2−x crystals [18] and cation diffusion in ThO₂, (U_x,Pu_y,Th_1−x−y)O₂ [19,20]. In [19,20], the primary mechanism for cation migration was the movement of cation vacancies formed in the near-surface layer.

2. Modeling Methodology

Molecular dynamics simulations were used to calculate the oxygen diffusion coefficients in ThO₂, PuO₂, (Pu_0.25Th_0.75)O₂ and (Pu_0.5Th_0.5)O₂ crystals. The selection of these model system compositions allowed us to track the correlation between oxygen mobility and the plutonium-to-thorium ratio. The model systems were nanocrystals comprising 5460 and 15,960 particles, isolated in vacuum. These crystallites possessed a free surface, predominantly faceted by {111} crystallographic planes. The presence of the free surface enabled the migration of point defects (such as interstitial ions and vacancies) from the surface into the crystal bulk.

At the start of the simulation, the nanocrystals had the shape of perfect octahedra (Figure 1). During the computational experiments, the model systems transformed into octahedra with truncated vertices. This result is consistent with the findings of [21], where a truncated octahedron was the equilibrium shape for model UO₂ nanocrystals.

To maintain a symmetric charge distribution on the nanocrystal surface, the structures were constructed from neutral PuO₂ and ThO₂ molecules as building blocks. Within each molecule, the two oxygen ions were positioned collinearly with the metal cation. Their specific coordinates were defined by the displacements (−0.25, −0.25, −0.25) a and (0.25, 0.25, 0.25)·a relative to the cation, where a is the lattice constant. This approach minimized the influence of surface charges on the bulk properties of the model crystals, which is essential due to the long-range character of Coulomb forces.

Despite the presence of surfaces, the aim of this study was to calculate the oxygen diffusion coefficient representative of the crystal bulk. The calculation of the diffusion coefficient was confined to a central region situated no less than 1.75 lattice constants from the surface to avoid the influence of surface migration (Figure 1). The evolution of the surface morphology during the computational experiments did not affect the central region of the crystals.

The interatomic interactions in the model crystals were described using the MOX-07 pair potentials (Potashnikov et al. [23]), which were supplemented with a compatible Th-O pair potential proposed later in the work [24]. In previous studies [14,16], the MOX-07 potentials yielded a quantitative agreement for the oxygen diffusion coefficient in UO₂ with experimental data [4,5]. Furthermore, the work [17] demonstrated that the MOX-07 potentials performed well in modeling radiation damage in (U,Pu)O₂.

A comparison of the MOX-07 potentials with the CRG many-body potentials proposed by Cooper et al. [25] was conducted in studies [12,17]. It was found that the MOX-07 and CRG potentials give close energies of intrinsic disordering of UO₂ and PuO₂ crystals. Moreover, the MOX-07 and [24] potentials have been shown to quantitatively predict the superionic transition temperatures for UO₂ [14] and ThO₂ [19], respectively.

The specified interatomic potentials were presented in the Born-Mayer form with dispersive attraction:

U_ij(R_ij) = K_E·q_iq_j/R_ij + A_ij·exp(−B_ij) − C_ij/R_ij⁶,

(1)

where K_E is Coulomb's constant; A_ij, B_ij and C_ij are the parameters describing both the repulsion between valence electron shells and the dispersion attraction. The parameter values are given in [23] and [24]. The effective charges for the cations and anions were set to q₊ = 2.745e, q_– = −1.3725e, respectively, where e is the elementary charge.

Molecular dynamics simulations were performed using original in-house software, previously employed in studies [20,26]. The particle equations of motion were integrated using the "leapfrog" scheme, which is mathematically equivalent to the Verlet integration [27]. The integration time step was set to ∆t = 3∙10⁻¹⁵ s. At each step, the center-of-mass displacement and rigid-body rotation of the entire crystal were compensated for.

At the beginning of each computational experiment, the particles were placed at the sites of a perfect crystal lattice. The components of the particle velocity vectors were assigned according to a Maxwell-Boltzmann distribution using the Box-Muller algorithm [28]. During the initial simulation stage, the system temperature increased due to the relaxation of the free surface. A Berendsen thermostat [29] with a time constant of τ = 200·∆t = 0.6 ps was used to bring the system to the target temperature.

For the production stage, during which the diffusion coefficient was calculated, the target temperature was maintained using a stochastic velocity rescaling thermostat [30] with a time constant of τ = 3000·∆t = 9 ps.

The diffusion coefficients D were calculated using the well-known relation

<a²(t)> = 6Dt,

(2)

where <a²(t)> is the mean squared displacement of the oxygen anions as a function of simulation time t. The total duration of the simulations varied from 1.5 ns to 120 ns, increasing as the temperature decreased. For statistical averaging, the displacements of all oxygen ions located within the central spherical region (Figure 1) were used. Furthermore, multiple independent simulation runs – from 10 to 50 for each temperature – were performed to improve the statistical accuracy of the results.

During long simulation times exceeding 20 ns, the model crystals could rotate from their initial positions despite the ongoing compensation for rigid-body rotation. These rotations led to a non-diffusive displacement of the oxygen ions, introducing an error into the calculated diffusion coefficient. This rotation was removed during the post-processing of the saved crystal configurations using a procedure detailed in [20]. This same procedure also enabled the tracking of individual diffusion jumps of specific anions. Consequently, it was used to analyze the underlying diffusion mechanisms, which are discussed in the Results and Discussion section.

To achieve the necessary computational performance, we employed an in-house code optimized for parallel execution on CUDA architecture Graphics Processing Units (GPUs).

The diffusion coefficient in PuO₂ crystals composed of 1500 particles with periodic boundary conditions (PBC) was calculated using the LAMMPS software package (Ver. 29August2024) [31] with a Nose-Hoover thermostat [32,33] and a time integration step of 3∙10⁻¹⁵ s. The LAMMPS simulations also employed parallel computing on GPUs.

To elucidate the mechanisms of oxygen migration, this work obtained the temperature dependence of the average concentration of anion interstitials in the central regions of model PuO₂ and ThO₂ nanocrystals consisting of 5460 particles. Interstitial anions were identified based on the fact that the Cartesian coordinates of the anion sublattice sites are multiples of half the lattice constant (0.5a), and these coordinates remained unchanged throughout the simulations. In PuO₂ and ThO₂ crystals, the cations are located in every other center of the cubes formed by the anion sublattice. The unoccupied centers of these cubes were defined as interstitial sites. Other types of interstitial positions, such as the centers of cube faces or edges, were not considered.

To identify interstitial ions, their coordinates were averaged over a sliding time window of 0.3 ps, which constituted approximately half of the characteristic time for an oxygen anion's diffusion jump between two equilibrium positions. A particle was considered to occupy an interstitial site if its time-averaged coordinates were closer to that interstitial position than to any node of the anion sublattice. The concentration of interstitial anions [O_I] was calculated as the ratio of their number to the number of anion sublattice sites. This concentration was then averaged over the entire simulation time. At sufficiently low temperatures, the appearance of even a single interstitial ion in the central region of the crystallite became a rare event. In such cases, the interstitial ion concentration was assumed to be proportional to the probability of detecting such an ion.

3. Results and Discussion

3.1. Oxygen Diffusion Coefficient as a Function of Model System Composition

The temperature-dependent oxygen diffusion coefficients calculated in this work for ThO₂, PuO₂, (Pu_0.25Th_0.75)O₂ and (Pu_0.5Th_0.5)O₂, nanocrystals are presented in Figure 2. The data are plotted in Arrhenius coordinates as ln(D) versus 1/kT, where k is the Boltzmann constant and T is the absolute temperature

In Figure 2 and throughout this paper, the abbreviation FBC refers to the free-surface boundary conditions employed in the simulations.

The lowest diffusion coefficients correspond to ThO₂. The oxygen diffusion coefficient increases with the plutonium fraction. Furthermore, increasing the plutonium content to 50% is sufficient for the oxygen mobility in (Pu_xTh_1−x)O₂ to reach values comparable to those in pure PuO₂.

The increase in the oxygen diffusion coefficient observed when transitioning from ThO₂ to PuO₂ correlates with a decrease in the melting temperature of the system. The experimental melting point of ThO₂ is 3665 K [40], whereas for plutonium dioxide it drops to 3017 K [41]. The model (Pu_xTh_1−x)O₂ crystals simulated using the MOX-07 potentials exhibit the same behavior, as demonstrated in a previous study [20].

3.2. High-Temperature Diffusion

The non-linearity of the lnD = f(1/kT) dependencies in Figure 2 is caused by changes in the oxygen transport mechanisms and/or the disordering of the oxygen sublattice with temperature. Values of lnD > −10 (D > 4.5 cm²/s) correspond to the superionic state of the system, where the oxygen sublattice resembles a melt. Consequently, the effective activation energies for oxygen diffusion E_D reach minimal values, close to 1 eV (Figure 3 and Figure 4). As the temperature decreases, the E_D values become linear combinations that include the formation energies of oxygen sublattice defects (either interstitial anions or anion vacancies) and the migration energies of these defects. Specifically, for the classical anti-Frenkel (AF) disorder of oxygen [4]:

E_D = E_M + E_AF/2,

(3)

where E_AF is the formation energy of an unbound interstitial anion–anion vacancy pair, and E_M is the migration energy, defined here as the height of the potential energy barrier separating the stable positions of the mobile defects, either interstitial anions or vacancies.

In order to discuss the mechanisms of oxygen diffusion, the lnD = f(1/kT) dependencies for PuO₂ and ThO₂ crystals are presented separately in Figure 3 and Figure 4. It can be seen that in both cases, as the temperature decreases, the superionic state gives way to a transitional region (Region I). The diffusion coefficients calculated using periodic and free-surface boundary conditions coincide in this region. This suggests that cation transport in this region is governed by bulk disorder within the crystallite, without significant influence from the surface.

Figure 3 and Figure 4 demonstrate the linear segments of the graphs, which are described by Equation (4)

ln(D) = −E_D/kT + const,

(4)

where E_D represents the effective activation energies for diffusion. In Region I, we obtained the following values for these energies:

-E_D(PuO₂, free-surface border conditions, 5460 particles) = 3.47 ± 0.05 eV;

-E_D(PuO₂, periodic border conditions, 1500 particles) = 4.2 ± 0.05 eV;

-E_D(ThO₂, free-surface border conditions, 5460 particles) = 5.24 ± 0.08 eV.

These values are overestimated compared to the experimental data presented in Figure 2, Figure 3 and Figure 4: E_D(PuO₂ [34]) = 1.8 eV, E_D(PuO2 [35]) = 1.95 eV, E_D(ThO₂ [37]) = 2.9 eV, E_D(ThO₂ [38]) = 2.2 eV, E_D(ThO₂ [36]) = 2.15 eV, E_D(ThO₂ [39a]) = 2.05 eV, E_D(ThO₂ [39b]) = 2.0 eV. However, the experimental data correspond to lower temperatures, where the oxygen transport mechanism may have changed. The present modeling demonstrates that such a change could indeed have occurred.

To clarify the diffusion mechanism, this work analyzed the trajectories of migrating anions. It was found that anion migration in Region I involved the collective motion of several particles. An example of such a process is illustrated in Figure 5 and Figure 6.

In the first stage of the migration (Figure 6a, purple arrows), ions #1526 and #1531 simultaneously move into adjacent sites of the anion sublattice. Concurrently, ion #1418 shifts into an interstitial position, while ion #878 moves into an anion vacancy.

The second stage of the migration corresponds to the 9.6 ps mark in Figure 5. At this stage, ions #1412 and #1426 move into adjacent interstitial sites, while ion #1406 begins to move towards a neighboring lattice site. In the third stage, ions #1412 and #1426 transition from their interstitial positions back into lattice sites. In this process, ion #1412 occupies the site vacated by ion #2290, which itself shifts into an adjacent vacant lattice site.

3.3. Ionic Conductivity Compared to Structural Analogs

It should be noted that some anions exchanged positions during the collective migration, such as ions #1412 and #2290 in Figure 5 and Figure 6. This exchange diffusion does not result in net charge transport. However, the movement of interstitial anions and anion vacancies in opposite directions was also observed. We propose that the diffusion mechanism described above for Region I is capable of supporting net charge transport. As indirect support for this hypothesis, Figure 7 compares the lnD = f(1/kT) dependencies obtained in this work with experimental data from [8] on the ionic conductivities of SrCl₂, CaF₂, SrF₂, and BaF₂ crystals. These materials share the fluorite structure with PuO₂ and ThO₂.

In order to compare the results of this work with the data from [8], we calculated the ionic conductivity of the model PuO₂ and ThO₂ crystals using the known formula

σ(T) = n·q²·D(T)/(kT),

(5)

charge, taken to be −1.3725e. For simplicity, it was assumed that charge transport could be mediated by any diffusive displacement of the anions that occurred during the simulation. According to Equation (5), the diffusion coefficient D(T) is effectively proportional to σ(T)·T, since the density n exhibits only a weak dependence on temperature.

The plots of σ(T)·T are presented in Figure 7, graphed as a function of the normalized coordinate T_m/T, where T_m is the melting temperature of the crystals. It is known that when using T_m/T coordinates, the transport coefficients for various fluorite-structure crystals are found to be similar [4]. The melting temperatures used for CaF₂, SrF₂, BaF₂, and SrCl₂ were T_m(CaF₂) = 1690 K, T_m(SrF₂) = 1740 K, T_m(BaF₂) = 1620 K, and T_m(SrCl₂) = 1148 K, as taken from [8].

All the ln(σ·T) curves in Figure 7 share a common nonlinear shape, indicative of a transitional region between the superionic state and the low-temperature crystalline phase. This region is characterized by an elevated effective activation energy for charge carrier diffusion. It can be inferred that the anion transport mechanisms in all systems considered in Figure 7 are similar. This supports the physical realism of the oxygen transport mechanisms identified in the present work. It is noteworthy that similar results were previously obtained in simulations of oxygen diffusion in uranium dioxide [14,16].

3.4. Interstitial Oxygen Anions

To analyze the mechanisms of oxygen transport, this work calculated the temperature dependence of the oxygen interstitial ion concentration [O_I] in PuO₂ and ThO₂ crystals, as shown in Figure 8. It can be seen that the shape of these dependencies closely resembles that of the diffusion coefficient graphs. The linear segments highlighted in Figure 8 can be described by the relations

ln(D) = −E_f,eff/kT + const,

(6)

where E_f,eff is the effective formation energy of an interstitial anion. Specifically, in the case of classical anti-Frenkel disorder, E_f,eff = E_AF/2, where E_AF is the energy required to form an unbound interstitial anion – anion vacancy pair within the crystal bulk.

The plots in Figure 8 demonstrate a transitional region between the superionic and crystalline phases, which corresponds to Region I discussed previously. Within this region, the effective formation energies of interstitial anions reach their maximum values:

- E_f,eff(PuO₂, I) = 2.52 ± 0.02 eV;

- E_f,eff(ThO₂, I) = 2.9 ± 0.2 eV.

These values are close to the experimental formation energy of an oxygen interstitial in uranium dioxide, E_f,eff(UO₂, exp) = 2.3 eV measured in [6].

The E_f,eff values obtained in Region I should correspond to the formation of interstitial anion – vacancy pairs within the crystal bulk, as the influence of the surface on oxygen diffusion in this Region is small. The formation energies of unbound Frenkel pairs E_AF, calculated for the MOX-07 potentials using the static lattice method, are E_AF(ThO₂) = 4.5 eV [19] and E_AF(PuO₂) = 3.9 eV [17,20]. Halving these values gives E_AF(ThO₂)/2 = 2.25 eV and E_AF(PuO₂)/2 = 1.95 eV, which are lower than the dynamic values from Region I.

The discrepancy between the dynamic E_f,eff values and the E_AF/2 values can be attributed to the formation of short-lived, close interstitial-vacancy pairs during the simulation. The recombination rate of such correlated pairs is not proportional to the product of the vacancy and interstitial concentrations. This leads to the disappearance of the 1/2 factor that links E_f,eff and E_AF in the classical model of Frenkel disorder.

Based on the plots shown in Figure 8, the oxygen migration energies in the model nanocrystals can be estimated as the difference between the effective activation energies for diffusion and the formation energies of interstitial anions:

E_M= E_D − E_f,eff.

(7)

This yields E_M(PuO2, I) = 3.47 − 2.52 = 0.95 eV and E_M(ThO2, I) = 5.24 − 2.9 = 2.34 eV. For plutonium dioxide, the E_M value is close to the experimental estimate of the potential barrier for interstitial diffusion in UO₂ (1.0 eV [5]). However, the significance of a direct comparison between these energies is questionable, since anion diffusion in Region I involved a complex collective motion of multiple particles.

3.5. Low-Temperature Region

As can be seen from Figure 2, Figure 3 and Figure 4, at oxygen diffusion coefficients below ~5·10⁻⁸ cm²/s (Region II), the slope of the lnD = f(1/kT) plots decreases significantly. Consequently, the effective activation energies for diffusion E_D are reduced. This change results in the low-temperature segment of the lnD(PuO₂) = f(1/kT) plot showing quantitative agreement with experimental data [34,35] (Figure 3). In the case of ThO₂, the agreement is substantially improved, with the calculated E_D = 3.15 eV aligning much more closely with experimental value [37] (2.9 eV). Figure 8 indicates that the primary reason for this change is a decrease in the effective formation energy of interstitial anions E_f,eff. In Region II, these energies become significantly lower than the corresponding E_AF/2 values.

The data obtained in this work reveal a fundamental divergence in the calculated oxygen diffusion coefficients for PuO₂ below 1800 K, depending on whether periodic (PBC) or free-surface (FBC) boundary conditions are used. The data presented in Figure 4 suggest that a similar discrepancy exists for ThO₂ below 2300 K. Therefore, the observed change in the E_f,eff values must be caused by the presence of the free surface.

Figure 9 illustrates a characteristic change in the coordinates of oxygen anions involved in a diffusion event within the bulk of a 5460-particle PuO₂ nanocrystal at T = 1500 K, which lies within the discussed low-temperature Region II. For brevity, only the motion along the Cartesian y-axis is shown; the other coordinates exhibited analogous changes. The values y = −a, −0.5·a, 0, and 0.5·a correspond to the sites of the anion sublattice, whereas y = −1.25·a, −0.75·a,−0.25·a, 0, and 0.25·a, correspond to interstitial positions.

The change in the y-coordinate shows that the anions moved in a chain, sequentially displacing each other into different interstitial positions – a process known as the interstitialcy mechanism. In this event, anion #1274, under the influence of ion #2287, moved into an adjacent lattice site very rapidly, although a brief displacement of ion #1274 into an interstitial position was still recorded.

The anion migration mechanism shown in Figure 10 is simplified compared to the collective particle motion characteristic of the high-temperature Region I (Figure 5 and Figure 6). Nevertheless, the diffusion is still facilitated by the presence of interstitial anions. We propose that in the low-temperature Region II, the primary source of these interstitial anions is the surface, or anti-Frenkel disorder within the near-surface layer. This would account for the observed decrease in the defect formation energy E_f,eff.

The migration energies of interstitial anions in Region II are calculated as follows:

- E_M(PuO₂, II) = 1.87 − 0.56 = 1.31 eV;

- E_M(ThO₂, II) = 3.15 − 1.22 = 1.93 eV.

These values are close to those in Region I, indicating the similarity of the anion diffusion mechanisms in both regions. Notably, the E_M(PuO₂, II) value for PuO₂ approaches the experimental estimate of 1.58 eV for the activation energy of low-temperature oxygen interstitial diffusion from Ref. [5].

Figure 3 and Figure 4 compare the simulations performed with free-surface boundary conditions for PuO₂ and ThO₂ nanocrystals consisting of 5460 and 15,960 particles. It can be seen that for the larger crystallites, the position of the kink in the lnD = f(1/kT) dependencies remains unchanged, and the diffusion activation energies E_D are also unaffected.

This result indicates that the formation of Region II may be independent of the crystallite size, a finding that requires further investigation. If this holds true, a similar oxygen transport mechanism could be applicable to real crystals, where the role of the surface may be played by grain boundaries or other extended defects.

It should be noted that in Region II, we did not observe the migration of anion vacancies through the internal region of the model nanocrystals. The methodology for processing the simulation data allowed for detecting the vacancy migration, since the movement of cation vacancies was detected in a previous study [20] investigating cation diffusion in (U_xPu_yTh_1−x−y)O₂ nanocrystals (at temperatures above 2700 K). We can hypothesize that at low temperatures anion vacancies are bound to the nanocrystal surface. Consequently, their migration into the central region is a rarer event than the appearance of interstitial anions.

4. Conclusions

In this work, oxygen diffusion in (Pu_xTh_1−x)O₂ mixed oxide crystals was investigated using molecular dynamics simulation. The primary model systems were isolated nanocrystals consisting of 5460 and 15,960 particles, featuring a free surface. The oxygen diffusion coefficient in the model systems increased with decreasing thorium content, which correlates with the decrease in the melting temperature of (Pu_xTh_1−x)O₂ crystals as x varies from 0 to 1. For plutonium contents exceeding 50%, the oxygen diffusion coefficients were virtually identical to those obtained for pure PuO₂.

For all thorium-to-plutonium ratios, the temperature dependence of D(T) exhibited non-linearity in the Arrhenius coordinates lnD = f(1/kT). The three linear segments of these dependencies corresponded to the superionic state, a transitional region, and the low-temperature crystalline phase. This shape of the calculated lnD = f(1/kT) curves is similar to the experimental temperature dependencies of ionic conductivity for structural analogs – CaF₂, SrF₂, BaF₂, and SrCl₂ crystals – reported in [8].

In the superionic state, the effective activation energy for diffusion was independent of the (Pu_xTh_1−x)O₂ composition. The obtained values, close to 1.0 eV, are in good agreement with existing computational results for uranium dioxide [14,16].

The transitional region between the superionic and low-temperature phases corresponded to diffusion coefficients exceeding 5·10⁻⁸ cm²/s. In this region, the results obtained from simulations using free-surface and periodic boundary conditions coincided. The effective formation energies of interstitial ions obtained in this region – E_f,eff(PuO₂, I) = 2.52 eV and E_f,eff(ThO₂, I) = 2.9 eV – are close to the experimental value of 2.3 eV measured for UO₂ at high temperatures (2000 K to 3000 K) in [6]. The effective activation energies for oxygen diffusion E_D in this region are elevated compared to low-temperature experimental data, which may be attributed to a difference in the mechanisms of interstitial anion generation.

For the first time, we investigated the low-temperature region in modeled PuO₂ and ThO₂ nanocrystals with free surfaces, where the formation energy for bulk anion interstitials is significantly reduced. Due to this reduction in energy, the lnD = f(1/kT) plot showed quantitative agreement with experimental data. For ThO₂, the calculated E_D value converged significantly with experimental results.

The existence of a similar temperature regime in real PuO₂ and ThO₂ crystals is supported by the agreement between the modeled lnD = f(1/kT) curves and the experimental ionic conductivity data of their structural analogs, CaF₂, SrF₂, BaF₂, and SrCl₂ [8].