First-Principles Study of Oxygen in ω -Zr

: Zirconium alloys, which are widely used as cladding materials in nuclear reactors, are prone to react with oxygen (O). Furthermore, the ω -Zr in zirconium alloys can signiﬁcantly increase the strength and hardness of these alloys, but there is a lack of reports on the behavior of oxygen in ω -Zr in the current literature. To investigate their interactions, we have studied the behavior of O in ω -Zr using the ﬁrst-principles approach. In this work, we examined the effects of vacancy and alloying elements (Nb, Sn) on the behavior of O in ω -Zr. The results show that O with a formation energy of − 5.96 eV preferentially occupies an octahedral interstitial position in ω -Zr. A vacancy reduces the formation energy of O in a tetrahedral interstitial position in ω -Zr. Nb and Sn decrease the formation energy of O in the octahedral interstitial position by 6.16 eV and 5.08 eV. Vacancy effectively reduces the diffusion barrier of O around it, which facilitates the diffusion of O in ω -Zr. Nb and Sn preferentially occupy the 1b and 2d substitution sites in ω -Zr, respectively. Nb makes the diffusion barrier of O in ω -Zr lower and promotes the diffusion of O in ω -Zr. Moreover, Sn makes the diffusion of O around Sn difﬁcult. It was further found that O is less prone to form clusters in ω -Zr and tends to independently occupy interstitial positions in ω -Zr. In particular, a single vacancy would make the binding energy between O atoms to be further reduced.


Introduction
Nuclear energy is currently considered to be the most ideal alternative energy source, having the benefits of being efficient, and renewable. Zirconium (Zr) alloys can maintain good physical and chemical properties at high-temperature and high-pressure environments. At the same time, these alloys have good compatibility with nuclear fuel and are the most commonly used fuel cladding materials for water-cooled reactors [1,2].
At ambient conditions, Zr is stable in the hcp structure (the α phase). At high temperatures, the α phase transforms to the β phase of the bcc structure. Both α phase and β phase can transform into a simple hexagonal structure (the ω phase) under pressure [3]. The ω phase of Zr has caught the interest of researchers [4,5]. It was shown that the Debye temperature of the ω phase is lower and the magnetization rate is lower compared to the α phase [6]. The yield strength and hardness of Zr alloys are influenced by the ω phase. A single α phase in Zr usually has a relatively low yield strength, while the ω phase will substantially increases the strength of the materia under consideration l. Catledge et al. [7] found that the hardness of the ω phase increases by 80% compared to the α phase. The ω phase is commonly considered to be the phase that strengthens zirconium alloys. The ω phase was first discovered by Silcock in 1958 in the structure of high-temperature β-Titanium alloys after quenching [8]. Subsequently, numerous static high-pressure experiments confirmed that the ω phase is a stable structure at high pressure [4,9], and Zong et al. [10] further confirmed this conclusion through impact experiments.
Due to the high reactivity between Zr and oxygen [11,12], it is important to control the oxygen content during the preparation of Zr alloys for nuclear application. Previous studies have shown that oxygen has a large solubility in α-Zr, which can form an interstitial solid solution and plays a strengthening role. At high temperatures, the solubility of oxygen in β-Zr reaches 10.5 at%, while the maximum solubility of oxygen in α-Zr can reach 28.6-35 at% [13]. Yamanaka et al. [14] found that an increase in oxygen content can increase the solubility of hydrogen in Zr alloys at a certain pressure, but the hydrogen solubility starts to decrease at a ratio of 0.177 of oxygen to Zr. Therefore, oxygen is often considered an alloying element in Zr or Ti alloys, but it causes the fracture toughness to decrease [15]. Paton et al. [16] found that oxygen interferes with atomic displacements during the phase transformation of the alloys and, thus, is able to inhibit the generation of non-thermal ω phases in the alloys. Williams et al. [17] found that oxygen can inhibit the formation of the isothermal ω phase in Ti alloys during aging. In contrast, Niinom et al. [18] found that the oxygen content in the ω phase is four times higher than in the β phase after aging the TNTZ alloy twice, and the high concentration of oxygen in the ω phase may be a stabilizer of the ω phase. This is in contrast to the conventional effect of oxygen on ω phase formation.
On the atomic scale, there is a lack of research on the interaction between ω-Zr and oxygen. More in-depth knowledge and understanding of the mechanism on the microscopic scale is needed to reduce harmful effects. In addition, one way to develop high-performance Zr alloys is to add other kinds of alloying elements. Considering that the oxidation behavior of Zr alloys is influenced by the content and distribution of alloying elements, a comprehensive understanding of the role of alloying elements in the oxidation process, especially Nb and Sn, is required. This is because Nb and Sn are mainly used in commercial zirconium alloys for alloying, which greatly improves the oxidation and hydrogen resistance of Zr alloys. The addition of niobium gives Zr alloys better resistance to creep [19] and corrosion. Jiang et al. [20] demonstrated that Nb can improve the oxidation resistance of alloy systems by reducing the concentration of oxygen vacancies and intermetallics in the oxide. Michael et al. [21] showed that the space charge in the oxide layer of zirconium alloys affects the parabolic of the Zr-Nb oxidation kinetics. Cui et al. [22] combined transmission electron microscopy (TEM), scanning electron microscopy (SEM), and Raman spectroscopy to analyze the breakaway oxidation of Zr-1 Nb-1Sn-0.1Fe in high-temperature steam. The results indicate that the structural changes at the O/M interface play an important role in kinetic acceleration and transition. Sn is present in Zr alloys in the form of a substitutional solid solution, which acts as a solid solution strengthening. Xu et al. [23] showed that the effect of O on corrosion enhancement is somewhat reduced when the Sn content increases from 1.0% to 1.2%. Zhao et al. [24] showed that the antioxidant properties of Zr-Sn-Nb alloy in steam and oxygen environments are superior to those of Zr-4 alloy. Therefore, the interactions between these alloying elements (Nb, Sn), and O atoms were also investigated.

Computational Methods
In this study, we used the Vienna Ab-initio Simulation Package (VASP) [25,26] code to perform density functional theory (DFT) [27,28] calculations. The electron-ion interaction was described using the projector augmented wave (PAW) [29,30] method, and the exchange-correlation potential was modeled using the Perdew-Burke-Ernzerhof (PBE) [31] form of generalized gradient approximation (GGA). The PAW pseudopotentials were taken from the VASP library. Throughout the calculation, the supercell with periodic boundary conditions was used to calculate defects properties. Based on the equilibrium time and accuracy considerations, the cutoff energy of the plane wave basis was set to 450 eV, and the same parameters were applied as in a previous study of Zr alloys. Moreover, the first Brillouin zone was sampled using a k-point mesh generated by the Monkhorst-Pack (MP) [32] scheme. The k-point spacing was approximately 0.2 Å −1 . The selected calculation parameters were tested to ensure that the energy convergence was less than 1 meV per atom. To avoid oxygen atom behavior being limited by the cell size, a 3 × 3 × 4 supercell containing 108 atoms was used. Before calculating the properties, all structures were fully optimized, including the relaxation atom positions and supercell volumes, until the force between two atoms was less than 0.01 eV/Å to ensure the accuracy of the calculation. The diffusion barriers of oxygen were calculated using the climbing image nudged-elastic band (CI-NEB) [33,34] method. Table 1 shows our calculated lattice parameters of ω-Zr. The experimental values and other theoretical values are listed for comparison. Our calculated lattice parameters of ω-Zr are a = 3.238a and c/a = 1.599. The experimental values [35] are a = 5.039 Å and c/a = 0.658, while the calculated values of Greeff et al. [36] are a = 5.050 Å and c/a = 0.623. It can be seen that our calculated lattice parameters are in good agreement with the experimental data and other theoretical calculations. It is well known that oxygen occupies mainly the interstitial positions in Zr [37], so we mainly consider the case where oxygen (O) is located at interstitial positions in ω-Zr. As shown in Figure 1, there are five main high-symmetry interstitial positions in ω-Zr: 1b(0, 0, 1 2 ), 2c( 2 3 , 1 3 , 0), 2e(0, 0, z), 3f(0, 1 2 , 0), and 3g( 1 2 , 1 2 , 1 2 ). The formation energy of an O atom at the interstitial position is calculated as follows:   Here E (Zr n +O) is the total energy of the supercell with an O atom in ω-Zr, E (Zr n ) is the total energy of a perfect cell including n Zr atoms, and E O 2 is the energy of an isolated O 2 molecule. According to our calculations, we found that when an O atom is located at the 2e positions, it relaxes to the nearby tetrahedral (Tetra) interstitial position, while the 3f position is the octahedral (Octa) interstitial position. As shown in Table 2, the most stable occupancy for an O atom in ω-Zr is an Octa position, corresponding to a formation energy of −5.97 eV, and the second most stable occupancy is a Tetra position, corresponding to a formation energy of −4.62 eV. Our calculations are in agreement with the values obtained by You et al. [37]. We investigated the diffusion behavior of O atoms in ω-Zr using the CI-NEB method. We assumed two types of intra-plane and inter-plane diffusion. For intra-plane diffusion, we assume four diffusion paths: O 1 → O 2 , with an O atom jumping from an Octa position into a first-nearest neighbor (1 nn) Octa position; O 1 → O 3, with an O atom jumping from the Octa position into a 2 nn Octa position; T 1 → T 2 , with an O atom jumping from a Tetra position into a first-nearest neighbor (1 nn) Tetra position; T 1 → T 3, with an O atom jumping from a Tetra position into a 2 nn Tetra position. For inter-plane diffusion, we assumed two diffusion paths: O 1 → T 1 , with an O atom jumping from an Octa into a Tetra position; and T 1 → O 1, with an O atom jumping from a Tetra position into an Octa position. This is shown in Figure 2.   Table 3 shows the diffusion energies of the above paths. For inter-plane diffusion, the diffusion barrier of an O atom jumping from an Octa to a Tetra (O1→T1 path) position is 2.00 eV, while the diffusion barrier of jumping from a Tetra to an Octa position (T1→O1 path) is 0.66 eV. For intra-plane diffusion, wefound a diffusion barrier of only 0.014 for the T1→T2 path, which seems to indicate a very fast diffusion of the O atom involved. A similar situation was found in ZrO2 [38], where the diffusion barrier between two O-sharing bonds with an Zr atom is 0.06 eV. Moreover, for the O1→O2 path, the diffusion barrier of O is up to 2.39 eV. The diffusion potentials of the other two intra-plane diffusion paths, O1→O3 and T1→T3 are 2.40 eV and 0.014 eV, respectively, which are the same as the paths O1→O2 and T1→T2. The initial diffusion path is set by the interpolation points between the initial and final states. However, the CI-NEB calculation shows that an O atom cannot jump directly to the 2 nn site. The O atom always jumps to a 1 nn site first and then to a 2  Table 3 shows the diffusion energies of the above paths. For inter-plane diffusion, the diffusion barrier of an O atom jumping from an Octa to a Tetra (O 1 → T 1 path) position is 2.00 eV, while the diffusion barrier of jumping from a Tetra to an Octa position (T 1 → O 1 path) is 0.66 eV. For intra-plane diffusion, wefound a diffusion barrier of only 0.014 for the T 1 → T 2 path, which seems to indicate a very fast diffusion of the O atom involved. A similar situation was found in ZrO 2 [38], where the diffusion barrier between two Osharing bonds with an Zr atom is 0.06 eV. Moreover, for the O 1 → O 2 path, the diffusion barrier of O is up to 2.39 eV. The diffusion potentials of the other two intra-plane diffusion paths, O 1 → O 3 and T 1 → T 3 are 2.40 eV and 0.014 eV, respectively, which are the same as the paths O 1 → O 2 and T 1 → T 2 . The initial diffusion path is set by the interpolation points between the initial and final states. However, the CI-NEB calculation shows that an O atom cannot jump directly to the 2 nn site. The O atom always jumps to a 1 nn site first and then to a 2 nn site, as shown by the solid line path in Figure 2b.

Oxygen-Vacancy Interaction in ω-Zr
In nuclear reactors, the envelope material produces vacancy defects due to irradiation [40]. Thus, we investigated the interaction between a single vacancy and an O atom in ω-Zr, based on the removal of a Zr atom from the ω-Zr supercell to form a vacancy model. There are two different symmetries of Zr atoms in ω-Zr, 1a(0, 0, 0) and 2d( 1 3 , 2 3 , 1 2 ), respectively. The vacancy formation energy can be calculated by Equation (2): Here, E Zr n−1 is the total energy of the ω-Zr supercell containing a vacancy. The results show that the vacancy formation energy for the 2d position is 1.26 eV, while the vacancy formation energy for the 1a position is 2.89 eV, indicating that the 2d vacancy is easier to form than the 1a; therefore, we do not consider the case of the 1a vacancy in subsequent discussion. The formation energy of an O atom in the Octa position increases of 0.47 eV from −5.96 eV to −5.49 eV under the effect of vacancy, while the formation energy in the Tetra position decreases of 0.91 eV from −4.62 eV to −5.53 eV. This indicates that the O atom is more stable in the Tetra position around the vacancy compared to the Tetra position in the pure ω-Zr bulk. As shown in Figure 3b, we found that all Tetra positions with a first-nearest neighbor of the vacancy position relax to a new Tetra position (T vac ). Therefore, we mainly consider the following diffusion paths for O atoms in the ω-Zr-vacancy system: T vac → T 3 , O 1 → T vac , and O 1 → O 2 . The diffusions pathways are shown in Figure 4.
As shown in Table 4, compared to ω-Zr, the diffusion barriers of O atoms in the O 1 → O 2 and O 1 → T vac paths decrease to 1.50 eV and 0.50 eV, respectively, while the the value for T vac → T 3 path increases to 1.03 eV. This indicates that the vacancy reduces the diffusion barriers of the O atom near itand hinders the diffusion of the O atom the distant interstitial position from the vacancy. This indicates that the vacancy favorably reduces the diffusion barriers for the O atoms near a the vacancy and increases the diffusion difficulty of the O atoms to interstitial positions farther away from the vacancy. In Ni [41], the diffusion barrier of O is also significantly influenced by the vacancy due to the interaction between O and vacancy. from −5.96 eV to −5.49 eV under the effect of vacancy, while the formation energy in the Tetra position decreases of 0.91 eV from −4.62 eV to −5.53 eV. This indicates that the O atom is more stable in the Tetra position around the vacancy compared to the Tetra position in the pure ω-Zr bulk. As shown in Figure 3(b), we found that all Tetra positions with a first-nearest neighbor of the vacancy position relax to a new Tetra position (Tvac). Therefore, we mainly consider the following diffusion paths for O atoms in the ω-Zr-vacancy system: Tvac→T3, O1→Tvac, and O1→O2. The diffusions pathways are shown in Figure 4.  from −5.96 eV to −5.49 eV under the effect of vacancy, while the formation energy in the Tetra position decreases of 0.91 eV from −4.62 eV to −5.53 eV. This indicates that the O atom is more stable in the Tetra position around the vacancy compared to the Tetra position in the pure ω-Zr bulk. As shown in Figure 3(b), we found that all Tetra positions with a first-nearest neighbor of the vacancy position relax to a new Tetra position (Tvac). Therefore, we mainly consider the following diffusion paths for O atoms in the ω-Zr-vacancy system: Tvac→T3, O1→Tvac, and O1→O2. The diffusions pathways are shown in Figure 4.

Effect of Alloying Elements (Nb, Sn) on Oxygen in ω-Zr
Similarly to the case of the α-Zr [42], the alloying elements Nb and Sn tend to occupy substitutional positions in the ω-Zr more than interstitial positions. The formation energy of the alloying atoms at interstitial positions and substitutional positions are defined as follows: Here E (Zr n +X) is the total energy of the supercell with an alloying element (Nb, Sn) in interstitial positions, and E (X) is the reference energy for the alloying elements. E (Zr n−1 +X) is the total energy of ω-Zr with an alloying element in substitutional positions. According to the previous description, there are two different substitutional positions in ω-Zr, 1a and 2d, and for the interstitial position, we mainly considered Octa and Tetra positions.   [43]. It can be seen that the formation energy of Nb and Sn in ω-Zr and α-Zr is close. By comparing the formation energy of Nb and Sn in different positions in ω-Zr, it was found that the alloying elements Nb and Sn prefer substitutional positions to interstitial positions. As shown in Figure 5, when the alloying elements (Nb, Sn) occupy the interstitial positions, the bond length of Sn to an adjacent Zr atom is 2.69 Å, which produces a significant deformation, corresponding to the same degree of deformation of the Zr planes on both adjacent sides. A similar deformation is obtained for the interstitial Nb atom, with a bond length of 2.57 Å to the adjacent Zr atom. When the alloying element (Nb, Sn) occupies the substitutional position, Nb remains within the base plane of the host lattice but deviates from the center of the substitutional position, and Nb causes a contraction of the Zr lattice. This indicates the formation of a relatively strong Zr-Nb bond. In contrast, Sn is located at the center of the substitutional position and has almost no effect on the Zr lattice.    If Nb occupies the substitutional position, the O atom is most stable in the Octa position near the Nb atom, with a formation energy of −11.04 eV. For Sn, the O atom also tends to occupy the Octa position (−12.12 eV). Compared to pure ω-Zr, the alloying elements (Nb, Sn) make the formation energy of the O atom in the Tetra position reduce to −9.84 and −9.77 eV for Nb and Sn, respectively. The results are shown in Table 6. It is noteworthy that, as occurs in vacancy, O relaxes to a new Tetra (T Nb ) position when it is located in a traditional Tetra position near Nb. We further investigated the diffusion of O atoms in ω-Zr with alloying elements. For ω-Zr with Nb, we assumed three paths: For ω-Zr with Sn, we constructed O 1 → T 1 , O 1 → O 2 , and T 1 → T 2 diffusion paths.  Table 7 shows the diffusion barriers of O in ω-Zr containing the alloying elements Nb and Sn. For ω-Zr with Nb, the diffusion of O atoms is similar to that in ω-Zr with a single vacancy, where the Nb atom effectively reduces the diffusion barrier for the O atom to diffuse around it, but the effect is not as strong as that of the vacancy on the O atom. In particular, when an O atom is located in an Octa position, the diffusion barrier for the diffusion path O 1 → T 1 is 1.32 eV, and for O 1 → O 2 is 1.05 eV, which is reduced of 0.68 eV and 1.34 eV, respectively, compared to the diffusion barrier in pure ω-Zr. Our calculated diffusion potential of O in the ω-Zr -Nb system is close to that of O in pure Nb calculated by Chen et al. [44]. This may be due to the stronger Nb-O bond formed in both the ω-Zr-Nb system and pure Nb. For ω-Zr with Sn, compared to pure ω-Zr, the diffusion of an O atom near Sn becomes difficult, and the diffusion barrier of the O 1 → O 2 path increases from 2.40 eV to 3.36 eV. For O 1 → T 1 , it increases from 2.00 eV to 2.86 eV. Thus, the Sn atom hinders the diffusion of O atom. The experiments [23] on the effect of oxygen content on the corrosion resistance of Nb-Sn-Zr alloys showed that an increase in Sn content can slower down the oxidation of zirconium alloys, which is consistent with our findings on Sn. Our results can also provide a microscopic explanation for their conclusions.

O Atoms Clustering in ω-Zr
We have investigated the formation of small O x clusters in ω-Zr. In pure ω-Zr, the most stable occupancy for an O atom is Octa. Therefore, we placed one O atom in the Octa position and the others in different near-neighboring Octa positions in turn. The clusters' binding energy E b (x) can be obtained as follows: Here E (Zr n +1O) is the total energy of ω-Zr supercell with an O atom, and E (Zr n+xO ) is the total energy of ω-Zr supercell with x O atoms. Positive values represent attraction  Table 8 shows the results of our calculations. In pure ω-Zr, when the second O atom is located at 3 nn, the binding energy of the O-O pair is 0.01 eV, indicating a weak attraction between the two O atoms. When the number of O atoms increases to three, the binding energies of all three near-neighboring positions are negative, meaning that the three O atoms repel each other and are not conducive to clustering in this region. In ω-Zr with a vacancy, considering that an O atom has the lowest formation energy in the Tetra position, we placed the first O atom in the Tetra position and the other O atom in the nearest Tetra position (1 nn to 3 nn). From Table 8, we can see that the binding energies are negative at all near-neighboring positions, indicating that the O-O pairs in ω-Zr with a vacancy are repulsive to each other and this repulsion is much greater than in pure ω-Zr. This seems to indicate that vacancy prevents O-O pairs from forming near a single vacancy.
We further considered the effect of the alloying element on the formation of O clusters in ω-Zr. In ω-Zr with Nb and ω-Zr with Sn, the Octa position where an O atom is located near Nb and Sn is the most stable. Therefore, we placed one O atom in the Octa position near the alloying element and the rest of the O atoms in the Octa positions from 1 nn to 3 nn away. For Nb, the maximum binding energy is 0.04 eV when the O atom is located in the third nearest Octa position, while the binding energy becomes negative at all positions when the number of O atoms is three. This situation is similar to that in pure ω-Zr. For Sn, the binding energy is negative at all three positions and O-O pairs are not easily formed near Sn. In all systems, the binding energy increases with increasing distance from the near-neighboring position, suggesting that O atoms in ω-Zr prefers to be independent in interstitial positions.

Conclusions
These results demonstrate that O atoms are more likely to occupy Octa interstitial position in pure ω-Zr and ω-Zr containing alloying elements (Nb, Sn), and they are more inclined to occupy Tetra interstitial position in ω-Zr with vacancy. In ω-Zr, a single vacancy is more likely to form at position 2d, and alloying elements (Nb, Sn) tend to occupy substitutional positions. The presence of a single vacancy facilitates the reduction of the diffusion barrier of an O atom between Octa positions and increases the diffusion barrier of an O atom between Tetra positions. Nb plays a similar role, but the effect is not as strong as that of vacancy. Sn increases the inter-plane and intra-plane diffusion barrier of O atoms in ω-Zr, making the diffusion of O atoms around Sn difficult. Meanwhile, the calculations show that only O-O pairs exist in pure ω-Zr and ω-Zr-Nb with relatively weak binding energies, the O cluster is not easily formed in ω-Zr, and O atoms prefers to occupy interstitial positions in ω-Zr independently. The results of this study provide a basis for further microscopic studies of the behavior of O atoms in ω-Zr.
Funding: This research was funded by the Fund of Science and Technology on Reactor Fuel and Materials Laboratory (No. 6142A06200211).

Data Availability Statement:
The data that supports the findings of this study are available within the paper.