A Mechanistic Study of Clustering and Diffusion of Molybdenum and Rhenium Atoms in Liquid Sodium

: Liquid Na is widely used as the heat transfer medium in high-temperature heat pipes based on Mo-Re alloys. In this study, ab initio molecular dynamics are employed in order to understand the interactions between the Na solvent and Mo or Re solute in the liquid phase. Both the temperature and concentration effects on the clustering and diffusion behaviors of solute atoms are investigated. It is found that Mo 2 and Re 2 dimers can be stabilized in liquid Na, and the higher temperature leads to a stronger binding force. Pure Re and Mo-Re mixed solutes can form tetramers at the highest concentration. However, for the pure Mo solute, Mo 4 is not observed. The diffusivities of a single solute atom and clusters are calculated. It is found that the Mo species diffuse faster than the Re species, and the diffusivity decreases as the cluster size increases.


Introduction
In recent years, photovoltaic cells, electrochemical energy storage devices and wind turbines have been greatly improved to reduce the energy risk and environmental problems caused by the utilization of fossil fuels [1][2][3][4]. However, these sustainable energy technologies have not met the demand of the fast developed modern society. Going beyond the technologies discussed above, nuclear fission has been delivering green and reliable energy for half a century, and nuclear energy is a competitive candidate to mitigate energy risks [5][6][7]. In addition, nuclear energy is also expected to power spacecraft [8,9].
Heat pipes (HPs) are the key device in the nuclear energy system for achieving a high efficiency and safety. In nuclear reactors, the high-temperature HPs usually use liquid sodium (Na) as the heat transfer medium, because liquid Na possesses a high latent heat of vaporization, low saturated vapor pressure and outstanding heat conductivity [10]. Molybdenum (Mo) alloys with a low rhenium (Re) content can be used as structural materials in high-temperature HPs for nuclear application due to its high melting point, good mechanical properties and adequate irradiation resistance [11,12]. The Mo-based HPs accompanied with liquid Na working fluid can be operated at the temperature window of approximately 1000 K to 1700 K [13]. Comparing with other refractory metals, Mo also exhibits a relatively higher corrosion resistance to liquid alkali metal. Inoue et al. [14] reported that Mo alloys showed a better corrosion resistance than Nb alloys in a liquid Na environment. Saito et al. [15] studied the corrosion of niobium (Nb)-based alloys and Mo-based alloys, and found that the weight of the corrosion product of Mo alloys was ten times smaller than that of the Nb alloys in liquid lithium (Li). In addition, their results also demonstrated that metal elements in Nb alloys are dissolved more easily.
The dissolution, migration and precipitation of alloy elements in HPs can change the properties of the material surface, which is harmful for the performance of HPs [16][17][18]. For Mo-based HPs using Na as the working fluid, it is still not well understood how Mo and Re atoms diffuse and accumulate in the Na solvent at the atomic scale. In the last few years, ab initio molecular dynamics (AIMD) have been successfully employed to investigate the properties of liquid alkali metal at extreme conditions [19][20][21]. In our present study, AIMD simulations are performed to reveal the interactions between liquid Na and Mo or Re solute atoms.

Computational Methods
All simulations in this study were performed on Vienna ab initio Simulation Package (VASP) [22,23] based on the density functional theory (DFT) [24,25]. The projector augmented-wave (PAW) method [26] was employed for describing the ion-electron interactions and the Perdew−Burke−Ernzerhof (PBE) functional [27] was used to describe the electron-electron exchange correlations. All AIMD simulations were carried out using NVT ensemble with the 400 eV energy cutoff of plane wave basis sets. AIMD simulations were run for at least 60 ps with a timestep of 2 fs. Only the Γ point was sampled in the first Brillouin zone. The liquid metal model was constructed by randomly distributing Na atoms in a 15 × 15 × 15 Å 3 box with a three-dimensional boundary condition. The number of Na atoms was determined by the liquid Na density reported by Argonne National Laboratory [28]. In this study, three temperature conditions of 700 K, 1100 K and 1600 K are considered, and corresponding number of Na atoms in the model are 75, 67 and 56, respectively, which correspond to the liquid Na density of 852 kg/m 3 , 756 kg/m 3 and 626 kg/m 3 , as reported in Ref. [28].

Results and Discussion
In this study, the pair correlation function g(r) is generated by the VASPKIT code [29] in order to characterize the structure of the liquid metal system. The function g(r) is formalized as where r is the radial distance, N is the total number of atoms in the modeling system and ρ is the average density of the system. The purpose of normalizing the g(r) function by the density is to ensure that the radial distribution approaches unity for the long radial distance. Following Equation (1), the partial pair correlation function between two elements A and B can be written as Figure 1 shows the total pair correlation function g(r) of liquid Na at 700 K, 1100 K and 1600 K after a 60 ps simulation. The primary peak of g(r) is located between 3 Å and 4 Å, which agrees well with Li et al's AIMD results [21]. It is interesting that our present simulation shows that there is a small peak before the primary peak at a relatively low temperature, and that this small peak disappears at 1600 K. This small peak was also experimentally observed in the g(r) function of liquid cesium when the temperature was below 400 K [30]. Figure 1 also demonstrates that the primary peak decreases and broadens obviously as the temperature increases. This phenomenon was also reported by Bickham and his collaborators [31]. The self-diffusion of Na atoms in the liquid phase is also evaluated. The diffusivity (D) can be calculated based on the Einstein-Stokes equation: here, < R 2 (t) > is the averaging mean-square displacement of Na. When recording the data, the system is equilibrated for 30 ps first, and then another 30 ps of simulation is run for collecting data. According to the present study, the Na diffusivities are 13.2 × 10 −5 cm 2 /s, 21.4 × 10 −5 cm 2 /s and 32.6 × 10 −5 cm 2 /s at 700 K, respectively. According to the result of Li et al, the Na self-diffusion coefficient at 723 K was 12.8 × 10 −5 cm 2 /s [21], which is close to our present study. Based on experimental data, Meyer fitted an equation to predict the Na self-diffusion coefficient as [32] It is worth noting that, in Meyer's work, the experimental data were collected below the temperature of 500 K. According to Equation (4), the Na self-diffusion coefficient is 19.2 cm −2 /s at 700 K, which is higher than the present and previous AIMD results. However, all computational results yield a magnitude order of 10 −5 cm 2 /s. The liquid Na system with Mo or Re atoms is modeled by replacing Na atoms with solute atoms. The effect of the solute concentration on the diffusion and clustering behaviors are investigated in the present study. Here, replacing one Na atom with a solute atom in the 15 × 15 × 15 Å 3 computational domain corresponds to a concentration of 0.5 mol/L, replacing two Na atoms with solute atoms corresponds to 1 mol/L and replacing four Na atoms with solute atoms corresponds to 2 mol/L. Figure 2 shows the partial pair correlation function of one solute atom in the computational domain (0.5 mol/L). It is found that the primary peak of the Mo-Na pair is lower than that of the Re-Na pair. Compared with the Re-Na pair, the primary peak position of the Mo-Na pair also shifts to the right. Therefore, it can be inferred that the solute Re atom can coordinate with more Na atoms than the solute Mo atom. The Bader charge analysis [33] is also performed in order to understand the interaction between the solute and the solvent from the aspect of the electronic structure. It is found that the net charge on the Mo atom is approximately −2.2~−3.0 |e|, and the net charge on the Re atom is approximately −2.8~−3.4 |e|. There is more negative charge transferred from the solvent atoms to the Re atom, which leads to the stronger Re-Na bonds. At the relatively higher concentration, the clustering of solute atoms is observed. Snapshots in Figure 3 show the final structures of the two solute atoms in liquid Na at different temperatures after 60 ps AIMD simulations. It is interesting to find that the solute atoms can stabilize as a dimer in the liquid Na. Figure 4 demonstrates the Mo-Mo distance evolution during the AIMD simulation. The dimer is dynamically stable once it forms, and the average Mo-Mo bond length is 1.79 Å, 1.85 Å and 1.91 Å at 700 K, 1100 K and 1600 K, respectively. Here the Mo-Mo bond length is calculated based on the last 1000 AIMD steps. At 700 K, the dimer forms after 38 ps, and the increase in temperature can speed up the formation of the dimers. It should be noted that Figure 4 cannot reflect the dimerization rate because each case is only run once.
The thermodynamic stability of the dimer is also estimated by the binding energy, which is expressed as Here, E Na N is the energy of the pure liquid Na without solute atoms, E M 1 Na N−1 is the liquid Na system with one solute atom, E M 2 Na N−2 is the energy of the system with two solute atoms and the subscript "N" represents the number of atoms in the model. All energetic terms are the average value of the last 1000 simulation steps. Table 1 shows the binding energy of dimers at different temperature conditions. The positive binding energy indicates that forming dimers is an exothermic reaction and that dimers are thermodynamically stable at the temperature range of 700 K to 1600 K.   It is worth noting that the binding energy increases as the temperature increases. The Bader charge analysis is performed to calculate the net charge on dimers shown in Figure 5. The net charge on Mo 2 is −3.6 |e| at 700 K, −2.4 |e| at 1100 K and −2.1 |e| at 1600 K. As discussed in Figure 2, the solute atom captures electrons from surrounding Na atoms, and fewer Na atoms coordinate with the solute atom when increasing the temperature. Therefore, the Mo 2 dimer at a higher temperature gains less of a net charge. It is worth noting that the binding energy of the dimer is proportional to the net charge ( Figure 5). For the Mo 2 molecule, the lowest unoccupied molecular orbital (LUMO) is empty doubly degenerate δ * 4d antibonding orbitals. The migration of electrons from the Na atoms to the Mo 2 dimer will occupy these antibonding orbitals. As the occupation of the antibonding orbitals increases, the Mo-Mo interaction weakens, which leads to a lowering of the binding energy. The net charge of the Re 2 dimer is −4.1 |e| at 700 K, −3.9 |e| at 1100 K and −2.3 |e| at 1600 K. As with the Mo 2 dimer, the binding energy of the Re 2 dimer is also proportional to the net charge. It is worth noting that Re 5d orbitals have one more electron than Mo 4d orbitals. When forming a neutral Re 2 dimer, its δ * 5d antibonding orbitals are halfly occupied, which indicates that the bond of neutral Mo 2 is weaker than Re 2 . In addition, the Re 2 dimer can capture more electrons from the Na solvent than the Mo 2 dimer. Hence, the occupation of the antibonding orbitals of Re 2 in the liquid Na is much higher than Mo 2 , which results in the much lower binding energy. Taking Re 2 in liquid Na at 700 K, for instance, the net charge on this diatomic molecule is −4.1 |e|. In this case, the antibonding δ * 5d orbitals are fully occupied, and the antibonding π * 5d orbitals are also halfly occupied, which lead to the lowest binding energy of 3.83 eV. The clustering behaviors of four solute atoms (corresponding to the concentration of 2 mol/L) in the liquid Na are investigated. The solute atoms can be four Mo atoms, four Re atoms or a mixture of two Mo atoms and two Re atoms. Figure 6 shows the total pair correlation function of the liquid system after a 60 ps simulation. For models including four Mo atoms, they still exhibit the typical liquid characterization, as shown in Figure 1. However, obvious peaks located around 2 Å can be found in Figure 6e-i. These peaks can be attributed to the segregation of solute atoms [21]. The atomic configurations of liquid systems after the 60 ps AIMD simulation are shown in Figure 7. It is interesting that the Mo solute and Re solute exhibit relatively different behaviors. Four Mo atoms form two Mo 2 dimers, but cannot form a stable Mo 4 tetramer. It should be noticed that these two dimers interact with each other and diffuse together (see videos in the Supplementary Material).
For the Re solute, Re 4 is observed at 1100 K and 1600 K after a 60 ps simulation. Li et al. also found a cerium tetramer (Ce 4 ) in the liquid Na using the AIMD simulation [21]. In their study, forming Ce 4 was attributed to the rejection of the liquid phase to the solute atoms. However, our simulation has demonstrated that there are attractive interactions between solvent Na atoms and Re atoms. It is known that the outmost shell of a Na atom is 3s orbital, whereas the outmost shell of a Re atom is 5d orbital. The energy level of the former is much lower than the latter. According to the linear combination of the atomic orbitals-molecular orbitals (LCAO-MO) theory, molecular orbitals are always formed by atomic orbitals with a small energy difference. Therefore, the Re-Re attractive interaction is stronger than that of the Re-Mo pair, which leads to the forming of the large Re cluster in the liquid Na. For the mixed solute, the formation of Mo 2 Re 2 is also observed, as shown in Figures 6g-i and 7g-i. The spin states of dimers and tetramers are also investigated. It is found that the Mo species does not show spin polarization. For the Re species, the magnetic momentum of Re2 is around 0.8 µ b per Re atom, whereas Re 4 does not show an obvious spin polarization.  In the present work, the accumulation behavior of six Mo or Re atoms is investigated at 1600 K. It is found that Mo atoms will form three interacting Mo 2 dimers, whereas Re atoms form an octahedral-like Re 6 cluster.
Understanding mass transport is essential for evaluating the performance and service life of alkali metal high-temperature heat pipes. In the present study, the diffusivities of single solute atoms and small clusters are calculated and shown in Figure 8. As expected, the diffusivity of any species is dependent on the temperature, and a higher temperature leads to a faster diffusion. In addition, the increase in cluster size can reduce the diffusivity. The diffusivities of the Mo x (x = 1, 2, 4) clusters are approximately two times as high as those of the Re x clusters, because the atomic weight of Mo (95.94 a. u.) is only half of the atomic weight of Re (186.2 a. u.). The diffusivity of a single Re atom is even lower than the Mo 2 cluster. As discussed above, there are not any stable Mo 4 clusters formed during the AIMD simulation, but two Mo 2 clusters can interact with each other and move together. The interaction between the two Mo 2 dimers can enhance the diffusion effects at low-temperature conditions, as shown in Figure 8. At 700 K, the diffusivity of one Mo 2 dimer is 2.53 × 10 −5 cm 2 /s, whereas the diffusivity of two interacting Mo 2 dimers is 2.81 × 10 −5 cm 2 /s. At 1100 K, the diffusivity of two interacting Mo 2 dimers is also approximately 0.3 × 10 −5 cm 2 /s higher than that of a Mo 2 dimer.

Conclusions
In this study, AIMD simulations are performed in order to understand the clustering and diffusion behaviors of Mo and Re atoms in liquid Na. For a single solute atom, the Re-Na interaction is stronger than the Mo-Na interaction. For models with two solute atoms, both Mo 2 dimers and Re 2 dimers are found in the Na solvent. The binding energy of the dimer is proportional to the occupation of the antibonding orbitals of the dimer. The higher occupation always leads to a lower binding energy. For models with four solute atoms, four Mo atoms form two interacting Mo 2 dimers, but a stable Mo 4 tetramer is not observed. However, for the pure Re solute or mixed solute, a stable tetramer is observed. The diffusivity of the Mo species is much faster than that of the Re species, and the diffusivity decreases as the cluster size increases.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/met11091430/s1, Video S1: Trajectory of four Mo atoms in Na at 1600 K; Video S2: Trajectory of four Re atoms in Na at 1600 K.