{'fraction_non_alphanumeric': 0.025499593091599604, 'fraction_numerical': 0.015191246948186996, 'mean_word_length': 4.742211838006231, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular dynamics study on planar clustering of xenon in UO2

Interatomic potentials of uranium dioxide are investigated on their applicability to model structural stabilities beyond fluorite phase by comparing with ab-initio results. A high pressure cotunnite phase and loosely stacking virtual crystal are involved in order to get a primary confidence for large-scale deformation modelings of UO2 under irradiation damages. The behavior of Xe atoms in UO2 fuel is studied with molecular dynamics simulations. Besides the ground state bubble geometry, we find that a planar distribution is also (meta-) stable for xenon under thermodynamic perturbations. The Xe atoms with a planar configuration are in a liquid state at a typical reactor temperature of 1000 K, which presents a modulated layer-structure near the interface with solid UO2. A planar defects loop remains after these Xe atoms are released out.


Introduction
Uranium dioxide (UO 2 ) is a standard nuclear fuel in pressurized-water reactors. In operation in power plants or in the context of direct disposal of spent fuel, a clear understanding of its structural, thermodynamical and kinetical properties is very important. Therefore many studies have been undertaken about the behavior of this material under irradiation. The crystallographic, elastic and thermodynamic properties, static and high frequency dielectric constants, point defects formation energies, activation energies for lattice migration have been studied very well today and can be easily found in abundant literature [1,2,3,4,5,6]. A number of work on molecular dynamics (MD) simulations for UO 2 also has been published in recent years with classical interatomic potentials to get various thermodynamic and transport properties [7,8,9,10,11]. However, most of these studies focused on the stable fluorite phase, and bears a limitation that actual UO 2 might deviate from this structure locally due to irradiation damages. It is the case especially in cascade simulations [12,13]. This provokes a question about the applicability of classical potentials to model the behavior of UO 2 with lager-scale deformations with defects.
On the other hand, as a kind of major fission product, the behavior of xenon (Xe) atoms in uranium dioxide also has attracted considerable experimental and theoretical attentions. Concern has been focused particularly on its important role in fuel swelling and this has accordingly led to a desire to obtain a greater understanding of the basic processes governing its migration and trapping within the fuel. The behavior of single isolated atoms of Xe in stoichiometric UO 2 and its interaction with Schottky vacancies have been well studied [14,15,16,17]. In this paper we prefer to investigate a large-scale collective behavior of Xe atoms.
From a point of view of experiment, as a kind of nuclear fuel, uranium dioxide has extensive sub-grain formation (grain sub-division) under energetic fission fragment irradiations. The restructuring starts after accumulating certain amount of rare gas as fission products, such as xenon in the matrix. Observing experimental results, one can assume that the process of sub-grain formation is initiated by bubble or planar structure of nanometer scale through decomposition and reformation of U-O-Xe complex. Sub-grains develop from a repeated cycle of damage (ballistic excitation and local heating), recovery (relaxation) and quenching. To investigate this dynamical process, molecular dynamics (MD) study will be employed using interatomic potentials with long range Coulombic interactions. In this paper, we show that a planar defect on [111] plane in UO 2 fluorite structure seems (meta-)stable when a nanometer scale of the defect is created. It is reasonable because the relaxation to the ground state bubble geometry might be prevent by an energy barrier originated from the rearrangement of local atomic environment. The mechanism, which we find by MD simulations, could be useful for understanding the nuclear fuel restructuring process at high burnup irradiations.
In following sections, we compare the results of empirical shell model and pair potential model with first-principles calculations. The variation of relative energies and structural stabilities of different phases of UO 2 is studied, to examine the transferability of these potentials beyond the standard fluorite phase. In section 3, we will discuss the molecular dynamics simulations of a planar clustering of xenon on [111] plane of fluorite UO 2 in details, followed by a conclusion.
2 Transferability of interatomic potential

Method of calculation
Interatomic potentials are checked by calculating the energy variation of bulk materials with density functional theory. Plane-wave pseudopotential method implemented in the code VASP [18,19] is used. Pseudopotentials model the potential yielded by nuclei and core electrons whereas the valence electron wave functions are expanded in a plane-wave basis. Spin-polarized local density approximation (LSDA) with Hubbard U term correction of the energy functional (LSDA+U) [20,21] is employed. The parameters for Hubbard term are taken as U = 4.5 eV and J = 0.51 eV, which have been checked carefully by S. L. Dudarev et al. [24,22,23]. All calculations employ projector-augmented wave (PAW) pseudopotentials [25,26] with a cutoff of kinetic energy for plane waves of 400 eV. Integrations in reciprocal space are performed in the first Brillouin zone with 18 irreducible k-points for fluorite structure and at least 28 irreducible k-points for other phases generated with the Monkhorst-Pack scheme [27]. The convergence of this parameter set is well checked. The energy tolerance for the charge self-consistency convergence is set to 1 × 10 −5 eV for all calculations.
For the part of interatomic potential models, newly parameterized shell model by C. Meis [7] and pair potential model by C. B. Basak [8] are used. The former provides the most accurate static and thermodynamic description for fluorite UO 2 whose details are listed in the table 1 of Ref. [7], and the latter also gives a reasonable result whose parameters are given in table 1 of Ref. [8]. When xenon is presented, a short ranged pair potential model is adopted to describe interactions between Xe atom and UO 2 . A Lennard-Jones (LJ) model is employed for Xe-Xe and Xe-O pairs whose parameters are determined by fitting to first principles calculated total energies of FCC-Xe and XeO 3 in CuAu 3 phase, respectively [28]. The Xe-U pair follows Jackson's adoption [14] which is a Born-Meyer (BM) function All parameters of these potentials are listed in table 1.
It is necessary to point out that the Xe-U potential parameterized by Jackson is as a part of a shell model where he studied the Xe-UO 2 system [14]. Therefore the performance of this potential may have some ambiguity. We improve it by fitting to ab initio energies and list its parameters in table 1 as BM-fitted.
To calculate the energy with classical potentials, GULP [29] code is employed where Ewald summation is involved to treat the long-ranged Coulomb interactions among uranium and oxygen atoms. To ensure the comparability of calculated energies, all structures have their geometry fixed both for ab initio and classical potential calculations. The cutoff radius of short-ranged Van der Waals potentials are set to 10.0Å and fixed for all calculations.

Results and discussions
To check the transferability of interatomic potentials beyond fluorite phase, cohesive energy in cotunnite phase with a space group of P nma that stabilizes under high pressure is calculated. An experiment showed the pressureinduced phase transition from fluorite to cotunnite phase starts at ∼42 GPa [30], whereas a first-principles calculated setout pressure is about 38 GPa [31]. Local pressure due to irradiation effects is usually higher than this threshold value and the excess energy transferred from irradiation tracks to lattice can overcome the possible barrier easily. Therefore it is very possible that a portion of cotunnite phase of UO 2 appears in irradiation damaged fuel.
The cell volume dependence of UO 2 cohesive energy is shown in figure 1, where each unit cell contains 12 atoms (U 4 O 8 ). The cohesive energy of fluorite phase at a volume of 166.375Å 3 calculated by VASP is taken as the reference point. The energies predicted by interatomic potential models are shifted so as to have the same value as ab initio results. It is evident from figure 1 that in high-compression and stretching regime the pair potential (Basak) and cotunnite (P nma) phases calculated by ab initio (VASP) compared with those predicted by interatomic potential models (shell model [7] and pair potential [8]). Note the latter curves have been shifted to have the same reference energy as ab initio results.
model becomes worse although a reasonable variation of energy is predicted around the equilibrium volume neighborhood compared with ab initio results.
On the other hand, the shell model (Meis) provides a rather accurate energy variation for fluorite phase (F m3m) within the whole studied volume region. However a big deviation is observed for cotunnite phase and the predicted energy variation is much poor by comparison with pair potential model. In this sense the transferability of shell model is unacceptable to treat UO 2 beyond fluorite structure.
A series of intermediate structures along one possible transition path from F m3m to P nma phase are also investigated 1 . The variation of energy along this path is shown in figure 2. The energy of cotunnite phase has been set as the reference point and all curves are shifted to have the same value at this point. The shell model parameterized by Meis has some deviations, the largest deviation from first-principles result occurs at the middlemost phase where an energy difference of ∼6.8 eV for 12 atoms is obtained. The shell model has a large error in description structural phase transitions of UO 2 . Obviously this fault originates from the classical treatment of electric polarization in the semi-empirical atomic model and damaged its transferability drastically. It seems to be a congenital problem for shell model and we can expect a similar situation when apply it to other ionic materials. Comparison of the relative energy of UO 2 from P nma to F m3m phase along a designed path calculated by VASP with interatomic potential models. The cohesive energy of P nma phase is set as the reference point.
When Xe atoms congregate to a local region and form, for example, bubbles, the structure of UO 2 near the interface boundary will departure from standard fluorite phase greatly and become loosely stacked. To investigate the applicability of interatomic potential models to this kind of circumstance, we employ a loose virtual structure in which one uranium atom occupies at (0, 0, 0) and two oxygen atoms occupy (1/4, 1/4, 1/4) and (3/4, 3/4, 3/4) sites in a cubic cell. Comparing with fluorite phase, this virtual structure can be viewed as derived from it by removing three uranium atoms on the face centers as well as the bonded oxygen atoms. Based on above discussions, we can expect a fail- ure of the shell model on this phase. The variation of cohesive energy plotted in figure 3 proves that it is true. The shell model parameterized by Meis [7] predicted a much deeper binding energy by comparing with ab initio result. The increment of energy under compression is slow while it is too fast under stretching. Oppositely, the performance of pair potential model is acceptable at least, although a shallower binding energy, a smaller lattice parameter and bulk modulus are predicted.
When one Xe atom is inserted into the (1/2, 1/2, 1/2) site of the previous loose virtual cell, as figure 4 showing, the performance of shell model deteriorated further. Its huge deviation from ab initio results bereaves its applicability to this circumstance. Under high compression region, the pair potential model also becomes poor, whose interaction of Xe-U pair is taken after Jackson that labeled as BM in table 1. This potential was not well parameterized because it was for a shell model [14]. To improve its quality in pair potential model, we reparameterize it by fitting to the ab initio cohesive energy curve. The resulted energy variation (associating with pair potentials of UO 2 parameterized by Basak) is shown as the dashed line in figure 4. It is necessary to point out that this treatment is expedient to get qualitative understanding of xenon behavior in UO 2 .  3 Classical molecular dynamics simulation

Simulation method
To model the collective behavior of xenon in UO 2 , the classical molecular dynamics method implemented in DL POLY 2 [32] code is employed, as well as the Basak's potential set for UO 2 [8]. When Xe atoms are involved, potentials listed in table 1 is added to the Basak's set where the parameters of (BM-fitted) rather than (BM) are adopted for Xe-U pairs. Previous section shows this set of potentials gives an acceptable description of energy variation for UO 2 -Xe system around large scale deformations, which provides us a confidence that the following simulations are qualitatively correct leastwise.
It is well known that Xe does not dissolve in UO 2 , which means phase separation is inevitable in nuclear fuel when xenon is presented as a kind of fission product. Usually this will lead to bubbles of Xe atoms. In nuclear reactor, however, Xe atoms (or clusters) may have other meta-stable structures around grain boundaries. Since the behavior of bubbles has been studied extensively [16,15], in this paper we prefer to focus mainly on a possible planar distribution of xenon in UO 2 .
Previous studies had suggested that Xe atoms diffuses in trivacancies (that is, bound Schottky trios consisting of a uranium vacancy and two oxygen vacancies) [33,34,35,17]. On the other hand, first-principles calculation using pseudopotentials approach revealed that a vacancy trio with a geometry lining up along (111) direction is the most stable one [36]. Therefore we intermingled UO 2 matrix by replacing bonded UO 2 trios (the yellow balls in figure 5) with Xe atoms on [111] planes randomly when prepare the initial configuration for MD simulations.
[111] plane is chosen just because it has the largest layer interval distance which can mostly preserve the UO 2 matrix when Xe is doped.
It is quite clear that the conventional cell unit of fluorite phase showed in figure  5 is inconvenient when simulating layer-related behavior along (111) direction. A three times larger cell unit with (111) direction aligning with Z axis is used instead. A 6 × 6 × 6 supercell based on this enlarged unit cell is employed for following MD simulations, corresponding to 648 conventional unit cells of fluorite structure and more than 7000 atoms. To check the size effects, we also use some supercells of 12 × 12 × 6 and 6 × 6 × 12, having nearly 30000 and 14000 atoms, respectively. In all simulations, a NPT ensemble is employed with zero pressure and 1000 K. To remove the artificialness introduced when preparing initial geometries, all configurations are heated up to 3000 K and then cool down firstly. 5000 simulation steps are used to equilibrate system before gathering data. The timestep is set as 1 fs; the cut-off radius is 7.5Å for short-ranged Van der Waals interactions and 12Å for long-ranged Coulomb interaction in real-space part of Ewald summation, with a Ewald precision of 1.0 × 10 −5 .

Results and discussions
Although it is well known that bubble is the most stable geometry for Xe in UO 2 , our following simulations show that if the initial distribution of xenon is far from a bubble geometry and close to, for example, a planar distribution, then the bubble structure can not be recovered just by usual thermodynamic perturbations. It is perhaps due to a high energy barrier between these two    Figure 6 illustrates that a layer-like structure is stable for a range of density of Xe in UO 2 .  figure 7(B) demonstrates a perfect planar distribution of Xe in UO 2 , and it is surprising that no any tendency towards bubble geometry was observed. Even increment temperature to 3500 K, quite similar results are still obtained.
Since the attractive interaction between Xe atoms is very small, this kind of aggregation evidently drives by uranium and oxygen atoms near the interface which attempt to conserve the fluorite structure of UO 2 as perfect as possible. The Xe-UO 2 boundary actually forms a solid-liquid interface here. The radial distribution function (RDF) shown in figure 8 illustrates that Xe atoms are in a liquid state while UO 2 matrix is in a solid state (see the RDF of Xe-Xe pairs in the inset). It can be seen clearly from figure 8 that UO 2 is in a perfect solid state: the nearest neighbor (NN) of O-U and O-O pairs comprise the first and second peaks of the total RDF, the next nearest neighbor (NNN) O-O pair and NN pair of U-U, as well as Xe-Xe NN pair contribute to the third peak. The total RDF of all atoms also exhibits a characteristic of solid phase. But the inset graph gives an essential characteristic of liquid for xenon itself. Therefore the behavior of Xe in UO 2 is quite analogous to a system with hard sphere liquid in a solid.
It is well known that near solid-liquid interfaces, for example a hard-sphere and hard-wall systems or a crystal and its melt, the liquid phase will present some oscillation of the density along the distance from the interface [37,38]. This behavior results from the geometrical constraining effect of the solid surface on the atoms or molecules of the liquid, which causes them to order into quasi-discrete layers, similar to the oscillation of RDF of liquid. For liquid Xe having an interface with FCC aluminum, similar behavior was observed [39]. It is very analogous to what was shown in figure 6, except here parts of solid fragment of UO 2 are mingled with Xe liquid on the same layer.
It would be very interesting to study a case with a large amount of Xe atoms in UO 2 , where some layers are occupied by xenon exclusively. We made a sand-   figure 9, which predicts a layer-like (or planar) behavior of liquid Xe clearly. The layer structure fades away along distance from the interface, but it still maintains beyond the fourth layer. We also found that when UO 2 cannot preserve a good crystal structure, the planar structure of Xe becomes blurry. Different from previous  studies that mainly focused on geometry constrained arised from solid surface, our simulation suggests that the planar distribution of Xe atoms depends on internal state of solid phase directly (interaction among atoms) and indirectly (structure relaxation and surface restructure). The most interesting is that although different number density of xenon will leads to different interval between layers, the interval distance of liquid Xe always equals to that of UO 2 .
It is difficult to completely recover the damaged structure after Xe is released out at about 1000 K. This temperature is far below the melting point and cannot provide enough kinetic energy to overcome recrystallization barriers. It is quite different from the situation under irradiations, where excess energy from fission tracks can melt or vaporize UO 2 locally and transport enough energy to matrix to conquer the barriers. Figure 10 compares the layer structure along (111) direction of uranium before and after Xe interstitials are released out. In spite of in a solid state, uranium atoms within the defected region still rearrange the structure from 3 layers into 2 layers. However it cannot recover completely and many vacancies and interstitials still remained, as figure 11 shown. Far from the defects left by released Xe atoms, UO 2 keeps a perfect fluorite structure except few oxygen interstitials. The residual defects loop inherits the planar geometry of foregoing Xe interstitials, partly due to the mechanical stability of solid UO 2 . The most remnant defects are uranium and oxygen vacancies. It is reasonable because a plenty of UO 2 molecules have been repulsed out of the simulation cell from their positions when fission gas Xe is introduced.
To eliminate the open space left by them, melting and annealing the material is necessary.

Conclusion
The transferability of classical potentials of widely used shell model and pair potential model beyond fluorite phase was investigated by comparing their predicted energy variations along compression and stretching of various structures with ab initio results. Some virtual structures and intermediate phases along a structural transition were involved. Results showed that despite shell model gives a precise description of fluorite phase of UO 2 , its applicability to other phases is questionable. The performance of pair potential model suggested that the fault might originate from the classical treatment of charge polarization. It arises a problem for shell model on how to control the balance between precision and transferability.
Using pair potential model, we studied the planar clustering of Xe aggregating in UO 2 . Providing a quasi-layer distribution of xenon initially, our simulations showed it would be stable even after the melting of UO 2 . The Xe interstitials in UO 2 are in a liquid state and have a layer modulation of density perpendicular to [111] interface, which (1) has a tendency to connect all Xe atoms together, except the case with a low Xe density, no isolated Xe cluster was observed; (2) has almost the same layer distance as the solid part of UO 2 ; (3) maintains beyond four layers at least. After releasing these Xe out, a planar defects loop left in the UO 2 matrix. We believe it has some effects on the formation of rim structure of nuclear fuel materials. In addition, it is necessary to point out that in this paper we did not consider the planar defects perpendicular to other direction (say, (001) or (110)) that we prefer to leave for later studies, where a similar conclusion can be expected.