Low compressible $β$-BP₃N₆ Open Access

Hardness is an important property that determines many of the technological applications of materials.¹ Design of materials based on the first principles approach, synthesis, and then characterization of these materials are of great interest to both theoretical and experimental material scientists. Well known superhard materials such as diamond and boron nitride display strong covalent bonds, low compressibility, and high wear resistance. The applicability of these materials at high-pressure is limited, and therefore, there is a need for new materials that would work perfectly under those harsh conditions. For many years now, search for new superhard materials has been an order of the day. This search is unlikely to end any time soon due to the fact that diamond remains to be the hardest known material made of a single light element called carbon. The combination of carbon and other light elements such as nitrogen has shown promising results from the theoretical perspective,^2–6 but most of the predicted materials have not been realized experimentally.

The present study has been prompted by the recent explorative investigation of phosphorus nitrides by Vogel et al.,^7,8 who synthesized a high-pressure polymorph of boron phosphorus nitride (β-BP₃N₆). This phase adopts an orthorhombic structure with space group Pna2₁ (no. 33). The β-BP₃N₆ phase is characterized with octahedral coordinated phosphorous (P) atoms in a distorted hexagonal closed-packing of nitride anions. This phase was obtained at a high-pressure experiment having transformed from α-BP₃N₆ at a pressure of about 42 GPa. To our knowledge, no previous first-principles calculations have been carried out to investigate the mechanical stability and the hardness of this new phase of β-BP₃N₆ at different pressures.

The mechanical stability is a concept based on the elastic constants of a single crystal sample of a material. The elastic constants determine the response of the crystal to external forces that can be characterized in forms of bulk modulus, shear modulus, Young’s modulus, and Poisson’s ratio. All these properties can be obtained from density functional theory (DFT),^9,10 which has become an essential tool for designing novel superhard materials.

In this paper, we first benchmarked our study with previous theoretical studies on the lattice parameters, elastic constants, and bulk and shear moduli of hypothetical superhard carbon mononitride (Pnnm-CN).² This compound adopts an orthorhombic crystal structure of space group Pnnm (no. 58), and its Vickers hardness has been predicted to be above 60 GPa. Once satisfied with Pnnm-CN results, we went ahead to investigate the lattice parameters, elastic constants, bulk modulus, shear modulus, Poisson’s ratio, Young’s modulus, Vickers hardness, and mechanical stability of β-BP₃N₆ at different pressures. Both Pnnm-CN and β-BP₃N₆ adopt an orthorhombic structure with a dense network of covalent bonds, which are precursors for low compressibility.

The remaining parts of this paper are as follows: in Sec. II, we give a brief outline of the computational details. Results are shown and discussed in Sec. III. Finally, in Sec. IV, we provide the conclusion.

Calculations in the present work were performed using Density Functional Theory (DFT)^9,10 within the Perdew–Burke Ernzerhof (PBE)¹¹ generalized gradient approximation and plane waves basis sets as implemented in the Quantum ESPRESSO package.¹² All pseudo-potentials were obtained from the pslibrary1.0 maintained by Dal Corso.¹³ The cell volume and ions were optimized at three different pressure values, i.e., 0 GPa, 20 GPa, and 42.4 GPa. The k-point sampling¹⁴ was performed on a 6 × 4 × 4 grid while the kinetic energy cutoff was set to 43 Ry. It is well known that the PBE functional overestimates the lattice parameter and consequently affects not only the volume and the density of the compound but also the mechanical properties of the compound. An external pressure exerted on the crystal structure also has profound influence on the mechanical properties. Therefore, to get physically meaningful results from the calculations of elastic constants, we did not allow the ions to relax. This was necessary to maintain the strain on each ion at different pressures. The thermo_pw package^15,16 was used to run the elastic constant calculations.

In Table I, we have listed the calculated physical properties of Pnnm-CN at zero pressure predicted using the PBE functional. The calculated lattice constants produced from fitting total energy as a function of volume to the Murnaghan¹⁷ equation of state are compared to the corresponding theoretical results. In general, our calculated lattice constants and other properties shown in Table I are in agreement with the results of previous studies. These results gave us confidence to study lattice parameters a (Å), b (Å), c (Å), volume V (ų), and density ρ (g/cm³) of β-BP₃N₆ at zero and elevated pressures.

Figure 1 shows the crystal structure of the primitive cell of β-BP₃N₆, where P are the large atoms labeled P1 and B are the medium-sized atoms labeled B1, whereas the smaller ones are N atoms labeled N1. The unit cell contains 40 atoms and has a space group of Pna2₁ (no. 33).

Table II shows the response of the lattice parameter to different external pressures exerted on β-BP₃N₆. As the exerted pressure increases from zero to 20 GPa, the unit-cell contracts by 1.87%, 2.04%, and 1.79% in a, b, and c axes, respectively. We see a lot more contraction in the b axis than in both a and c axes. This is a clear indication of the existence of different bonding behaviors along the axes. One would expect strong covalent bonding between boron and nitrogen atoms and weaker bonding character in phosphorus-nitrogen bond formation. On increasing the external pressure to 42.4 GPa, the unit-cell contracts by 3.51%, 3.78%, and 3.39% in a, b, and c axes, respectively. Again, the b axis displays slightly more collapse in its distance than in both a and c axes. At 42.4 GPa, the structure of the real crystal obtained from the work of Vogel et al.⁸ is not completely known to support a direct comparison with current theoretical data. Further discrepancies are from a well known concept that the PBE functional overestimates the lattice constants and consequently the volume of the crystal structure. In this case, the predicted volume of β-BP₃N₆ is slightly higher than the experimental values by 3.5%. Furthermore, the volumetric density of a material changes linearly with pressure and inversely with volume. At 42.4 GPa, the predicted density is lower than the experimental value by 165 g mm⁻³. Generally, good agreement is reached between the computed and measured values.

Since 1954, the work of Born¹⁸ has been used to evaluate the elastic stability of crystal structures. The cubic crystals are easy to handle due to their high symmetry. Efforts to apply the Born¹⁸ elastic stability criteria on lower symmetry crystals are reported in the works of Ravindran et al. and those of Mouhat and co-author.^19,20 These authors gave a precise description of the form it should take for lower symmetry crystal classes. They^19,20 particularly considered the orthorhombic crystal structure as one such system with a lower symmetry. They described it as having nine independent elastic constants, namely, C₁₁, C₁₂, C₁₃, C₂₂, C₂₃, C₃₃, C₄₄, C₅₅, and C₆₆. The mechanical stability of such orthorhombic crystals is achieved whenever the elastic constants satisfy the following necessary and sufficient conditions as prescribed by Born:¹⁸ 

More work on the orthorhombic crystal system has been illustrated by Wen et al.²¹ These particular authors applied elastic stability conditions shown in Eq. (1) on orthorhombic TiAl alloys. They underscored the usefulness of Eq. (1) in checking the stability of the orthorhombic crystal system.

In this work, the calculated values of the nine independent elastic constants at different values of pressures (0, 20, and 42.4 GPa) are listed in Table III, and they satisfy the necessary and sufficient conditions prescribed for an orthorhombic system. All nine independent elastic constants increase monotonically with pressure. Already at zero pressure, β-BP₃N₆ depicts a material with large elastic constants that has a strong correlation with the compressibility and hardness of a material. At elevated pressures, shorter covalent bonds will resist change in volume, a response that can be captured as enhanced elastic constants. At 42.4 GPa, β-BP₃N₆ has very large elastic constants comparable to that of diamond. We anticipate that such large elastic constants will attract wide industrial application of β-BP₃N₆ and where possible be a substitute to already existing hard materials. It is noted that the computing properties of a low symmetry system are tricky, and therefore, one has to be careful with the values of C₁₁, C₁₂, C₁₃, C₂₂, C₂₃, C₃₃, C₄₄, C₅₅, and C₆₆ as will always depend on the orientation of the unit cell. To our knowledge, there are no data on elastic constants to compare with our results, and hence, our work lays a foundation for future references.

From the single crystal elastic constant data,²¹ the polycrystalline bulk modulus B and shear modulus G can be calculated by using the Voigt approximation,²² which describes the upper bound, and the Reuss approximation,²³ which describes the lower bound. The average of the two bounds is often described by using the Hills approximation.²⁴ Below are the equations that relate the nine independent elastic constants to B_V and G_V in Voigt notations, whereas B_R and G_R represent Reuss notations,

For simplicity, the averages of Voigt and Reuss approximations are represented as B and G for bulk and shear moduli, respectively,

The derived quantities from B and G are presented as Young’s modulus (E) and Poisson’s ratio (n),

Young’s modulus is a measure for the stiffness of a material, while Poisson’s ratio is often used to classify materials as either ductile or brittle.²⁵ This concept by Haines et al. will normally treat a brittle material as the one having n below 0.33, while a ductile material as the one that has n greater than 0.33,

On the other hand, Pugh’s²⁶ definition of a ductile/brittle material is based on the G/B ratio. Ductility is measured as a value of G/B less than 0.5, while a material is rated brittle if G/B goes beyond 0.5. Superhard phases of hypothetical C₃N₄ and diamond are good examples of brittle materials as demonstrated in the previous work.⁴ 

For this work, the calculated bulk modulus (B), shear modulus (G), Young’s modulus (E), and Poisson’s ratio (n) of β-BP₃N₆ in Voigt’s and Reuss’s approaches are listed in Table IV. The bulk modulus, which is also referred to as incompressibility, is a ratio that relates an applied constant stress to the fractional volumetric change. β-BP₃N₆ has a bulk modulus of about 330 GPa comparable to many other low-compressible compounds, such as SiO₂ with a bulk modulus of 305 GPa, HfN with a bulk modulus of 306 GPa, and OsB₂ with a bulk modulus of 297 GPa.^27–30 The low compressibility of β-BP₃N₆ can be attributed to a highly condensed network built up from the tetrahedra BN₄ and octahedra PN₆ coordination. Although it has been argued that incompressibility is a good indicator of the hardness of a material, it is the shear modulus that is strongly associated with hardness. The shear modulus relates to how a material responds to the change in shape at constant volume. The calculated shear modulus of β-BP₃N₆ is about 313 GPa and is in the same range as that of Pnnm-CN, which is reported as 322 GPa; see Tables IV and II, respectively. The current value of shear modulus (313 GPa) does not compete that of diamond (546 GPa) as reported elsewhere.^3,4 It is worth noting that shear modulus increases with increasing pressure. This implies that a material has high resistance to shape change when it is under the influence of an external pressure. As indicated earlier, Young’s modulus is a measure for the stiffness of a material. It is derived from both shear and bulk moduli. Computed values of β-BP₃N₆’s Young’s moduli for different pressures are presented in Table IV. The calculated Poisson’s ratio suggests that β-BP₃N₆ should be classified as a brittle material since the value of n is less than 0.33 even at high pressures. Similarly, the G/B ratio calculated as 0.94 at zero pressure, 0.86 at 20 GPa, and 0.81 at 42.4 GPa predicts β-BP₃N₆ as a brittle material with G/B greater than 0.5. The values of B/G for β-BP₃N₆ are shown in Table V. There are tendencies of β-BP₃N₆ to transform from being brittle to ductile when exposed to extremely high pressures, as depicted by the diminishing G/B ratio at elevated pressures. This trend can also be interpreted to mean that BP₃N₆ will undergo pressure-induced phase transformation to lower symmetry crystal structures such as triclinic at extremely high pressures. The isotropy of the crystal structure is very critical in the characterization of any compound. The single crystal shear anisotropy factors A_{100} in {100} planes, A_{010} in {010} planes, and A_{001} in {001} planes are defined as²⁰ 

The percentage anisotropy in compressibility A_B and shear A_G is defined as³¹ 

A zero value of A_B and A_G corresponds to elastic isotropy, while a value of 100% corresponds to the largest possible anisotropy. The universal anisotropy index (A^U) is defined as³² 

A nonzero value of A^U is a measure of the anisotropy. The shear anisotropy factors (A_{100}, A_{010}, A_{001}), percent anisotropy factors of bulk (A_B) and shear (A_G) moduli, and universal anisotropy factor A^U of β-BP₃N₆ are presented in Table V. For the isotropic case, we expect A_B, A_G, and A^U to be zero and any deviation from zero is a clear indicator for anisotropy. Therefore, β-BP₃N₆ is anisotropic with highest deviations observed in both A_{100} and A_{001} directions. The A_{100} direction is more isotropic than other directions. In the Introduction, it was stated clearly that hardness is an important property that determines many of the technological applications of materials.¹ The Vickers hardness in this work was estimated on the basis of the so-called Chen model³⁰ that has been widely used to predict Vickers hardness of a variety of crystalline metals, insulators, and semiconductors. Chen’s model is given as

where G and B are the shear and bulk modulus, respectively. The Vickers hardness of β-BP₃N₆ is given in Table V. Clearly, β-BP₃N₆ has high Vickers hardness close to that of Pnnm-CN and hence a potential superhard material. The strong resistance to change in volume and shape as a result of highly condensed orthorhombic structures explains this unexpected behavior in β-BP₃N₆. Right from large elastic constants to large bulk and shear moduli, small values of Poison’s ratio (ν), and large Vickers harness, β-BP₃N₆ portrays similar behavior to well known superhard materials. However, pressure frustrates the hardness of β-BP₃N₆.

β-BP₃N₆ has been investigated at different pressures starting from zero pressure to a maximum of 42.4 GPa. Clearly, the high value of bulk modulus indicates that β-BP₃N₆ is a low compressible material. This character emanates from the highly condensed orthorhombic structure built up from the tetrahedra BN₄ and octahedra PN₆ coordination. A comprehensive analysis of elastic constants and its derived properties besides bulk modulus have been presented. β-BP₃N₆ displays high resistance to change in shape, which greatly contributes to its high Vickers hardness of about 51 GPa at zero pressure. Overall, β-BP₃N₆ has a similar behavior as those of well known superhard materials, and hence, the conclusion is that β-BP₃N₆ is a potential superhard material. These results provide guidance for further exploration of thermodynamic stability and many other properties of β-BP₃N₆ at different pressures.

This work was financially supported by the Kenya Education Network (KENET) through the Computational Modeling and Materials Science (CMMS) Research mini-grants 2019. We acknowledge the Centre for High Performance Computing (CHPC), Cape Town, South Africa, for providing us with computing facilities.