Atomic Insights into the Structural Properties and Displacement Cascades in Ytterbium Titanate Pyrochlore (Yb2Ti2O7) and High Entropy Pyrochlores

1. Introduction

Despite the environment-friendly and high-efficacy energy source nuclear energy has few distinctive concerns [1,2,3]. Commercial use of nuclear power is only evident when radioactive waste management is done under a proper strategy.[4]. Therefore, the reprocessing of spent fuel is an important part of ensuring nuclear safety and environmental safety. Nuclear wastes are in different forms depending on the sources and radioactive concentration [5]. Solidification of liquid waste is a common process for managing nuclear waste disposal. In this regard, the most common processes are vitrification and synroc methods [6]. Although borosilicate glass is being used frequently as nuclear waste substrate matrix but having low solubility of actinide elements restrict their usage [7]. An alternate substitute to solve this limitation are ceramics pyrochlore due to their superior durability, better potential at high temperatures and humid environments [8,9]. Moreover, ceramic has high values of waste loading as compared to glass [10]. The ceramic-based substrate with the minor addition of ionic concentration has higher radiation stability and excellent chemical and physical properties [11]. Long term radiation damage is critical for nuclear waste disposal strategy [12]. Pyrochlore has been developed rapidly as high-entropy ceramics waste substrate materials in recent years [13,14]. Future state of the art nuclear reactors involving the recycling of nuclear fuels and burning of minor actinide series and decontamination of fission fragments [15]. Titanate pyrochlore (${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$) has been demonstrated as an effective waste substrate material for pyrochemical reactions [16,17].

Minervini et al.; reported disorder in different pyrochlore oxides using the atomic scale simulation method [18]. Additionally, pyrochlore demonstrated stability range dependence on the relative size of cations [19]. Brendan et al.; have investigated the relationship between the structural and bonding energy in lanthanide pyrochlore oxides (${\mathrm{S}\mathrm{n}}_2{\mathrm{O}}_7$) and found that the position parameter of the oxygen vacancies is inversely proportional to the lattice parameter. Recently, Zhang et al, have reported the machine learning (ML) methods to determine lattice constants of different multi-substitutional pyrochlores in the range of 9-11 Å [20]. Moreover, the pyrochlore stannate (Ce₂Sn₂O₇) demonstrated temperature-dependent anisotropic nature [21]. Superconductivity of ternary pyrochlore oxide (${\mathrm{Cd}}_2{\mathrm{Re}}_2{\mathrm{O}}_7)$ reported by Sakai et al., at 1.1K , and the lattice constant was found to be 10.22 Å at room temperature [22]. Liyuan et al.; performed molecular dynamics to study structural and elastic properties of different pyrochlores under numerous combinations of A and B cations. It was found that displacement energy in pyrochlores strongly depends on the energy of the incident PKA’s as well as their atomic masses [23]. They also reported that lattice parameters and atomic radius predominantly affect structural as well as thermal properties. Chartier et al., performed molecular dynamics simulations to study displacement cascades in lanthanum zirconate pyrochlore (Ln₂Zr₂O₇) with uranium ions bombarded at 350 K with 6 keV along different orientations . It was observed model doesn’t lose its crystallinity [24]. It is reported in earlier studies that displacement cascade is a determinant of material radiation stability [25]. Atomic scale radiation damage studies in pyrochlores (Gd₂Zr₂O₇) have demonstrated healing mechanism during radiation [12,26]. Moreover, a combination of different anion and cations result in a radiation-resistant response of pyrochlores [27]. Pyrochlores have disordered structures due to irradiation and alteration of thermodynamics [28]. Pyrochlores are known for their radiation-resistant properties, primarily because of their inherent structural stability.. This structure allows for the migration and annihilation of defects, minimizing their impact on the material's overall properties.

High entropy pyrochlores refer to a class of materials that exhibit high configurational entropy due to the random distribution of multiple cations on the crystal lattice sites. These materials have been studied for various applications, including radiation tolerance. High entropy pyrochlore (HEPy) oxides synthesized with ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ are obtained by substituting the cations in pyrochlore structure and they have displayed higher radiation resistance than their individual pyrochlore [29,30,31]. Atomic scale simulation results have also displayed that radiation resistant are affected by addition of Zr content in HEPy [32,33].

Defects dynamics in pyrochlore type structures is quite difficult to determine in experimental studies due to complex structure and associate irradiation associated phenomena [23,34]. The displacement cascades involve the initiation of damage in the materials that alternately determine the long term possible outcomes. Atomic scale simulation methods can shed light on determining the fundamental material properties and interaction mechanisms in pyrochlores [35,36,37,38,39,40,41,42,43,44]. It is also important to mention that a limited number of experimental and simulation studies have been performed on the behavior of different pyrochlore under radiation [7,12,29,45,46,47,48,49,50,51,52,53,54]. Molecular dynamics simulations have been applied to study the influence of displacement cascades on the microstructural properties of different pyrochlores [7,12,20,48,49,54,55]. Computer simulations interpret the radiation damage as a predictive tool for processing the experiments. The simulation models describe the physics behind the phenomena responsible for radiation damage mechanism in irradiated material by providing valuable tools observed in nuclear power plants (NPP). Numerous radiation damage studies have been reported on metals and alloys [36,38,39,40,41,42,43,56,57,58,59,60,61,62,63,64,65,66,67,68]. This study will serve as a predictive model for providing insight into expected behavior in radioactive environments.

This article studies three aspects of pyrochlores oxides . Primarily, interatomic potential was established by parametrizing existing literature and later validated through calculation for equilibrium lattice, atomic radius threshold energy of different combination of A and B in pyrochlore structure and was performed. Furthermore, displacement damage cascades simulations of ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ were performed with each constituent elements with corresponding E_PKA~1, 2, 5, 10 keV in different incident directions. Moreover, high entropy alloys have been constituted with compositional adjustment of rare earth elements. The relation of irradiation stability was evaluated to determine the nature of ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ was compared with high entropy pyrochlores.

2. Computational Details and Methodology

2.1. Crystal Structure

Pyrochlore have cubic fluorite-type structures containing 5 or more elements with oxygen-deficient vacancies having ${\mathrm{A}}_2{\mathrm{B}}_2{\mathrm{O}}_7$ type structure with A (rare-earth) and B (transition metal) being different cations [69,70]. The ‘A’ position coordinated with 3⁺ cation while ‘B’ position coordinated with 4⁺ cation atoms, e.g. 3⁺ cation (La³⁺, Nd³⁺, Gd³⁺, Sm³⁺, Y³⁺) and 4⁺ cation (Zr⁴⁺, Ti⁴⁺, Mo⁺⁴) [71]. Moreover, 1/8 of the oxygen atoms are vacant to balance the charge [5]. The different ion combinations are used to compute different structural properties. The crystal structure is developed by Visualization of Electronic and Structural Analysis (VESTA) [72]. The lattice parameter for pyrochlore oxide was set as 10.194 Å. The initial model comprised of 8 unit cells with 88 atoms per unit cell with Fd$\overline{3}$m group having composition Yb₁₆Ti₁₆O₅₆. Periodic boundary conditions were applied along axes. The initial structure was equilibrated using steep decent method (SD) since this method is considered a quick approach for an optimal efficient solution [73].

Figure 1 displays the computational model and unit cell of pyrochlore, where (a) MD snapshot of a relaxed computational model under study and (b) a unit cell of generic representation of pyrochlore [74]. The characteristic crystal structure of pyrochlores consists of a network of corner-sharing tetrahedra, which provides an open framework that can accommodate defects induced by radiation damage [75,76]. Molecular dynamics simulations were performed by a large-scale atomic/molecular massively parallel simulator (LAMMPS) developed by Plimpton [77] and structure analysis and visualization was done by OVITO developed by ppt [78].

2.2. Interatomic Potential

In MD simulation realistic interatomic potentials are very crucial for reliable results [79]. The interatomic potential was modified by coupling potential developed in earlier studies [18,80,81]. This potential was test for replicating equilibrium lattice constants and lattice constants reported in earlier studies. The potential function is divided into two parts and all pair interactions are described separately. The short-range -combination of ZBL (Ziegler-Biersack-Littmark) potential and Buckingham potential are adjoined with fitting function for smooth truncation. The long-range potential- coulomb potential, and the algorithm uses the particle-particle particle-mesh solver (PPPM) summation algorithm with an accuracy up to 10^-4 and with r_cutoff =10 Å [82]. The Buckingham potential result in negative infinity due to short range whereas ZBL potential is more realistic for charged modes at short ranges interactions [83]. The piecewise interatomic potential of ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ is expressed in Eq. (1), where ${\mathrm{r}}_{1,2}$ is active range of potentials. The ZBL potential and Buckingham potential is fitted with fourth order exponential function ($V_{\mathrm{Fit}. \mathrm{fun}}\left({\mathrm{r}}_{\mathrm{ij}}\right)=c_1+c_2{r}_{ij}+c_3{r_{ij}}^2+c_4{r_{ij}}^3$) for smooth interconnection, where $c_{1-4}$ are fitting constants.

$$V（r_{ij}）=\left\{\begin{array}{ll}{V}_{ZBL}\left({\mathrm{r}}_{\mathrm{ij}}\right)& r_{ij}<{ r}_1\\ V_{\mathit{Fit}. \mathit{fun}}\left(r_{ij}\right)& r_1<r_{ij}<{ r}_2\\ V_{Buck.}\left(r_{ij}\right)& r_2<r_{ij}<{ r}_{cut}\\ 0& r_{\mathit{cut}}<r_{ii}\end{array}\right.$$

(1)

$\mathrm{V}\mathrm{ZBL}=\frac{1}{4\mathsf{\pi }{\mathsf{\epsilon }}_0}\frac{{\mathrm{Z}}_1{\mathrm{Z}}_2{\mathrm{e}}^2}{{\mathrm{r}}_{ij}}\mathrm{G}(\mathsf{\alpha })$ where$\mathrm{G}\left(\mathsf{\alpha }\right)=0.18118{\mathrm{e}}^{-32\mathsf{\alpha }}+0.5099{\mathrm{e}}^{-0.9423\mathsf{\alpha }}+0.2802{\mathrm{e}}^{-0.4029\mathsf{\alpha }}+0.02817{\mathrm{e}}^{-0.2016\mathsf{\alpha }}$ and $\mathsf{\alpha }=\frac{{\mathrm{r}}_{ij}({\mathrm{Z}}_1^{0.23}+{\mathrm{Z}}_2^{0.23})}{0.8854a_0}$ . The notation in above expression such as ${\mathsf{\epsilon }}_0$ represents dielectric constant $a_0$ is the Bohr’s radius, e is the electronic charge and ${\mathrm{Z}}_1{,Z}_2$ are atomic number of elements. The atomic numbers of Yb₇₀ , O₈ and Ti₄₈ . The potential function of Buckingham $V_{Buck.}\left({\mathrm{r}}_{\mathrm{ij}}\right)$potential is represented as function of distance in Eq. 2, where ${\mathrm{A}}_{\mathrm{ij}}$, $\rho$ and ${\mathrm{C}}_{\mathrm{xy}}$ are element dependent potential parameters.

$${\mathrm{V}}_{\mathrm{Buck}.}(r_{\mathrm{ij}})=A_{\mathrm{ij}}{\mathrm{e}}^{-\frac{r_{ij}}{p}}-\frac{{}_{{\mathrm{C}}_{\mathrm{XY}}}}{r\frac{6}{ij}}$$

(2)

The parameters for Buckingham potential as function interatomic distance for Yb and O satisfying condition in Eq. 1. The values of the parameters are set for each computed pairs in Table 1 reported by earlier studies [18,23]. For Ti-O, the value of ${\mathrm{C}}_{\mathrm{xy}}$ is adjusted for smooth truncation. Before testing the potential further separately , interatomic potentials were coupled through splined function to ZBL repulsive and Buckingham functions . The resulting functions are plotted in Figure 2. The connected potentials for cation-anion/anion-anion are computed through Table 1 , while cation-cations are defined by pair style automatically computed by LAMMPS mentioned in the previous studies[84,85].

2.3. Displacement Cascades Simulations

Displacement cascade simulations were performed through molecular dynamics code LAMMPS [86]. A system having dimension 10×10×10 with a total of 88000 were constructed . Our calculation for lattice parameters resulted in 10.09 Å with periodic boundary set along all axes. Interatomic potential function was employed as mentioned in previous section . After completing the preliminary settings, thermal relaxation was performed initially using canonical (NVT) ensemble at 300K with 0.001ps for a simulation time of 30ps. After the thermal relaxation, the model was relaxed again using microcanonical ensemble (NVE) at a time step of 0.0001ps for another 10ps to ensure the system has attained sufficient equilibrium before proceeding displacement cascade simulation. For damage cascade simulations , the recoil energies were chosen as 1 keV, 2 keV, 5 keV and 10 keV. The direction of PKA’s were set at 45° along xy plane at 7° to avoid channeling effects. The selected PKA were Yb atoms and PKA’s angles are all included with the positive Z axis to prevent channeling effect.

3. Results and Discussion

3.1. Calculation for Equilibrium LATTIce Constant

The sites on A and B in structure have different ions and their combination affects the properties of the pyrochlores. An important parameter of the unit cell is the lattice constant and its equilibrium properties . Figure 3 displays total energy as function of lattice constant of ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$. The relation concave curve indicates that energy of the system will reach the minimum at a certain equilibrium lattice constant. The curve fitting is obtained by a polynomial function. The equilibrium lattice constant was found to be 10.05 Å which agrees with results reported earlier [18,74].

3.2. Threshold Displacement Energy (E_d)

The incident energy of primary knock-on atom (E_PKA) exceeds to overcome the constraints of the entire force field and leave a vacancy. This minimum energy (E_min ) required to displace an atom from its lattice site is known as threshold displacement energy (E_d) [48]. The theoretical calculations displayed that off-site threshold energy depends on crystallographic direction, temperature, orientation of incident PKA, energy, and type of incident particles.

We have computed threshold energy of ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ initially by selecting Yb, Ti with O atoms as incident PKA along [100], [110] and [111] axes. The Ed is displayed along the orientation in Figure 4. It can be seen that E_{d (Yb)}>> E_{d (Ti)} > E_{d (O)}. We have seen that E_d is anisotropic. Our finding are in agreement with other studies [92,93].

The calculation of the threshold energy (E_d) of high entropy pyrochlore $\mathrm{YbYTiZr}{\mathrm{O}}_7$ formed by doping Y element at A site and Zr element at B site. The atomic number ratio for simulated pyrochlore with substitution set as Yb:Y:Ti:Zr is 1:1:1:1. Later A-site is doped with three elements as Y, Eu, and Gd while B- site is doped with total of four elements in equal atomic proportion of seven constituent elements, i.e. Yb: Y: Gd: Eu: Ti: Zr as 1: 1: 1: 1 :1 : 2: 2 result in formation of ${\mathrm{Yb}}_{0.5}{\mathrm{Y}}_{0.5}{\mathrm{Eu}}_{0.5}{\mathrm{Gd}}_{0.5}\mathrm{TiZr}{\mathrm{O}}_7$. Table 4 compares the average E_d of two simulated models with incident PKA’s directions. It is established that ${\mathrm{E}{}_d(\mathrm{Yb}}_{0.5}{\mathrm{Y}}_{0.5}{\mathrm{Eu}}_{0.5}{\mathrm{Gd}}_{0.5}\mathrm{TiZr}{\mathrm{O}}_7$) > $\mathrm{E}{}_d$ (${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ ). This establishes that HEPy are less prone to defects [92].

3.3. Effect of Displacement Cascades on Energy and Direction of PKA

Irradiated materials are associated with displacement cascades. The displacement cascades are divided into different subcascade during bombardment of incident neutrons and ion implantation. The displacement cascades with different incident PKA’s are repress ted in Figure 5.

As the energy of the initial PKA increases, the total number of defects of various atoms increases due to an increase in mean free path between atoms . Higher E_PKA produces larger area cascades which have higher possibility of reaching E_d [48]. The total number of defects will have a characteristic of rapid growth and then annihilate with time. This situation is generally referred to as defect recombination and permanent defect after radiation. In Figure 6, it can be noted that during the initial stage of the simulation cascade collision initiates a rapid accumulation of defects occurred in the pyrochlore lattice. The number of defects has reached a peak at 0.3ps and afterwards continuously reducing the size of cascade. The defects are stabilized at 0.9ps leaving permanent defects.

During the cascade simulations , incident PKA was located at the lower left corner and moved along the diagonal. The effect of incident direction is investigated by Yb incident with 5keV along X-axis inclined at 45° at different at an interval of 5° with 7° channeling angle along Z axis. The relation between incident angle as a function of time is expressed in Figure 7. It can be seen that the least number of defects smaller inclination angles have lower number of defects. This is due to the fact that higher inclination angles are associated with lower threshold energies of incident PKA’s as well as successive overlapping cascades [94,95,96].

In order observe the effect of ions on inclination angle ,we compared the effect of defect evolution at lower and higher inclination angle by 5keV Yb atom. It is also known from Table 4 that Ti has higher values of E_d as compared to other constituents. The relation between number of defects as function of time are displayed in Figure 8 for different incident direction of PKA . It is clearly seen that Ti atoms and Yb atoms are almost the same in 55° (Figure 8 (a)) whereas Yb are double in Figure 8 (b). In both cases O has a higher number of defects than other constituents and higher angle has twice the number of defects .Therefore, it is more difficult to produce more defects in 55° inclination under the same energy of incident PKA.

3.4. Displacement Cascades in HEPy

This effect of displacement cascades on HEPy are studies by doping different elements on A and B sites of pyrochlore ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ structure. We have selected Zr at B sites while A site is changed forming three kinds of HEPy i.e., $\mathrm{Yb}\mathrm{Y}\mathrm{TiZr}{\mathrm{O}}_7$ , $\mathrm{Yb}\mathrm{Gd}\mathrm{TiZr}{\mathrm{O}}_7$ , ${\mathrm{Yb}}_{0.5}{\mathrm{Y}}_{0.5}{\mathrm{Eu}}_{0.5}{\mathrm{Gd}}_{0.5}\mathrm{TiZr}{\mathrm{O}}_7$. The displacement cascades formed by Yb with 5keV and defects morphology in HEPy are compared with ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$. The comparison is based on the analysis and defect morphology.

Figure 9 displays the displacement cascades at its peak and after defect annihilation with Yb incident PKA of 5keV. The defect forms of ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ and other three high-entropy pyrochlores are compared. The subscript for each displays the peak defects N_max and after annihilation. It is noted that the area of displacement cascades is larger for ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ whereas number of defects formed by displacement cascades in HEPy are higher than pyrochlore whereas ${\mathrm{Yb}}_{0.5}{\mathrm{Y}}_{0.5}{\mathrm{Eu}}_{0.5}{\mathrm{Gd}}_{0.5}\mathrm{TiZr}{\mathrm{O}}_7$ has lower number of surviving defects . For all the models the direction of PKA coincides with the defects cluster direction. We have reasonable to believe that high-entropy pyrochlore has a better radiation resistance compared with ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$. It is also seen from the displacement cascade snapshots that characteristic of limiting the circumferential development of defects and reducing the damage volume in HEPy.

Higher energies incident PKA are reported to have higher number of surviving defects [17,97,98]. The number of defects as function of time are compared for HEPy and pyrochlore are plotted in Figure 10. It can be seen from Figure 10 that the number of defects corresponding to YGdEu doped HEPy are lesser in all cases. The surviving defects in undoped pyrochlore are relatively larger than doped. The defect growth rate and defect disappearance rate of and the defect number annihilation is the faster in HEPy specially ${\mathrm{Yb}}_{0.5}{\mathrm{Y}}_{0.5}{\mathrm{Eu}}_{0.5}{\mathrm{Gd}}_{0.5}\mathrm{TiZr}{\mathrm{O}}_7$ . It is preliminarily due to the fact that HEPy provides a large number of available atomic configurations, which can accommodate point defects and compensate for local structural changes induced by radiation. This configurational flexibility enhances the material's ability to self-heal and recover from radiation damage, further contributing to its lower radiation damage radius. This also establishes that HEPy shows a promising characteristic for radiation tolerance,. Thus, MD results in this study can be used to optimize the radiation resistant nature of HEPy by testing other models with orientation and other constituent elements in the composition.

4. Conclusion

Molecular dynamics simulations were performed to study the effect of displacement cascades in ytterbium titanate oxide ${(\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7)$ and high entropy alloys. The calculations for equilibrium lattice constant pyrochlore oxides. The displacement cascades in ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_{7 }$ with Yb, Ti, O of 1 keV, 2keV, 5 keV, 10keV along [100], [110], [ 111] were used to determine the threshold energy (E_d) on the radiation damage response of pyrochlore. The effect of displacement cascades on the formation of high-entropy pyrochlore $\mathrm{YbYTiZr}{\mathrm{O}}_7$ , $\mathrm{YbGdTiZr}{\mathrm{O}}_7$ , ${\mathrm{Yb}}_{0.5}{\mathrm{Y}}_{0.5}{\mathrm{Eu}}_{0.5}{\mathrm{Gd}}_{0.5}\mathrm{TiZr}{\mathrm{O}}_7$ and ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ was compared to analyze the radiation stability. The displacement cascades of different incident PKA display that number of defects are proportional to the incident energy . Higher incident angles have a lower number of defects. The defects morphology of HEPy was compared with ${\mathrm{Yb}}_2{\mathrm{Ti}}_2{\mathrm{O}}_7$ .