## Abstract

Boron nitride (BN) is a material with outstanding technological promise due to its exceptional thermochemical stability, structural, electronic, and thermal conductivity properties, and extreme hardness. Yet, the relative thermodynamic stability of its most common polymorphs (diamond-like cubic and graphite-like hexagonal) has not been resolved satisfactorily because of the crucial role played by kinetic factors in the formation of BN phases at high temperatures and pressures (experiments) and by competing bonding and electrostatic and many-body dispersion forces in BN cohesion (theory). This lack of understanding hampers the development of potential technological applications and challenges the boundaries of fundamental science. Here, we use high-level first-principles theories that correctly reproduce all important electronic interactions (the adiabatic-connection fluctuation-dissipation theorem in the random phase approximation) to estimate with unprecedented accuracy the energy differences between BN polymorphs and thus overcome the accuracy hurdle that hindered previous theoretical studies. We show that the ground-state phase of BN is cubic and that the frequently observed hexagonal polymorph becomes entropically stabilized over the cubic at temperatures slightly above ambient conditions (*T*_{c→h} = 335 ± 30 K). We also reveal a low-symmetry monoclinic phase that is extremely competitive with the other low-energy polymorphs and that could explain the origins of the experimentally observed “compressed h-BN” phase. Our theoretical findings therefore should stimulate new experimental efforts in bulk BN and promote the use of high-level theories in modeling of technologically relevant van der Waals materials.

## INTRODUCTION

Despite the tremendous technological interest of bulk boron nitride (BN) (*1*, *2*), fundamental knowledge of its phase diagram remains contentious to this day. The two most common BN polymorphs have hexagonal (h-BN) and cubic (c-BN) symmetries and are structurally analogous to the graphite and diamond phases of carbon, respectively (Fig. 1). On the basis of empirical observations, and in analogy to the carbon phase diagram, h-BN is generally regarded as the most stable BN polymorph at ambient conditions (*3*–*5*): c-BN does not exist in nature, and its synthesis in laboratories requires high-temperature and high-pressure conditions. Notably, experimental phase diagrams based on thermodynamic and in situ x-ray diffraction measurements strongly suggest that c-BN is more stable than h-BN under normal conditions (*6*–*8*). The reported c-BN↔h-BN transition temperatures as extrapolated to ambient pressure, however, vary by as much as (*8*). The cause of this huge variation is the critical importance of kinetic effects on the c-BN↔h-BN transformation, which depends strongly on difficult-to-control parameters such as grain size, defect concentration, and purity of the starting material (*7*, *8*).

Calculations based on quantum mechanics are free of the abovementioned kinetic factors affecting experiments; however, weakly bound layered materials, such as most BN polymorphs, are known to pose serious challenges to standard first-principles methods, which do not include dispersion (van der Waals) interactions [e.g., density functional theory (DFT) based on the local density approximation (LDA) (*9*) and generalized gradient approximation (GGA) (*10*) to the exchange-correlation energy]. DFT estimations based on LDA and GGA reach contrary conclusions on the relative stability of the c-BN and h-BN polymorphs [see, for instance, (*11*, *12*)], thus adding further confusion to the BN phase diagram puzzle. The last decade has seen extraordinary progress in the development of dispersion-corrected theories that overcome the limitations of standard theories; modern approaches such as the D3 empirical correction (*13*–*15*) and many-body dispersion (MBD) methods (*16*, *17*) are able to reliably optimize complex layered polymorph structures and predict energy differences between them with fair accuracy (*18*). However, even these recent developments might not be relied upon to determine the relative thermodynamic stability of low-energy BN polymorphs since the involved energy differences can be below ~ 1 kJ/mol (~ 10 meV per formula unit), that is, the characteristic scale of nonsystematic errors in most dispersion approximations.

The challenge of capturing these small energy differences is made more difficult by the importance of many-body interactions in dispersion-bound systems. Low-dimensional systems, such as the layered polymorphs of BN, exhibit collective dispersion interactions (*16*, *19*, *20*) identified as “Type B” nonadditivity by Dobson (*21*). Type B effects cannot be represented as a sum over pairwise interactions, as is done in many dispersion correction schemes or even in more sophisticated perturbation approaches such as the Møller-Plesset theory. This makes BN polymorph ranking especially challenging, as any theoretical model must incorporate type B terms in its underlying physics.

Here, we use the adiabatic-connection fluctuation-dissipation theorem in the random phase approximation (RPA), a method that has been shown to reliably and seamlessly treat both strong and weak interactions in a wide variety of systems and has a full treatment of many-body interactions (*22*–*26*). We thereby calculate the energies of BN polymorphs with unprecedented accuracy to overcome previous theoretical bottlenecks. We show that the ground state of BN is c-BN and that the h-BN polymorph becomes thermodynamically most stable at temperatures close to ambient, namely, *T*_{c→h} = 335 ± 30 K. By using a multistage approach (Methods, Supplementary Methods, and table S1) to determine and rank low-energy structures from an initial list of 15 structures, we reveal a low-symmetry monoclinic phase that turns out to be energetically very competitive with other known polymorphs and could explain the origins of previously overlooked experimental observations. We discuss the causes of the phase stability phenomena revealed in bulk BN and argue the need for using high-level theories in computational studies of technologically relevant materials.

## RESULTS

Figure 1 shows the zero-temperature energies of most popular BN polymorphs relative to that of h-BN(*AA*′), calculated with several dispersion-corrected DFT methods and RPA (corresponding space groups and polymorph notation used throughout the text are explained therein). These zero-temperature energies include quantum nuclear effects (*27*) via zero-point energy (ZPE) corrections (Methods), although the ZPE effects do not lead to qualitative changes in the relative stability of BN polymorphs (table S2). Our results demonstrate two key points. First, the RPA ordering of low-energy states confirms that c-BN is the ground state (*6*–*8*), which is about 1 to 2 kJ/mol lower in energy than the h-BN(*AA*′), m-BN (discussed below), h-BN(*AB*), and h-BN(*ABC*) (also known as r-BN) polymorphs. Second, the MBD method (*16*) agrees quite well with RPA, this being the only semiempirical method that predicts correctly the energy ordering among all the most competitive phases. Meanwhile the D3 method with Becke-Johnson damping [D3(BJ)] (*15*) and the fractionally ionic (FI) MBD (*17*) method come close and the energy results obtained with the LDA and GGA methods turn out to be of unsatisfactory quality (table S2).

Our study reveals a low-symmetry monoclinic phase, denoted here as m-BN [space group *Cm*, different from a previously predicted monoclinic phase (*28*), hence the “2” in Fig. 1], that is energetically very competitive with respect to the h-BN polymorphs (Δ*E*^{RPA} < 1 kJ/mol). This new phase is vibrationally stable and presents a reduced two-atom unit cell with equilibrium parameters *a*_{m} = 4.34 Å, *b*_{m} = 2.51 Å, *c*_{m} = 3.56 Å, and β = 113°, as obtained with MBD and FI methods (Fig. 2, table S3, and fig. S1). The predicted m-BN polymorph is structurally very similar to a previously reported monoclinic phase that was experimentally observed during the c-BN↔h-BN transformation occurring at high-*P* high-*T* conditions (, , to 3.3 Å, and β^{exp} = 92*°* to 95°) and that was named as “compressed h-BN” (*29*); hence, we tentatively identify the two phases as the same (up to the non-negligible effects of pressure and temperature disregarded in our simulations).

In Table 1, we show the elastic constants, *C*_{ij} (given in Voigt notation), of the low-energy polymorphs calculated with the RPA and LDA methods. The elastic properties of the monoclinic phase are very similar to those of h-BN(*AA*′), a result that, along with the minute energy difference among the two, might explain the causes of the experimentally observed compressed h-BN/h-BN coexistence (table S4) (*29*). We note that the series of *C*_{ij} values estimated with the LDA and RPA methods are in very good agreement, with key LDA values systematically 5 to 6% above their RPA counterparts. This outcome supports the use of LDA for accurately assessing the vibrational properties of BN polymorphs, a task that has a prohibitive computational cost in RPA.

Figure 3 shows the Gibbs free energy *G* of several BN polymorphs estimated at zero pressure with the quasi-harmonic approach (QHA) (Methods). Static/vibrational contributions to *G* are calculated with the RPA/LDA method (Methods). Temperature-induced volume expansion effects are fully accounted for in our results to provide precise phase transition data (Methods and fig. S2). We find that the h-BN(*AA*′) polymorph becomes entropically stabilized over c-BN at *T*_{c→h} = 335 K, a temperature relatively close to ambient conditions that falls substantially below the corresponding experimental data (*6*–*8*). Taking into account a possible numerical error of 0.2 kJ/mol in the RPA calculations and an additional 0.06 kJ/mol in the vibrational free energies (Methods) leads to *T*_{c→h} = 335 ± 30 K, still below the lowest experimental result. We note, however, that the mass densities that we estimate (using FI + LDA corrections) for h-BN(*AA*′) and c-BN at *T*_{c→h} agree almost perfectly with the experimental measurements carried out at room temperature (Fig. 3) (*30*, *31*). The low-energy phonon excitations in the h-BN(*AA*′) phase present much lower frequencies than those in the c-BN (fig. S3); hence, the vibrational entropy of the former polymorph becomes increasingly more favorable as the temperature is raised. Meanwhile, the Gibbs free energy of the wurtzite polymorph, w-BN, falls out of the energy range considered in Fig. 3 due to its extreme vibrational stiffness, which translates into destabilizing entropy contributions under increasing temperature. By contrast, the Gibbs free energy of the m-BN, h-BN(*AB*), and h-BN(*ABC*) (also known as r-BN) polymorphs follows closely that of h-BN(*AA*′), all falling within an energy difference range of about 1 kJ/mol at 300 ≤ *T* ≤ 400 K.

## DISCUSSION

Our results confirm that the c-BN polymorph is most stable at low pressures; however, the transition temperature that we predict for the c-BN↔h-BN transformation lies relatively close to ambient conditions. Consequently, the phase diagram of BN, if not analogous, is certainly more similar than previously thought to that of carbon. Actually, our *T*_{c→h} estimation appears to be consistent with the general belief based on empirical observations that h-BN is most stable at ambient conditions. The likely reason for the substantial difference between theory and measurements, , may be the high-temperature high-pressure conditions and kinetic factors involved in the experimental synthesis and analysis of BN samples (*8*). We hope that our theoretical results will motivate new experimental activity in bulk BN. On the theory side, we have shown that (i) LDA is a good method for calculating elastic properties of materials and, hence, is probably good for estimating vibrational free energies, (ii) but for accurate prediction of energy ordering among van der Waals polymorphs, one must use methods that include MBD interactions, ideally at a high level using RPA, but certainly approximately when RPA is infeasible.

## METHODS

### DFT and phonon calculations

First-principles calculations based on DFT were performed to analyze the energy and structural and vibrational properties of BN polymorphs. We performed these calculations with the VASP (Vienna Ab initio simulation package) code (*32*) using projector augmented-wave method potentials (*33*). The electronic states 1*s*-2*s* of B and 2*s*-2*p* of N atoms were considered as valence. Wave functions were represented in a plane-wave basis truncated at 650 eV. By using these parameters and dense **k**-point grids for the integration within the first Brillouin zone (IBZ), energies were converged to within 1 meV per formula unit (0.1 kJ/mol; fig. S4). In the geometry relaxations, a tolerance of 0.01 eV Å^{−1} were imposed in the atomic forces.

Ab initio phonon frequencies were calculated with the direct method to assess the vibrational stability of the analyzed BN polymorphs and estimate their Gibbs free energies as a function of temperature and pressure within the QHA (*27*). In the direct method, the force-constant matrix was calculated in real space by considering the proportionality between atomic displacements and forces (*34*). The quantities with respect to which our phonon calculations were converged include the size of the supercell, the size of the atomic displacements, and the numerical accuracy in the sampling of the IBZ. We found the following settings to provide quasi-harmonic free energies converged to within 0.1 kJ/mol: 4 × 4 × 4 supercells (where the figures indicate the number of replicas of the unit cell along the corresponding lattice vectors; fig. S4), atomic displacements of 0.02 Å, and **q**-point grids of 14 × 14 × 14. The value of the phonon frequencies were obtained with the PHON code developed by Alfè (*34*). By using this code, we exploited the translational invariance of the system to impose the three acoustic branches to be exactly zero at the center of the Brillouin zone and applied central differences in the atomic forces.

### RPA calculations

To overcome the accuracy barrier of semiempirical theories, we carried out RPA calculations, which give comparable results to high-level coupled-cluster and related wave function theories but are valid for bulk systems with small or zero gaps (*22*–*26*). Because of its high numerical cost, our RPA calculations were performed using structures optimized at the FI level (*17*). Two sets of calculations were carried out. First, exact exchange (EXX) and RPA correlation energy calculations(1)were carried out for all low-energy structures. These used an energy cutoff of 480 eV, a dense 12 × 12 × 12 (or equivalent) **k**-point grid for EXX, and a coarser 7 × 7 × 7 one for RPA, evaluated on self-consistent PBE orbitals as per standard practice (*22*). To further refine the energies of the c-BN, h-BN(*AA*′), and m-BN structures, we carried out additional calculations using a more accurate cutoff of 550 eV and the dense grid for both EXX and RPA; this procedure yields energy difference results within 2 meV per formula unit (0.2 kJ/mol) of the initial calculations, which we use as our numerical error bar. Errors in energy differences between various h-BN phases and the m-BN phase are expected to be much smaller (< 0.2 meV), due to the similarity of the systems and consequent additional error cancellations.

### Estimation of thermodynamic quantities

We used the QHA (*27*) to calculate the Gibbs free energy *G* of BN polymorphs as a function of temperature and pressure. Anharmonic effects beyond the QHA have been shown to be negligible for bulk BN at temperatures close to ambient conditions (i.e., below 0.1 kJ/mol) (*11*); hence, we disregarded them here. [We should note that in the unlikely case that anharmonicity played a role at *T* ~ 300 K, it would probably tend to further stabilize the hexagonal polymorph over the cubic (*11*), thus additionally reducing *T*_{c→h}.] In the QHA approximation, the vibrational free energy of a crystal *F*_{vib} with volume *V* and at temperature *T* is(2)where *N*_{q} is the total number of wave vectors used for integration within the first Brillouin zone, the summation runs over all wave vectors ** q** and phonon branches

*s*, and ω

_{qs}are the vibrational frequencies of the crystal, which depend on volume. In the zero-temperature limit,

*F*

_{vib}becomes(3)which is usually referred to as the “zero-point energy“ (ZPE). The Gibbs free energy of a crystal then reads(4)where

*E*

_{el}is the energy of the system when all atoms rest immobile in their equilibrium positions, and the hydrostatic pressure

*P*is estimated via the volume derivative(5)

Last, by using the thermodynamic constraint *P*(*V*_{0}, *T*) = 0 and performing *E*_{el} and *F*_{vib} calculations over a dense grid of volume points, it is possible to account precisely for *T*-induced volume expansion effects at zero pressure (fig. S2). The zero-temperature energies reported in this study account for possible quantum nuclear effects by means of the expression(6)where *m* denotes the method of calculation, denotes the resulting equilibrium volume, and denotes the ZPE as given by Eq. 3. We checked that the value of ZPE differences between BN polymorphs is practically independent of the used method (fig. S5), hence the reason for our fixed choice of *E*_{ZPE} in Eq. 6. In the RPA case, given the huge computational expense associated with this method, the Gibbs free energies have been estimated by using both *F*_{vib} and hydrostatic pressure values obtained with the LDA method.

This particular choice is justified by the fact that LDA often performs similarly to RPA for stress tensors (*35*), as we have explicitly corroborated in this study (Table 1). On the basis of numerical errors of 6% for the elastic constants of LDA, versus RPA, we assigned a corresponding numerical error to the vibrational free energies of ±0.06 kJ/mol (calculated as 6% of the ZPE difference between the cubic and hexagonal BN polymorphs; Supplementary Methods).

We must note, however, that this good agreement is a result of LDA’s ability to predict the energies of small crystal perturbations, as required for phonon calculations and elastic coefficients. It does not transfer to the energies of structurally distinct systems required for accurate polymorph prediction, which explains LDA’s failures in that regard (table S2).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/1/eaau5832/DC1

Supplementary Methods

Fig. S1. Phonon spectrum of the new monoclinic phase m-BN (space group *Cm*) reported in this study as calculated with the LDA method at zero pressure.

Fig. S2. Gibbs free energy differences among several BN polymorphs calculated at zero pressure and expressed as a function of temperature.

Fig. S3. Phonon spectrum of the c-BN (space group ) and h-BN(*AA′*) (space group *P*6_{3}/*mmc*) polymorphs at zero pressure as calculated with the LDA method.

Fig. S4. Convergence tests of the electronic and vibrational free energies calculated with standard DFT methods (i.e., PBE-D3) in the h-BN(AA′) polymorph.

Fig. S5. ZPE corrections accounting for zero-temperature quantum nuclear effects (*27*) in several BN polymorphs as calculated with the LDA (*9*) and PBE-D3 (*13*, *14*) methods at zero pressure.

Table S1. Numerical tests performed for stage 3 in the RPA calculations.

Table S2. Zero-temperature electronic energies (i.e., neglecting zero-point energy corrections) of several BN polymorphs as compared to that of the c-BN (space group ) phase.

Table S3. Structural properties of the monoclinic phase m-BN (space group *Cm*) reported in this study as calculated with different methods based on DFT.

Table S4. Elastic constants associated with compressive deformations *C*_{ij}’s, bulk modulus, *B*_{VRH}, shear modulus, *G*_{VRH}, isotropic Poisson’s ratio, μ_{VRH}, and Young’s modulus, *Y*_{VRH}, calculated in the Voigt-Reuss-Hill approximation [“VRH,” as this is appropriate for polycrystalline samples; see previous work (*30*) for the corresponding analytical expressions] with the LDA method.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This research was supported under the Australian Research Council’s Future Fellowship funding scheme (no. FT140100135). Computational resources and technical assistance were provided by the Australian Government and the Government of Western Australia through Magnus under the National Computational Merit Allocation Scheme and the Pawsey Supercomputing Centre and by Gowonda high-performance computing facilities.

**Author contributions:**C.C. and T.G. designed the research. T.G. carried out DFT and RPA energy calculations. C.C. carried out the structural and thermodynamic analysis. C.C. and T.G. wrote the manuscript. Both authors contributed equally to the present work.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability**

**:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).