Realized stable BP-N at ambient pressure by phosphorus doping

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 N₂ 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.¹ 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.² 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,¹² layers,^13–16 clusters, chains,^17–19 cages,²⁰ ions,²¹ 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 N₅⁺ ion was first synthesized in 1999,² and there have been several breakthroughs in the synthesis of N₅⁻ since 2016.³³ 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.⁴ 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.³⁶ 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.³⁷ Similar decomposition phenomena at low temperature have been observed in other high-pressure-synthesized polymeric all-nitrogen forms, including cg-N,⁴ LP-N,³⁴ and HLP-N,³⁵ 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.³⁷ 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.⁶⁰ Moreover, BP exhibits promising properties and holds significant potential for various applications.³⁶ 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.

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).⁶¹ The generalized gradient approximation (GGA)^62,63 parameterized by the Perdew–Burke–Ernzerhof (PBE)⁶⁴ exchange–correlation functional and the projected augmented wave (PAW)⁶⁵ method with N-2s²2p³ and P-3s²3p³ 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⁻⁶ 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 Å⁻¹, 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, α′-P₃N₅⁶⁹ and α-N₂²⁷ were chosen as the ground-state phases. VAPKIT⁷⁰ and MULTIWFN⁷¹ were also adopted in our work.

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.

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 

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.

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.

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.

The decomposition enthalpies of 12.5% and 25% doped BP-N relative to α′-P₃N₅⁶⁹ and α-N₂²⁷ 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 N₂, 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 E_d to be ∼8.07 kJ/g for a 12.5% doping concentration (density 3.06 g/cm³) and 6.27 kJ/g for a 25% doping concentration (density 3.13 g/cm³). These values are significantly higher than those for TNT and HMX. Notably, the E_d for 12.5% doping concentration is approximately twice as high as that of TNT (4.3 kJ/g). Calculations of the volumetric energy density E_v further highlight the superior performance of 12.5% doping concentration, with a value of 24.71 kJ/cm³. This is more than three times higher than the E_v of TNT (7.05 kJ/cm³).

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.⁷² The validity of this equation is widely accepted. The detonation velocities V_d of 12.5% and 25% doping concentrations were calculated to be 18.07 and 15.02 km/s, respectively, while the corresponding pressures P_d were 1885.19 and 1314.39 kbar. Notably, the 12.5% doping concentration gave V_d values two to three times higher than that of TNT (6.90 km/s) and P_d values approximately ten times higher than that of TNT (190 kbar).

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 E_d and E_v of 12.5% doped BP-N are 8.07 kJ/g and 24.71 kJ/cm³, 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.

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.

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.

The authors have no conflicts to disclose.

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).

All data included in this study are available upon request from the corresponding author.