*2.1. Interatomic Potential*

Molecular dynamics (MD) simulation is an effective means of studying the microworld, providing insight into the atomic structure and the mechanism of bubble growth. An MD simulation program, LAMMPS [23], was employed for all simulations in this study. However, the accuracy of MD simulations depends largely on the potential function of the atoms used to support the simulation.

Yang et al. studied the influence of different potential functions on the growth of high-pressure Xe bubbles [3]. With the continuous addition of Xe atoms, the pressure of the bubbles is divided into two stages. In stage I, where the bubble pressure is monotonically increased, the bubble characteristics and the microstructure evolution of UO2 are relatively independent with the interatomic potential used; thus, the results of the five potential functions are approximately the same. In stage II, the bubble pressure is released and then fluctuates and the volume and pressure of the bubbles, as well as the evolution of the bubble configuration and the UO2 matrix as a function of the number of Xe atoms, highly depend on the potentials used [3]. In addition, the formation energy of Xe at the same location in UO2 varies significantly with the potential of UO2 and Xe–UO2 [3,24,25].

To ensure the accuracy of MD simulations herein, we evaluated different sets of potential functions and compared them with density functional theory (DFT) or DFT + U data to select a better comprehensive potential function for adding Xe to UO2 [26–30].

The first and second sets of potentials have the same UO2 potential and were developed by Basak et al. [31]. However, IPR [27] and Geng [32] potentials were used to interact with Xe–UO2. In the third and fourth sets of potentials, the potential reported by Morelon [33] was used for the UO2 matrix, whereas those reported by Chartier [29] and IPR [27] were used to describe the Xe–UO2 interaction. In the fifth and sixth sets of potentials, the CRG potential [25] was used as a potential function of UO2 and the interaction of Xe–UO2 was described by Cooper [34] and IPR [27]. Therefore, the six sets of potentials are called Basak/IPR, Basak/Geng, Morelon/IPR, Morelon/Chartier, CRG/Cooper and

CRG/IPR, respectively. Among them, three potential functions of UO2, are paired with different potential functions of Xe–UO2, because the interactions between Xe–U and Xe–O are determined in Xe–UO2. Thus, even in the case of the same UO2 potential function, we simulated it with different Xe–UO2 potential functions and the results are different.

We can determine the accuracy of using interatomic potentials by comparing them with the energetics of relatively small defects calculated using MD [3]. We developed several physically relevant benchmarks to compare the six sets of potentials on an equal basis. We placed Xe at different sites in UO2 (small, intermediate and large vacant sites). For small incorporation, we used the following defects: Xe in an interstitial (IntXe) and Xe on a uranium site with uranium in the neighboring interstitial (SubXe-IntU). The formation energies of the intermediate-sized incorporation sites are a Xe atom located at a U vacancy (Xe in Vu), an O vacancy (Xe in Vo), two nearest U and O vacancies (Xe in Vuo), or the three configurations of Schottky defect clusters (Xe in SD). Finally, we calculated the defect energetics of a single large incorporation site with Xe in a double Schottky defect cluster (Xe in 2SDs). To directly compare the MD values, a 2 × 2 × 2 supercell of the cubic fluorite unit cell (96 atoms for stoichiometric UO2) with a lattice constant of 5.454 Å was used in the molecular statics calculations. The simulation results are compared with the data in Table 1 [3]. The difference in energy is small, the maximum error value is 0.28 eV and the difference is 0.1 eV. The energy difference may be attributed to the difference in the machines used. However, in general, the results show that the potentials used are adequate. Then, through the MD simulations, we used a larger supercell (12,000 atoms) for better convergence to calculate the energies of Xe atoms at different sites in the UO2 lattice. In addition to the formation energies, we calculated the migration energy of interstitial Xe (E<sup>m</sup> IntXe ). The results are compared with selected DFT and DFT + U data in Table 1. The results obtained using the CRG/IPR potentials show the best consistency with the DFT results; thus, we chose CRG/IPR as the potential function for adding Xe to UO2 for our subsequent simulations.

**Table 1.** Comparison of energy of xenon (Xe) point defects in uranium dioxide (UO2) bulk calculated using various interatomic potentials at 0 K. Charge corrections for charge defects were not considered in the calculations herein.

