Quantum chemical study of reaction mechanism between plutonium and nitrogen

The study of the reaction between plutonium and nitrogen is helpful in further understanding the interaction between plutonium and air molecules. Currently, there is no research on the microscopic reaction mechanism of plutonium nitridation reactions. Therefore, the microscopic mechanism of the Pu with N2 gas phase reaction is explored in this study, based on density functional theory (DFT) using different basis functions. In this paper, the geometry of stationary points on the potential energy surface is optimized. In addition, the transition states are verified by frequency analysis and intrinsic reaction coordination (IRC). Finally, we obtained the reaction potential energy curve and micro reaction pathways. Analysis of the reaction mechanism shows that the reaction of Pu with N2 has two pathways. Pathway 1 (Pu + N2 → R1 → TS1 → PuN2) has a T-shaped transition state and pathway 2 (Pu + N2 → R2 → TS2 → PuN + N) has an L-shaped transition state. Both transition states have only one imaginary frequency. According to the comparison of the energy at each stagnation point along the two pathways, and the heat energy emitted by the two reaction paths, we found that pathway 1 is the main reaction pathway. The nature of Pu-N bonding evolution along the pathways was studied by atoms in molecules (AIM) and electron localization function (ELF) topological approaches. In order to analyze the role of the plutonium atom 5f orbital in the reaction, the variation in density state along the pathways was measured. Results show that the 5f orbital mainly contributes to the formation of Pu-N bonds, and the influence of temperature on the reaction rate is revealed by calculating the rate constants of the two reaction pathways.


Introduction
The gas phase study of actinides reveals reaction details which are difficult to obtain experimentally. References to this problem can be found in several papers [1][2][3][4][5][6][7]. Plutonium is a complex actinide element, widely applied in the field of energy. It locates at the boundary between light actinide elements containing delocalized 5f electrons and heavy actinide elements containing localized 5f electrons. Its electron structure [8] endows it with some complex physical and chemical properties [9,10]. When exposed to air, plutonium readily reacts with oxygen, nitrogen, water, and hydrogen to form a variety of compounds, which cause surface oxidation and corrosion [11][12][13][14]. These phenomena have significant impact on the usage and storage of plutonium and as a result have become the focus of research in recent years. Multiple studies [15][16][17][18][19][20] on the interaction between actinides and small gas molecules have provided a theoretical and experimental basis to resolve these issues. However, due to the high toxicity and radioactivity of plutonium, basic scientific plutonium research still faces many challenges, and it is therefore necessary to study the reaction mechanism of the corrosion process from a microscopic viewpoint. In recent years, quantum chemical calculations on the detailed reaction mechanism of gas actinides atoms and gas small molecules have been explored, including Pu + O 2 [18], U + (U 2+ ) + H 2 O [19], and Pu + H2O [20]. Multiple reaction pathways of these reactions have been obtained along with the geometric structures of each stationary point along the reaction pathway. Some experts have also conducted theoretical and experimental studies on the interaction between small gas molecules and the plutonium solid surface. Experimentally, the kinetics of plutonium reaction with oxygen, water, and moist air have been studied in literatures [21][22][23], and the mechanism of water-catalyzed oxidation has been proposed. In terms of theoretical calculation, some researchers [24] used a density functional method to simulate the interaction between C, N, and O atoms and the plutonium surface. However, no research currently exists on the microscopic reaction mechanism between plutonium atoms and nitrogen molecules.
The main aim of this paper was to study the microscopic reaction mechanism of plutonium atoms in a simulated gas phase reaction with nitrogen molecules. Two reaction pathways have been found, and the geometric structures of all equilibrium and transition states in the reaction pathways have been optimized and demonstrated. Electron localization function (ELF) [25][26][27] and atoms in molecules (AIM) [28][29][30] theoretical methods were used to analyze the bond evolution process. By calculating the density state of each stationary structure, the role of plutonium atom 5f orbital in the bonding process was compared and analyzed. Furthermore, by calculating the reaction rate constants of the two reaction pathways at different temperatures, the reaction rates of plutonium and nitrogen along the two pathways at different temperatures are discussed.

Theoretical calculation methods
All the theoretical calculations involved in this paper are based on density functional theory (DFT) using the Gaussian 16 program package [31]. We used the hybrid functional method B3LYP [32,33], which not only couples the Becke switching functional with gradient correction and the Lee, Yang, and Parr correlation functional with gradient correction, but also processes its local correlation functional with Vosko, Wilk, and Nusair (VWN) local spin density [34]. In previous calculations, the B3LYP method has demonstrated the advantages of handling strongly associated systems, especially when calculating the reaction between plutonium atoms and small molecules of common air gases. The calculated results through this function match the experimental data the best [35]. For the structural optimization of each equilibrium and transition state, nitrogen atoms adopt the polarization basis set 6-31 g (d,p), and plutonium adopts the Stuttgart small nucleus relativistic effective pseudopotential [36] (SDD), which replaces the 60 electrons in shells 1-4, leaving the n ≥ 5 shells as valence electrons. It selects the target data as the total electronic energy of hundreds of electronic atom states and their ionic states (rather than the orbital energy of each valence), and the average difference between the pseudopotential and the total electronic energy of these states is calculated by relativistic all-electron calculation kept as small as possible by constantly adjusting the pseudopotential parameters. To accurately calculate the single point energy and vibration frequency of the optimized structure, nitrogen atoms are calculated by def2-TZVP, a high angular momentum basis set, while plutonium atoms are still calculated via a SDD basis set.
By analyzing the energy and vibration frequency, we can confirm whether the equilibrium structure is located at the minimum point and whether the transition state structure has only one imaginary frequency. Meanwhile, the intrinsic reaction coordinate (IRC) is calculated to verify whether the transition state structure is connected to the two minimum points. Since the spin-orbit interaction has little effect on actinide molecular structure [37,38], the spin-orbit effect can be disregarded in all calculations.
AIM theoretical analysis of chemical bonds mainly relies on topological analysis to find the most representative bond critical point [39] (bcp) in the interatomic interaction region, whose properties can be used to investigate characteristics of the corresponding chemical bonds such as strength and nature [40][41][42][43]. It is believed that the greater the electron density ρ(r) at the bcp position, and the smaller the electron potential energy density V(r), the greater the homogenous bond strength. The corresponding bond can be considered covalent if the electron density Laplace function ∇ 2 ρ(r) is negative [43,45]. The value of |V(r)|/G(r) (H(r) = G(r) + V(r), where G(r) is the kinetic energy density of electrons) and H(r)/ρ(r) can also be used to determine the interaction type between the two corresponding atoms [46]. The ELF shaded surface projection map of each equilibrium and transition state shows the bond evolution process more intuitively, and the ELF analysis results do not affect the selection of basis sets. The density of state (DOS) figure is useful for analyzing the structural characteristics of chemical systems. By calculating the partial DOS (PDOS) of the 5f orbital of each equilibrium and transition state structure and comparing them with the total DOS (TDOS) of the structure and the overlap population DOS (OPDOS) of the interaction between plutonium and nitrogen atoms, the contribution of plutonium atom 5f orbitals to the bond evolution process can be attained. All the above calculations were performed by Multiwfn software [44].
Based on transition state theory (TST) [45], we calculated the reaction rate constants of the two reaction pathways at different temperatures and used them to investigate the reaction rates in the temperature range of 273.15 ~ 2273.15 K. It is worth emphasizing that the thermodynamic data needed for calculation was obtained via Shermo software [47].

Reaction mechanisms
There are two different reaction pathways in the Pu + N 2 reaction: pathway 1 (Pu + N 2 → R1 → TS1 → PuN2) and pathway 2 (Pu + N 2 → R2 → TS2 → PuN + N). The geometric structures and relative energies of all the equilibrium and transition states in the two pathways are shown in Fig. 1 and Fig. 2, respectively. As shown in Fig. 1, for pathway 1, the plutonium atom approaches the nitrogen molecule from the direction perpendicular to the N2-N3 bond, forming the initial reactant R1. In pathway 2, the plutonium atom approaches the nitrogen molecule from a direction parallel to the N2-N3 bond, forming the initial reactant R2, whose three atoms are approximately collinear. By analyzing the data in Table 1, it can be seen that the N2-N3 bond length in R1 is 0.13163 Å, which is longer than the N2-N3 bond in R2. Moreover, as shown in Fig. 2, the R1 energy is higher than that of R2, indicating that R2 is more stable than R1. This indicates that plutonium atoms and nitrogen molecules are more likely to be activated when they are close to each other perpendicular to the N2-N3 bond and separate to the two nitrogen atoms in the nitrogen molecule. Similarly, as shown in Fig. 1, Fig. 2, and Table 1, in pathway 1, when the system passes a 109.4 kJ/mol (E TS1 -E R1 ) barrier from R1 to product PuN 2 , the length of Pu1-N2 and Pu1-N3 bonds gradually becomes shorter, and the distance between the two nitrogen atoms gradually increases, and the N2-N3 bond in transition state TS1 has been completely broken.
In addition, for pathway 2, from R2 to product PuN + N, the system passes through a barrier as high as 365.7 kJ/ mol (E TS2 -E R2 ), and the length of the Pu1-N2 bond gets  The relative energies of all intermediates and transition states involved in the two reaction pathways composites progressively shorter, while the length of the N2-N3 bond gets longer and eventually breaks completely. TS1 and TS2 have only one imaginary frequency, which verifies the rationality of the structure. By analyzing TS1 and TS2 structural parameters and the vibration patterns corresponding to their imaginary frequency, it can be seen that the plutonium atom and the two nitrogen atoms in pathway 1 have the tendency to form bonds, and the two nitrogen atoms tend to gradually separate. The plutonium atom in pathway 2 tends to form bonds with the closest nitrogen atom, and the two nitrogen atoms also tend to gradually separate.
In general, the barrier of reaction pathway 2 is significantly higher than that of pathway 1. The energy of products in reaction pathway 2 is 341.8 kJ/mol, which is higher than its reactant, while the energy of product in reaction pathway 1 is 237.9 kJ/mol, lower than its reactants. The energy of the transition state in pathway 1 is much lower than in pathway 2, and these phenomena indicate that pathway 1 is more likely to occur than pathway 2. The two pathways both occur along the triplet spin surface. It is worth noting that the Pu + N 2 reaction mechanism is similar to that of Pu + O 2 [11] and Pu + H 2 O [13] where plutonium atoms interact with gas molecules to dissociate them, and after dissociation, the atoms form bonds with plutonium atoms.
For pathway 1, it can be seen from Table 2 that there is a (3, − 1) bcp between the plutonium atom and each nitrogen atom in R1, TS1, and PuN2 molecules. Additionally, the topological properties of the two bcps are identical. Figure 3 shows that the consistency of Pu1-N2 and Pu1-N3 interactions in the reaction process can also be obtained. These phenomena indicate that the interaction between the plutonium atom and the two nitrogen atoms is completely synchronous along reaction pathway 1. For the R1 structure, the ∇ 2 ρ(r) [48,49] at the bcps between the plutonium and nitrogen atoms is positive, the values of V(r) and H(r)/ρ(r) are negative, and the value of |V(r)|/G(r) is between 1 and 2. Combined with previous studies [48][49][50][51], we can conclude that these phenomena indicate that the interaction between plutonium and nitrogen atoms is a closed-shell interaction with covalent characteristics. Similarly, the interaction between the plutonium and nitrogen atoms in TS1 is also a closed-shell interaction with a stronger covalency. However, the interaction between plutonium and nitrogen atoms in the PuN 2 molecule is a thoroughly covalent interaction. By analyzing the bcp parameters between the two nitrogen atoms in R1, we can see that the values of ∇ 2 ρ(r), H(r), and H(r)/ρ(r) are all significantly less than 0 and the V(r)|/G(r) values are greater than 2. This indicates that the N1-N2 bond in R1 is a strong covalent bond. However, there was no bcp between the two nitrogen atoms in the TS1 and PuN 2 molecules, indicating that the N1-N2 bond was completely broken during the transition from R1 to TS1. This view is also supported by a marked decrease in electron localization between the two nitrogen atoms, as seen in the TS1 ELF projection map. Through a comprehensive analysis of the data in Table 2, we find that the bcps properties between the plutonium atom and the two nitrogen atoms changed significantly during the reaction. The value of ρ(r) increased gradually; the values of V(r), H(r), and H(r)/ρ(r) were less than 0 and gradually decreased; the value of ∇ 2 ρ(r) changed from positive to negative; and the value of |V(r)|/G(r) increased gradually and eventually exceeded 2. Based on these changes, we believe that the interaction between the plutonium and two nitrogen atoms gradually increases and changes from the initial closed shell interaction to a covalent interaction; that is, the two Pu-N bonds in the product PuN 2 molecule are covalent bonds. It can also be clearly seen from the ELF projection map that with the progress of the reaction, the degree of electron localization between the plutonium and two nitrogen atoms gradually increases, and the two nitrogen atoms in the PuN 2 molecule, the product, have separated on both sides of the plutonium atom, and there are two synaptic valence basins between the Pu1-N2 and Pu1-N3 bonds. By using the same method to analyze the data of reaction pathway 2 in Table 2, we can see that the bcp values of ∇ 2 ρ(r), H(r), and H(r)/ρ(r) between the two nitrogen atoms in R2 are negative and their absolute values are relatively large, while the value of |V(r)|/G(r) is positive and greater than 2. Thus, we can conclude that the N2-N3 bond in R2 is a strong covalent bond and stronger than the N2-N3 bond in R1, which also indicates that plutonium atoms are more likely to activate and break the N2-N3 bond when they approach nitrogen molecules from the direction perpendicular to the N2-N3 bond. Since the bcp values of ∇ 2 ρ(r), H(r), and H(r)/ρ(r) between atom Pu1 and atom N2 in R2 are positive, we can consider the interaction between these two atoms as a closed shell interaction. Furthermore, since the value of ρ(r) and the absolute bcp value of V(r) are relatively small, the interaction is weak. The absolute values of ρ(r) and V(r) at the bcp between the two nitrogen atoms in TS2 are very small, and thus, we can assume that the interaction between the two nitrogen atoms is extremely weak and the N2-N3 bond has been completely broken. However, the bcp values of H(r) and H(r)/ρ(r) between atoms Pu1 and N2 are negative, the value of |V(r)|/G(r) is approximately 2, and the absolute values of ρ(r) and V(r) are much larger than that of R2. These phenomena indicate that the Pu1-N2 interaction in TS2 is much stronger than in R2, and it has converted to a covalent interaction. Based on the same method, AIM analysis of the product structure of reaction pathway 2 shows that the covalent interaction between Pu1-N2 is stronger than in TS2, so that a covalent bond between Pu1 and N2 can be considered to be formed in the product. In general, as the reaction progresses, the interaction between N2 and N3 gradually weakens, while the interaction between Pu1 and N2 gradually strengthens, and a covalent bond is formed. As seen from the R2 ELF projection map in Fig. 3, there is no valence basin between the plutonium and nitrogen atoms, but there is a bisynaptic valence basin between the two nitrogen atoms. When TS2 forms, there is no valence basin between the two nitrogen atoms, suggesting that the N2-N3 bond has been broken and the interaction between Pu1 and N2 is greatly enhanced. As can be seen from the product ELF projection map, a valence basin is formed between Pu1 and N2, and the interaction between the two nitrogen atoms is further weakened. This is consistent with the AIM results.

Density of states
As shown in Fig. 4, the black curve represents DOS, the red curve represents the PDOS for the 5f orbital of the Pu1 atom, and the blue curve represents the PDOS of the two nitrogen atoms. The green OPDOS curve examines whether molecular orbitals in different energy regions are bonding or anti-bonding or have little or no effect on the Pu1 and nitrogen atoms. A study by Tian Lu [44] indicates that if the OPDOS curve of a certain region is significantly positive, the molecular orbital of this energy region has a promoting effect on the combination of the two segments, which will have a bonding effect. If it is significantly negative, the molecular orbital of this energy region is unfavorable to the combination and will have an anti-bonding effect. Furthermore, if the value is close to 0, the effect on the combination is relatively small. The dotted lines show the HOMO Fig. 4 The TDOS, PDOS, and OPDOS for all intermediates and transition states (highest occupied molecular orbital) positions, and the contours of the stationary points' HOMO are also shown in this figure.
By analyzing the density of states of R1, TS1, and the PuN2 molecule in pathway 1, it can be found that in R1, the orbitals with energy lower than − 0.28 a.u. are basically contributed by the orbitals of the two nitrogen atoms, while those with energy higher than − 0.28 a.u. are basically contributed by the 5f orbital of the plutonium atom and the orbitals of the two nitrogen atoms. At the HOMO position, the PDOS of the plutonium atom 5f orbital is very close to TDOS, indicating that the HOMO is mainly composed of the plutonium atom 5f orbital and the participation of the two nitrogen atoms is weak. Moreover, the OPDOS curve at the HOMO position is positive, denoting that the interaction between the plutonium 5f orbital at the R1 HOMO position and the two nitrogen atoms can promote the formation of a Pu-N bond. In TS1, the orbitals with energy lower than − 0.52 a.u. are basically contributed by the orbitals of the two nitrogen atoms, while those with energy higher than − 0.39 a.u. are mainly contributed by the plutonium atom 5f orbital and the orbitals of the two nitrogen atoms. The OPDOS curve is clearly positive in the range that contains the HOMO energy. Moreover, the state density of the plutonium atom 5f orbital and the orbitals of the two nitrogen atoms peaks, signifying that the interaction between the plutonium atom 5f orbital and the orbitals of the two nitrogen atoms also promotes the formation of the Pu-N bond, and the degree of promotion is much greater than in R1. In the product PuN2 molecule, the TDOS in the low energy region is also mainly contributed by the orbitals of two nitrogen atoms, while the TDOS in the high energy region is mainly contributed by the plutonium 5f orbital and the orbitals of the two nitrogen atoms. The degree of promotion is also much greater than in R1. In general, the plutonium atom 5f orbital plays a major role in the formation of the Pu-N bond. As the reaction progresses along reaction pathway 1, the tendency of the plutonium atom and two nitrogen atoms to form chemical increases, and the effect of the nitrogen atoms on the formation of chemical bonds also gradually increases.
In the following section, we will discuss the electronic density of states of R2, TS2, and the PuN + N structure in reaction pathway 2. It can be determined that TDOS in the low energy region is basically contributed by the orbitals of the two nitrogen atoms, while TDOS in the high energy region is contributed by both the plutonium 5f orbital and the orbitals of the two nitrogen atoms. In TS2, the OPDOS curve within an energy range of − 0.31 to − 0.19 a.u. is positive and has the highest peak value. Thus, the interaction between the plutonium atom 5f orbital and the two nitrogen atom orbitals in this energy range strongly promotes the formation of the Pu-N bond. In the product (PuN + N) structure, similarly, the OPDOS curve is positive in the energy range of − 0.33 to − 0.22 a.u., and the interaction between the plutonium 5f orbital and the two nitrogen atom orbitals is also a bonding interaction. Consistent with reaction pathway 1, the plutonium atom 5f orbital contributes more to the formation of the Pu-N bond than the orbitals of the two nitrogen atoms. Through comprehensive analysis of all the structures shown in Fig. 4, we found that the plutonium atom 5f electron played a significant role in the formation of the Pu-N bond, and the role of the two nitrogen atoms gradually increased. Therefore, it can be said that the formation of the Pu-N bond is mainly caused by the interaction between the plutonium atom 5f electron and the two nitrogen atoms. This phenomenon is similar to the formation of the Pu-O bond in the reaction between plutonium and oxygen [11].

Kinetic analysis
To calculate the reaction rate constant, based on the TST [45] method, and referring to relevant literature [52][53][54], we deduced the calculation formula to be as follows: where σ [55] is the degeneracy of the reaction path and the value in this paper is 1; κ is the tunnel effect correction coefficient, and the method of calculating it can refer to the relevant literature [56,57]; k B is the Boltzmann constant, and the value is 1.381E-23 JK −1 ; T is temperature, and its unit of measurement is (K); N A is the Avogadro constant, and the value is 6.022E + 23; h is the Planck constant, and its value is 6.626E − 34 Js; R is the ideal gas constant, and the value is 8.31447 Jmol −1 K −1 ; ΔG 0 represents the Gibbs free energy of the reaction path under gas phase standard configuration pressure (1 bar), which is obtained by Shermo software [47]. k (s −1 mol −1 ) is the reaction rate constant we wish to calculate.
As shown in Fig. 5 It can be seen from Fig. 5 that R1 → PuN2 progresses easily in the range of 273.15-2273.15 K, while R2 → PuN + N is relatively difficult. It is worth noting that, with the increase of temperature, the increasing reaction rate constant k gradually slows down, indicating that the reaction rate constant is significantly influenced by low temperature, while the degree of influence on the reaction rate constant gradually decreases with the increase of temperature.

Conclusions
In this paper, the gas phase reaction of plutonium and nitrogen was studied via DFT. The results show the microscopic reaction mechanism in detail. When the spin multiplicity of the system was set to triplet, two reaction pathways were found. In addition, the structure of each stationary point in the reaction pathways was optimized. Simultaneously, we obtained the energy of each stationary point and analyzed the energy changes along the reaction pathways. For the bond evolution process, we used two different topological analytical methods (ELF and AIM) to investigate the changes in electronic properties at bcps. Furthermore, in order to investigate the contribution of the interaction between the plutonium atom 5f orbital and the orbitals of the two nitrogen atoms in Pu-N bond formation, we analyzed the density of states of each stationary point. Additionally, in order to ascertain the reaction rate of each reaction pathway at different temperatures, we calculated the reaction rate constant (k) of each reaction pathway at different temperatures using the Shermo software [47] and transition state theory (TST) [45]. All results are shown below. 1) Plutonium atoms and nitrogen molecules approach each other in different ways to form the initial complexes (R1 and R2) with different structures, and two different reaction pathways are generated as the reaction progresses. pathway 1: Pu + N 2 → R1 → TS1 → PuN 2 . pathway 2: Pu + N 2 → R2 → TS2 → PuN + N.
2) The results of energy and kinetics analysis of each structure in the two reaction pathways show that reaction pathway 1 is more likely to occur than reaction pathway 2. At the same temperature, there is a significant difference in reaction rates between them.
3) The results of density state calculations show that the interaction between the plutonium atom 5f orbital and the orbitals of the two nitrogen atoms makes an important contribution to the formation of the Pu-N bond. The energy region in which the interaction promotes bonding is generally smaller than the energy at the HOMO, whereas in the energy region where the energy is larger than that at the HOMO, the interaction between them is generally adverse to Pu-N bond formation. 4) The results of AIM analysis provide a quantitative understanding of the bcps properties in each structure in the two reaction pathways and show more clearly the electronic property changes during the breaking of N-N bonds and formation of Pu-N bonds. From the ELF analysis, the bond evolution process can be comprehensively investigated. It is worth noting that both analytical methods reached the same conclusion; that is, the proximity of plutonium atoms to nitrogen molecules will activate and break its N-N bond, and the separated nitrogen atoms will then bond with the plutonium atoms. Furthermore, we can see that N-N bonds are easier to activate when the plutonium atom approaches the nitrogen molecule from the direction perpendicular to the N-N bond than from one end of the nitrogen molecule.