Diffusion of oxygen in hypostoichiometric uranium dioxide nanocrystals. A molecular dynamics simulation

A molecular dynamic simulation of diffusion of intrinsic oxygen anions in the bulk of hypostoichiometric UO2x nanocrystals with a free surface was carried out. The main diffusion mechanism turned out to be the migration of oxygen by the anionic vacancies. It is shown that in the range of values of the non-stoichiometry parameter 0.05 x  0.275 the oxygen diffusion coefficient D is weakly dependent on temperature, despite the uniform distribution of the vacancies over the model crystallite. The reliable D values calculated for the temperature T = 923 K are in the range from 310 to 710 cm/s, in quantitative agreement with the experimental data. The corresponding diffusion activation energy is in the range from 0.57 eV to 0.65 eV, depending on the interaction potentials used for the calculations. Keywords


Introduction
Recently, considerable efforts have been directed towards experimental research for the development of a nonaqueous technology for reprocessing irradiated nuclear fuel [1][2][3][4]. The key technology for highly efficient processing of such fuels is pyrolysis combined with molten salt electrolysis. This technology can be applied to highly burned fuel without long exposure time. The advantages of the electrochemical method are due to the high radiation resistance of inorganic molten salt, which makes it possible to dissolve highly radioactive nuclear fuel. The possibility of recycling the molten salt, and compactness of the equipment are also important.
This study is related to the problem of reduction of crystalline UO 2 to metallic uranium, which consists in the removal of oxygen from the crystal lattice. In the process of electrochemical reduction of uranium from UO 2 , oxygen is gradually released into the environment (for example, into the LiClLi 2 O molten salt surrounding the crystal surface). The crystal becomes hypostoichiometric, and the relative content of uranium increases in the course of uranium reduction up to the complete oxygen removal. The described technique has been well studied experimentally [1][2][3][4][5][6]; however, further optimization of the technology requires an understanding of the mechanisms of oxygen release at the atomic level. These mechanisms remain insufficiently studied to date.
The rate of oxygen release is limited by the coefficient of its diffusion in the hypostoichiometric crystal UO 2x . The experimental data on the diffusion of oxygen in UO 2x relate mainly to values of x within 0.1 [7][8][9]. The values of the oxygen diffusion coefficient obtained in [7] and [8] for the temperature used in the electrochemical reduction of uranium (923 K) differ by two orders of magnitude (being in the range from 210 6 to 110 8 cm 2 /s at x = 0.03). Higher values of x in UO 2x are achievable only with increasing temperature. The oxygen diffusion coefficient for x  0.2 has been measured at relatively high temperatures of 20002100 С [10].
To clarify the features of oxygen diffusion at significant deviations of UO 2x from stoichiometry, computational simulation of this process is of interest. To date, such studies have been performed mainly for hyperstoichiometric UO 2 that contains excess oxygen [11][12][13][14]. Oxygendepleted hypostoichiometric uranium dioxide has been studied much less [15][16][17].
Let us note that correct accounting for the ion charge transfer in a non-stoichiometric crystal requires the use of ab initio calculations. Nevertheless, the simulation of the crystals within the framework of classical molecular dynamics also remains relevant. An important advantage of the latter approach is the ability to study large systems over a long time evolution.
Molecular dynamic modeling of oxygen diffusion in UO 2 was previously carried out at temperatures significantly higher than 923 K [18][19][20], and in most of these studies, the periodic boundary conditions describing a quasi-infinite crystal were used. However, upon reduction in an electrochemical cell, the surface of uranium dioxide can become severely indented at the atomic level, which requires the presence of a free surface in the simulation. In this work, the presence of the surface was taken into account by modeling UO 2-x nanocrystals isolated in vacuum.
The results of classical modeling are fundamentally dependent on the choice of the particle interaction potentials. To simulate a UO 2 crystal, a number of sets of empirical potentials have been proposed, optimized for different problems (see, for example, reviews [20][21][22]). The oxygen diffusion coefficients calculated using different potentials could differ significantly. The applicability of the existing potentials under hypostoichiometric conditions has not been sufficiently studied. In the present work, three such sets of potentials [17,20,23] are compared.

The Model
The molecular dynamics model in this work was a UO 2x crystallite in the shape of an octahedron with a free surface. The specified stoichiometry was ensured by removing the required amount of oxygen anions, selected at random, from the crystal. In this case, the electroneutrality of the crystal was maintained by replacing part of the U q+ cations with U (q1)+ ions. Here, q is the effective charge of the uranium cation within a specific set of interaction potentials. At x = 0, the crystallite consisted of 2720 uranium ions and 5440 oxygen ions.
The x values for which the oxygen diffusion coefficients were calculated in this work varied from 0.003 to 0.275. This range is within the experimental existence of non-stoichiometric uranium dioxide (O/U from 1.65 to 2.25 [24]). Values of x above 0.1 in UO 2x are achieved only at temperatures above 1000 K [10][11]25]. Nevertheless, in the process of the electrochemical reduction of uranium at lower temperatures, the occurrence of local oxygendepleted regions at the boundary of the crystal with the electrolyte is not excluded. The temperature 923 K in this work has been chosen to match the common conditions of the electrolytic uranium reduction.
To ensure the symmetry of the charge distribution over the surface of the nanocrystal, neutral UO 2 molecules were used as "building" blocks. The two oxygen ions were placed on the same straight line with the uranium cation, using the displacements (0.25, 0.25, 0.25)a and (+0.25, +0.25, =0.25)a relative to the cation center, where a is the lattice constant. This made it possible to minimize the effect of surface charges on the bulk properties of the model crystals, which had to be taken into account considering the long-range character of the Coulomb force.
The interaction of the intrinsic ions was simulated by three sets of potentials [17,20,23], in order to compare the results of using each of the sets. Set I (MOX-07, S. Potashnikov et al., 2011 [20]) is characterized by a relatively accurate reproduction of both the mechanical and thermophysical characteristics of the stoichiometric UO 2 crystal and the energies of its own disordering [20][21][22]. Set II (Yakub-09, E. Yakub et al., 2009 [17]) provides the most accurate melting temperature of UO 2 . Set III (Busker-02, G. Busker, 2002 [23]), in contrast to potentials I and II, uses the whole charges of oxygen anions (2e) and uranium cations (+4e), which would allow, in a further study of the electrochemical cells, to simulate the transition of oxygen ions from the UO 2x crystal to the LiClLi 2 O molten salt without changing their charge.
Note that along with the potentials listed above, there are a number of other sets of potentials of similar quality, such as Basak-03 [26], Morelon-03 [11], Goel-08 [19]. Most of these sets are characterized by non-integral effective charges of the intrinsic ions. The lesser-known Busker potentials are studied in this work as an example of potentials with integral charges, which have acceptable quality at temperatures up to 1000 K, according to the review [20].
All the potentials were represented in general form that takes into account the repulsion of overlapping electron shells, the van der Waals attraction, and the possibility of chemical bond formation. The first term on the right-hand side corresponds to the long-range Coulomb interaction, the second term is the repulsive potential between the ion cores, the third term is the dispersive (van der Walls) attraction and the last term is the Morse potential accounting the covalent bonding energy attributed to U-O interaction. The quantities q i and q j are the effective ion charges, K E is the Coulomb's law constant, and the rest are fitting parameters that cannot be measured or computed separately. The parameter values are listed in Table  1. The non-Coulomb part of the potential UU was zero in the potential sets I and III. For the model nanocrystals in this work, the longrange Coulomb forces were calculated explicitly, due to the finite number of particles in the model. In order to avoid the non-physical effects of the long-range Coulomb interaction, the charge distribution over the nanocrystal surface was made as symmetric as possible.
The Newtonian equations of motion were integrated by the leapfrog method with a time step t = 310 15 s. The Berendsen thermostat was used, as well as an original procedure for correcting the crystallite rotation. Simulation times reached 200 ns.
To determine the oxygen diffusion coefficient in the bulk of the model crystals, we calculated the mean square of the displacement of anions that were separated from the surface by a distance no less than 1.75a, where a = 0.551 nm was the lattice constant of the UO 2.00 crystal at T = 923 K. There were 895 such anions at x = 0. The indicated lattice constant is the experimental value [27]. It coincides with the calculated values for all the three sets of potentials used in this work with an accuracy of 0.001 nm.
The simulations were carried out using an in-house developed software package, the same as in earlier works [2829]. The required high computing performance was obtained by parallelizing the calculations on CUDA architecture graphics processing units.  Fig. 2.

Results and Discussion
In this work, the oxygen diffusion coefficient was determined from the slope of the latter section using a generalization of the well-known Einstein relation in the form (2): Let us note that, at small deviations from stoichiometry x < 0.002, the equilibrium concentration of the vacancies in the bulk of the crystal turned out to be zero, so that in the equilibrium state of the crystallite diffusion of anions was not observed. An example of the corresponding graph <a 2 (t)> is shown in Fig. 3. The dependence <a 2 (t)>  is shown within 150 ns in order to show the region of its increase in more detail. Until the end of the simulation (200 ns), the mean square of the displacement remained constant. Fig. 4 illustrates the relationship between the equilibrium bulk concentration of the anionic vacancies (referred to the concentration of sites of the anionic sublattice) and the non-stoichiometric parameter x, which specified the composition of the model crystallites UO 2x as a whole, including the surface. These plots show that the distribution of the vacancies over the crystal was practically uniform, regardless of their concentration, for all the three sets of interaction potentials. Ideally, the concentration of vacancies and the parameter x would be related as [V] = 0.5x. Fig. 5 shows the results of calculating the oxygen diffusion coefficient at a temperature of T = 923 K depending on the stoichiometry of the UO 2x . It can be seen that the potentials I, II, and III give significantly different values of the diffusion coefficient, the maximum divergence of which reaches almost 500 times. The highest diffusion coefficients were obtained using potentials I, the lowest correspond to potentials III. Potentials II made it possible to obtain a better agreement between the calculation and experiment [8]. The calculation with potentials I also does not contradict the experimental data, taking into account the high values of the diffusion coefficient obtained in [7]. Potentials III clearly underestimate the oxygen diffusion coefficient. However, they reproduce the tendency observed in [7] towards a decrease in the diffusion coefficient with an increase in the deviation of the crystal composition from stoichiometry.
Regardless of the set of interaction potentials, the simulation predicts the absence of a strong dependence of the oxygen diffusion coefficient D on stoichiometry at x above 0.05. Experimental data [8] also indicate a weakening of the dependence of D on x with an increase in hypostoichiometry. Similarly, in the experiment [9], at values of x near 0.2 and temperatures of 20002100 С, the oxygen diffusion coefficient was practically independent on x. To explain these results, it can be assumed that simultaneous diffusion jumps of many anions in a confined space are hindered. Fig. 5 also shows the results of calculating the oxygen diffusion coefficient with periodic boundary conditions (6144 ions in the ideal supercell) for the potentials II. These diffusion coefficients are overestimated in comparison with the rest of the calculations of this work, although they are within the spread of the experimental data. More research is needed to clarify the reason for the discrepancy.

Conclusions
According to the comparison of the calculation results with experimental data, the pair potential approximation has provided a quantitative accuracy of modeling of the oxygen diffusion process in hypostoichiometric uranium dioxide nanocrystals, despite the simplicity of the model. Most reliable results were obtained using the sets of interaction potentials I and II, which are characterized by realistic effective charges of intrinsic uranium and oxygen ions. The diffusion coefficient values calculated using these potentials at a temperature of T = 923 K were in the range from 310 9 to 710 8 cm 2 /s. The calculated dependence of the diffusion coefficient on the stoichiometry of the crystal UO 2x at 0.05 x  0.275 is weak, in agreement with the experimental data [7][8][9]. A similar effect can take place in the process of electrochemical reduction of uranium from UO 2 dissolved in a molten salt such as LiClLi 2 O.