Black-phosphorus-structured nitrogen (BP–N) is an attractive high-energy-density material. However, high-pressure-synthesized BP-N will decompose at low pressure and cannot be quenched to ambient conditions. Finding a method to stabilize it at 0 GPa is of great significance for its practical applications. However, unlike cubic gauche, layered polymeric, and hexagonal layered polymeric nitrogen (cg-N, LP-N, and HLP-N), it is always a metastable phase at high pressures up to 260 GPa, and decomposes into chains at 23 GPa. Here, on the basis of first-principles simulations, we find that P-atom doping can effectively reduce the synthesis pressure of BP-N and maintain its stability at 0 GPa. Uniform distribution of P-atom dopants within BP-N layers helps maintain the structural stability of these layers at 0 GPa, while interlayer electrostatic interaction induced by N–P dipoles enhances dynamic stability by eliminating interlayer slipping. Furthermore, pressure is conducive to enhancing the stability of BP-N and its doped forms by suppressing N-chain dissociation. For a configuration with 12.5% doping concentration, a gravimetric energy density of 8.07 kJ/g can be realized, which is nearly twice that of trinitrotoluene.
I. INTRODUCTION
Metastable allotropes of molecular and nonmolecular nitrogen have been extensively investigated as potential high-energy-density materials (HEDMs) for future applications. This interest stems from the fact that the N≡N triple bond in N2 is one of the strongest chemical bonds found in molecules, with a bonding energy of 954 kJ/mol. Transforming this triple bond into a double bond (N=N, bonding energy of 418 kJ/mol) or a single bond (N–N, bonding energy of 160 kJ/mol) significantly decreases the bonding energy and increases the total energy.1 Conversely, transforming N–N or N=N into N≡N results in the release of a tremendous amount of energy. Consequently, forms of polymerized nitrogen containing N–N or N=N bonds has attracted significant interest as HEDMs.
Nonetheless, due to the strength of the triple bond and the fragility of the single bond, synthesizing polymerized nitrogen at ambient pressure is challenging.2 Direct compression of molecular nitrogen or nitrogen-containing compounds (e.g., azides) often leads to pressure-induced symmetry-breaking phases,3–11 thus expanding nitrogen’s structural diversity. Theoretically, numerous polymerized nitrogen structures, such as rings,12 layers,13–16 clusters, chains,17–19 cages,20 ions,21 and polynitrogen compounds,22–32 have been predicted. Scientists have been attempting to synthesize polymerized nitrogen and polynitrogen compounds since the first isolation of nitrogen gas in 1772 and the first synthesis of the azide anion in 1890. For instance, using chemical methods, the N5+ ion was first synthesized in 1999,2 and there have been several breakthroughs in the synthesis of N5− since 2016.33 As an alternative, the application of high pressure has proven to be an effective and popular method for synthesizing polymerized nitrogen. The first polymeric all-nitrogen material, cubic gauche nitrogen (cg-N), was synthesized at 110 GPa in 2004.4 Layered polymeric and hexagonal layered polymeric nitrogen (LP-N and HLP-N) were subsequently synthesized at pressures of 126 and 244 GPa, respectively.34,35 Remarkably, prior to 2020, nitrogen had never been observed in an isostructural crystal structure similar to black phosphorus (BP), which possesses unique puckered two-dimensional layers.36 However, in 2020, a significant breakthrough was achieved with the first successful synthesis of BP-N, at a pressure of 146 GPa.37–39 This groundbreaking achievement completed the set of BP structures for Group V elements.
However, high-pressure-synthesized BP-N decomposed when the pressure was decreased to 48 GPa.37 Similar decomposition phenomena at low temperature have been observed in other high-pressure-synthesized polymeric all-nitrogen forms, including cg-N,4 LP-N,34 and HLP-N,35 which decompose at 42, 52, and 66 GPa, respectively. The instability of polymeric all-nitrogen materials at low pressure greatly limits their potential applications as HEDMs. Strategies to enhance the stability of polymerized nitrogen include forming compounds with alkali or alkaline earth metals to not only discover new polymerized nitrogen types but also decrease the synthesis pressure.40–49 Doping is another effective method,50,51 since it can redistribute the electrons in the system, and the dopants can act as bonding agents. Moreover, employing nano-confinement effects and surface modification can lead to further improvements in stability.52–57 Our previous research demonstrated that the decomposition of cg-N at low pressure was caused by surface instability, and it has been found that surface modification can realize stable cg-N at 0 GPa.58,59
Unlike other experimentally synthesized polymeric all-nitrogen materials, BP-N exhibits an unstable configuration that decomposes into chains at low pressure according to structural relaxation calculations.37 To obtain BP-N at 0 GPa, first, we need to maintain the BP-N structure and prevent the formation of a chain-like structure at low pressure, realizing structural stability. Subsequently, it is crucial to ensure the absence of imaginary frequencies in the phonon dispersions, realizing dynamical stability. Phosphorus and nitrogen are both in Group V of the periodic table and have similar electronic structures and chemical properties. Interestingly, BP itself demonstrates remarkable stability under ambient conditions.60 Moreover, BP exhibits promising properties and holds significant potential for various applications.36 As a result, phosphorus doping introduces a promising avenue for achieving low-pressure, stability of BP-N, which shares phase similarities with BP.
In this work, we investigate the stability of P-atom-doped BP-N at high pressure with different doping concentrations, doping uniformity, and doping configurations, using a first-principles method. Our results shown that P-atom doping plays a crucial role in maintaining the stability of BP-N. Uniform distribution of dopants in the BP-N layer helps maintain its structural stability at 0 GPa, while interlayer electrostatic interactions due to N–P dipoles enhance dynamic stability by preventing interlayer slipping.
II. METHODS
First-principles calculations of the structure and electronic properties based on density functional theory with a plane-wave basis set were performed using the Vienna ab initio simulation package (VASP).61 The generalized gradient approximation (GGA)62,63 parameterized by the Perdew–Burke–Ernzerhof (PBE)64 exchange–correlation functional and the projected augmented wave (PAW)65 method with N-2s22p3 and P-3s23p3 as valence electrons were used for the calculations. The atomic positions, lattice parameters, and cell volume were fully optimized using the conjugate gradient (CG) algorithm, and the iteration relaxation of the atomic positions stopped when all forces were less than 0.001 eV/Å. The convergences of total energy and stress tensor were set to 1 × 10−6 eV and 0.01 GPa, respectively. The Monkhorst–Pack method was used to analyze the reciprocal space of all Brillouin zones. The k-grid spacing was set to 0.20 Å−1, ensuring a K-point density that converged the energy within 0.5 meV/atom. The phonon dispersions were calculated using the phonopy code based on the finite displacement method. We tested various van der Waals (vdW) corrections using the DFT-D2, DFT-D3, and optB86 methods, and found that DFT-D3 and optB86 gave similar results. Therefore, the DFT-D3 method was employed for all structures to account for interlayer interactions.66–68 To compare the enthalpies of the investigated polymerized nitrogen, supercells containing 32 atoms were used to ensure comparability among BP-N, cg-N, LP-N, and HLP-N. Other calculations for P-atom-doped BP-N were performed using a 1 × 1 × 2 supercell with 16 atoms. For 12.5% doping concentration, all possible configurations were investigated, while for the 25% doping concentration, we randomly selected 50 nonidentical configurations. For the calculation of the decomposition enthalpy and energy density of BP-N, α′-P3N569 and α-N227 were chosen as the ground-state phases. VAPKIT70 and MULTIWFN71 were also adopted in our work.
III. RESULTS AND DISCUSSION
To investigate the decomposition mechanism of BP-N at low pressure, we fully relaxed the BP-N structure, with pressure gradually decreasing from 100 GPa. In Fig. 1(a), the variations of the longest and average bond lengths of BP-N are shown as functions of pressure. Above 25 GPa, there is a gradual and smooth variation in bond length, and the BP-N structure is still maintained. Moreover, the average bond length is close to 1.50 Å at 25 GPa, which is greater than the bond length of cg-N at 0 GPa (1.40 Å).73,74 However, when the pressure falls below 25 GPa, the lengths of the bonds between nitrogen atoms increase dramatically, indicating that a phase transition occurs near this pressure. As shown in Figs. 1(b) and 1(c), the BP-N will decompose to nitrogen zigzag chains at a pressure of 23 GPa. According to the results presented in Fig. S1 (supplementary material), the occurrence of imaginary frequencies in zigzag chains at 23 GPa implies that they represent an intermediate stage in the decomposition of BP-N and will undergo further decomposition.
(a) Variation of longest (red circles) and average (black circles) bond lengths of BP-N with pressure. The relaxation proceeds from high to low pressure, leading to an increase in bond lengths as the pressure decreases. (b) and (c) BP-N structures at 25 and 23 GPa, respectively. Stability is maintained at 25 GPa, whereas decomposition into N-chains occurs at 23 GPa owing to breaking of the longest bonds. (d) and (e) Phonon spectra of BP-N at 40 and 30 GPa, respectively.
(a) Variation of longest (red circles) and average (black circles) bond lengths of BP-N with pressure. The relaxation proceeds from high to low pressure, leading to an increase in bond lengths as the pressure decreases. (b) and (c) BP-N structures at 25 and 23 GPa, respectively. Stability is maintained at 25 GPa, whereas decomposition into N-chains occurs at 23 GPa owing to breaking of the longest bonds. (d) and (e) Phonon spectra of BP-N at 40 and 30 GPa, respectively.
The phonon dispersions of the BP-N structure from 100 to 30 GPa were also calculated, and the results at 40 and 30 GPa are shown in Figs. 1(d) and 1(e), respectively. We find that beyond 40 GPa, BP-N exhibits dynamic stability, but its phonon dispersion at 30 GPa exhibits imaginary frequencies. The decomposition process of BP-N can be summarized as follows: with decreasing pressure, the BP-N will become dynamically unstable at 30–40 GPa and then decompose to nitrogen zigzag chains at 23 GPa.39,75
Since BP exhibits exceptional stability under ambient conditions, we seek to use P-atom doping to enhance the stability of BP-N, which shares phase similarities with BP. The effects of P-atom doping on BP-N stability were investigated, and the enthalpy–pressure curves at different P-atom doping concentrations for various phases, namely, cg-N, BP-N, LP-N, and HLP-N, are shown in Fig. 2. Among the pristine cases [Fig. 2(a)], cg-N has the lowest enthalpy at pressures up to 186 GPa, where a phase transition to LP-N occurs. Finally, HLP-N becomes the ground state at pressures higher than 280 GPa. The enthalpy of BP-N at high pressure is always higher than those of cg-N, LP-N, and HLP-N. BP-N is the metastable phase in the pressure ranging from 60 to 260 GPa, which is consistent with previous computational results.37–39
Variations of enthalpy with pressure at different P-atom doping concentrations for cg-N, BP-N, LP-N, and HLP-N: (a) without P-atom doping; (b) 3.125% doping; (c) 6.25% doping. The enthalpy values are referenced to those of cg-N. For 3.25% doping, there is one unique configuration for cg-N, one for BP-N, four for LP-N, and four for HLP-N (only stable configurations with minimal enthalpy were used). The same is found for 6.25% doping concentration.
Variations of enthalpy with pressure at different P-atom doping concentrations for cg-N, BP-N, LP-N, and HLP-N: (a) without P-atom doping; (b) 3.125% doping; (c) 6.25% doping. The enthalpy values are referenced to those of cg-N. For 3.25% doping, there is one unique configuration for cg-N, one for BP-N, four for LP-N, and four for HLP-N (only stable configurations with minimal enthalpy were used). The same is found for 6.25% doping concentration.
As shown in Fig. 2(b), with 3.125% P-atom doping, the phase-transition sequence does not change, and only the phase-transition pressure is modified; for example, the phase-transition pressure for cg-N to LP-N decreases from 186 to 140 GPa. Interestingly, after the P-atom doping concentration has increased to 6.25%, as shown in Fig. 2(c), the phase-transition sequence is changed significantly, and a sequence of cg-N to BP-N to HLP-N is observed, resulting in BP-N becoming the most stable phase at pressures ranging from 154 to 234 GPa. The results indicate that P-atom doping can enhance the thermodynamic stability of BP-N at high pressure. It should be noted that although BP-N exhibits a thermodynamically stable pressure region spanning 80 GPa at 6.125% P-atom doping, it still decomposes into N-chains at 0 GPa. In the following, we discuss the effects of P-atom doping on BP-N stability at higher doping concentrations.
We constructed configurations with higher doping concentrations (12.5% and 25%), as illustrated in Figs. S2 and S3 (supplementary material). Figure 3(a) illustrates the structural stability and uniformity of the various configurations at 0 and 50 GPa. We find that in the cases with higher doping concentrations, more than half of the P-atom doping configurations are able to stabilize the BP-N structure at 0 GPa, indicating that P-atom doping can prevent BP-N decomposition into N-chains. Furthermore, a distinct boundary separates stable and unstable structures, with uniformity serving as the key separator. Particularly noteworthy is the observation that the structures exhibit instability when their uniformity falls below 0.80, especially in the case of 25% doping concentration. The observed instability can be attributed to clustering of P atoms. Owing to the significant difference in radius between P and N atoms, a large stress will be introduced between P-atom clustering regions and N-atom clustering regions, which is not conducive to structural stability. More importantly, without P-atom doping to enhance stability, the regions with N-atom clustering will decompose. Moreover, pressure is also conducive to stability, as evidenced by the greater number of configurations that stabilize at 50 GPa. In addition, we conducted an analysis of the relationship between the energy and uniformity of configurations with 12.5% and 25% P-atom doping at 0 GPa. Our findings, depicted in Fig. 3(b), highlight that structures exhibiting higher uniformity tend to possess lower energy.
(a) Relationships between structural stability and uniformity at 0 and 50 GPa. The blue circles represent stable structures, and the orange circles depict unstable structures for structural relaxation. (b) Relationship between energy and uniformity at 0 GPa for 12.5% and 25% doping concentrations. The circles represent stable configurations, while the black lines are to guide the eye.
(a) Relationships between structural stability and uniformity at 0 and 50 GPa. The blue circles represent stable structures, and the orange circles depict unstable structures for structural relaxation. (b) Relationship between energy and uniformity at 0 GPa for 12.5% and 25% doping concentrations. The circles represent stable configurations, while the black lines are to guide the eye.
Since BP-N has a two-dimensional stacking structure, interlayer interactions may play a significant role in its dynamic stability, and so investigation of the positions of P atoms between the layers is necessary. A tendency toward uniform doping between layers is observed. In the case of 12.5% doping concentration, for all stable configurations, it was found that each layer contains an equal number of P atoms. Upon analyzing them, we found that their single-layer configurations were equivalent, with the only difference being in the relative positions between the interlayer atoms. We calculated the energy of these configurations and studied their relationship to the distance of P atoms between the layers, as shown in Fig. 4(c). The configuration labeled by the red triangle has the lowest energy, and the P atoms between the layers are arranged perpendicularly. Notably, this is the only configuration without imaginary frequencies. Its structure, along with phonon spectra, is shown in Fig. 4(a). The configuration indicated by a blue triangle in Fig. 4(c), which is one of those exhibiting imaginary frequencies in phonon dispersions, has the structure and phonon characteristics depicted in Fig. 4(b). Other structures with imaginary frequencies are similar to this one, with the imaginary frequencies occurring at the Z point. The presence of imaginary frequencies at the Z point indicates that interlayer sliding has occurred, as revealed by an analysis of vibrational modes [the black arrows in Fig. 4(b) indicate the direction of this vibration]. In the configuration with the lowest energy, the lowest frequency at the Z point is found a certain distance away from zero. This signifies that this particular doping method is beneficial in preventing interlayer sliding. The strengthening of interlayer interactions can be attributed to the contribution of electrostatic interactions.
Structures alongside their respective phonon dispersions for (a) the configuration possessing the lowest energy and (b) a representative dynamically unstable configuration with 12.5% doping concentration. The black arrows in (b) indicate the direction of vibration for the Z-point imaginary frequency. (c) Variation of energy as a function of interlayer P–P distances. Each circle represents a unique configuration, with the red triangle indicating the configuration with the lowest energy and the blue triangle one of the dynamically unstable configurations. (d) Schematic of interlayer electrostatic interaction for configurations with N–P dipoles. Atoms with charged labels (positive or negative) contribute to interlayer stability. The orange waves visually depict the electrostatic interaction between P and N atoms, highlighting its key role in stabilizing the system. For all structures, blue spheres represent N atoms and red spheres P atoms.
Structures alongside their respective phonon dispersions for (a) the configuration possessing the lowest energy and (b) a representative dynamically unstable configuration with 12.5% doping concentration. The black arrows in (b) indicate the direction of vibration for the Z-point imaginary frequency. (c) Variation of energy as a function of interlayer P–P distances. Each circle represents a unique configuration, with the red triangle indicating the configuration with the lowest energy and the blue triangle one of the dynamically unstable configurations. (d) Schematic of interlayer electrostatic interaction for configurations with N–P dipoles. Atoms with charged labels (positive or negative) contribute to interlayer stability. The orange waves visually depict the electrostatic interaction between P and N atoms, highlighting its key role in stabilizing the system. For all structures, blue spheres represent N atoms and red spheres P atoms.
Due to the difference in electronegativity, there is a charge transfer from P atoms to N atoms, resulting in P atoms carrying a positive charge and N atoms a negative charge. For example, in the configuration with the lowest energy [Fig. 4(a)], the average Bader charge of the P atom and N atom are 3.10 and 5.27 e/atom, respectively. In detail, the N atom nearest to a P atom has a Bader charge of 6.42 e, and on moving farther away from the P atom, the charges of the N atom decrease and approach 5.0 e. Furthermore, as shown in Fig. 4(d), the N atoms that bonded to P atoms form dipoles with the adjacent layer of P atoms, contributing to stable electrostatic interactions. However, in configurations with imaginary frequencies, the absence of these dipoles will cause interlayer sliding.
In the case of a 25% doping concentration, we conducted an extensive analysis of all phonons for stable configurations and found four configurations without imaginary frequencies. In these configurations, the P atoms between the layers are arranged perpendicularly, and they all exhibit interlayer N–P dipoles. For the other configurations, imaginary frequencies appear at the Z point. This result is consistent with the case of 12.5% doping concentration. Figure 5 depicts the two dynamically stable configurations, along with their respective phonon dispersions. In Fig. 5(a), a unit cell contains four dipoles, and the lowest frequency at the Z point is 2.524 THz, while in Fig. 5(b), there are two dipoles, and the lowest frequency at the Z point is 0.954 THz. When the number of N–P dipoles decreases, the lowest frequency at the Z point will also decrease. It is worth mentioning that for configurations in which the P atoms between the layers are not arranged perpendicularly, we observed a tendency for these P atoms to align perpendicularly after structural relaxation in order to form more dipoles.
Structures and corresponding phonon dispersions for 25% doping concentration. (a) Configuration with four N–P dipoles. (b) Configuration with two N–P dipoles. The orange waves visually depict the electrostatic interaction between P and N atoms. Blue spheres represent N atoms and red spheres P atoms, and atoms with charged labels contribute to interlayer stability.
Structures and corresponding phonon dispersions for 25% doping concentration. (a) Configuration with four N–P dipoles. (b) Configuration with two N–P dipoles. The orange waves visually depict the electrostatic interaction between P and N atoms. Blue spheres represent N atoms and red spheres P atoms, and atoms with charged labels contribute to interlayer stability.
The decomposition enthalpies of 12.5% and 25% doped BP-N relative to α′-P3N569 and α-N227 at ambient pressure are 1.34 and 1.18 eV/atom, respectively. The considerable decomposition enthalpy is attributed to the fact that P-doped BP-N is mainly composed of N–N single bonds, while the primary decomposition product is N2, which consists of N≡N triple bonds. There exists a substantial energy difference between N–N single and N≡N triple bonds. Therefore, doped BP-N is a potential HEDM. To further investigate their energy density and explosive performance, we conducted additional studies and compared the results with experimental values for the known trinitrotoluene (TNT), 1,3,5,7-tetranitro-1,3,5,7-tetrazocane (HMX),76–78 as depicted in Fig. 6. Subsequent calculations estimated the gravimetric energy density Ed to be ∼8.07 kJ/g for a 12.5% doping concentration (density 3.06 g/cm3) and 6.27 kJ/g for a 25% doping concentration (density 3.13 g/cm3). These values are significantly higher than those for TNT and HMX. Notably, the Ed for 12.5% doping concentration is approximately twice as high as that of TNT (4.3 kJ/g). Calculations of the volumetric energy density Ev further highlight the superior performance of 12.5% doping concentration, with a value of 24.71 kJ/cm3. This is more than three times higher than the Ev of TNT (7.05 kJ/cm3).
Calculated density ρ, energy density Ed, volumetric energy density Ev, detonation velocity Vd, and detonation pressure Pd of 12.5% and 25% P-atom-doped BP-N. For comparison, the values for TNT and HMX explosives are also listed. Superscript “Expt.” indicates experimental data.
Calculated density ρ, energy density Ed, volumetric energy density Ev, detonation velocity Vd, and detonation pressure Pd of 12.5% and 25% P-atom-doped BP-N. For comparison, the values for TNT and HMX explosives are also listed. Superscript “Expt.” indicates experimental data.
Furthermore, as illustrated in Fig. 6, we estimated the explosive performance of 12.5% and 25% doping concentrations by utilizing the Kamlet–Jacobs empirical equation.72 The validity of this equation is widely accepted. The detonation velocities Vd of 12.5% and 25% doping concentrations were calculated to be 18.07 and 15.02 km/s, respectively, while the corresponding pressures Pd were 1885.19 and 1314.39 kbar. Notably, the 12.5% doping concentration gave Vd values two to three times higher than that of TNT (6.90 km/s) and Pd values approximately ten times higher than that of TNT (190 kbar).
IV. CONCLUSION
A systematic investigation into the stabilities of BP-N using first-principles calculations has been carried out. It has shed light on the effects of various factors such as doping concentration, doping uniformity, and doping configuration. P-atom doping can induce BP-N to the ground state at high pressure. Furthermore, a tendency toward uniform dopant distribution can improve the structural stability of BP-N, which remains stable even at 0 GPa. By introducing P-atom dopants, the adhesion between layers is strengthened through interlayer electrostatic interactions of N–P dipoles, which reduces the tendency for interlayer slip and enhances the stability of the structure. Furthermore, pressure is conducive to enhancing the stability of BP-N and its doped forms by preventing decomposition into N-chains. The calculated Ed and Ev of 12.5% doped BP-N are 8.07 kJ/g and 24.71 kJ/cm3, which are respectively about two and three times higher than those of TNT. This highlights the significance of doping as an effective approach to achieve stable and high-performance BP-N.
SUPPLEMENTARY MATERIAL
The supplementary material encompasses a phonon spectrum diagram and structural diagrams of the doped structures. These figures provide additional support and visual explanations for the research findings.
ACKNOWLEDGMENTS
This work is supported by the National Natural Science Foundation of China (NSFC) under Grant No. U2030114 and the CASHIPS Director’s Fund (Grant No. YZJJ202207-CX). The calculations were partly performed at the Center for Computational Science of CASHIPS, the ScGrid of the Supercomputing Center and Computer Network Information Center of the Chinese Academy of Sciences, and the Hefei Advanced Computing Center.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Guo Chen: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Chengfeng Zhang: Investigation (supporting). Yuanqin Zhu: Investigation (supporting). Bingqing Cao: Investigation (supporting). Jie Zhang: Formal analysis (supporting); Writing – review & editing (supporting). Xianlong Wang: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).
DATA AVAILABILITY
All data included in this study are available upon request from the corresponding author.