Polymorphism of Bulk Boron Nitride

TL;DR

This study uses advanced first-principles calculations to accurately determine the stability and phase transitions of boron nitride polymorphs, revealing the cubic phase as the ground state and identifying a competitive monoclinic phase.

Contribution

It applies high-level many-body theories to resolve longstanding uncertainties about BN polymorph stability and phase behavior, providing new insights into its thermodynamic properties.

Findings

Cubic BN is the ground-state phase.

Hexagonal BN becomes entropically stabilized above 335 K.

A low-symmetry monoclinic phase is identified as energetically competitive.

Abstract

Boron nitride (BN) is a material with outstanding technological promise because of 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, 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_{\rm c \to h} = 335 \pm 30$ K). We also reveal a low-symmetry monoclinic phase that is extremely competitive with the other low-energy polymorphs and which could explain the origins of the experimentally observed "compressed h-BN" phase. Our theoretical findings therefore should stimulate new experimental efforts in bulk BN as well as promote the use of high-level theories in modelling of technologically relevant van der Waals materials.

Figures

Tables

Table 1. Table 1: Elastic constants and bulk modulus in the Voigt-Reuss-Hill approximation (“VRH”, as this is appropriate for polycrystalline samples ivanovskii12 ) of the most stable polymorphs calculated with the LDA lda and RPA rpa0 ; rpa1 ; rpa2 ; rpa3 ; rpa4 methods. LDA results are very close to those from RPA, showing the suitability of LDA for estimating second energy derivatives. Results are in units of GPa.

		$C_{11}$	$C_{22}$	$C_{33}$	$C_{12}$	$C_{23}$	$C_{13}$	$B_{VRH}$
c–BN	LDA	997	997	997	101	101	101	402
c–BN	RPA	968	968	968	81	81	81	378
h–BN($A{A}^{\prime }$)	LDA	923	926	28	174	3	3	138
h–BN($A{A}^{\prime }$)	RPA	910	915	29	153	-4	-5	131
m–BN	LDA	751	839	45	156	100	68	160
m–BN	RPA	745	830	43	137	90	58	151

Equations

$$E_{\mathrm{RPA}}=$$

$$F_{\mathrm{vib}}(V,T)=\frac{1}{N_q}{k}_{B}T\sum _{\mathbf{q}s}\ln \left[2\sinh \left(\frac{\hslash {\omega }_{\mathbf{q}s}}{2k_{\mathrm{B}}T}\right)\right],$$

$$E_{\mathrm{ZPE}}=\frac{1}{N_{\mathrm{q}}}\sum _{\mathbf{q}s}\frac{1}{2}\hslash {\omega }_{\mathbf{q}s},$$

$$G(V,T)=E_{\mathrm{el}}(V)+F_{\mathrm{vib}}(V,T)+PV,$$

$$P(V,T)=-\left(\frac{\partial E_{\mathrm{el}}}{\partial V}+\frac{\partial F_{\mathrm{vib}}}{\partial V}\right).$$

$$E^m(V_0^m)=E_{\mathrm{el}}^m(V_0^m)+E_{\mathrm{ZPE}}^{\mathrm{LDA}}(V_0^{\mathrm{LDA}}),$$

Full text

Polymorphism of Bulk Boron Nitride

Claudio Cazorla

School of Materials Science and Engineering, UNSW Australia, Sydney NSW 2052, Australia

  

Tim Gould

Queensland Micro- and Nanotechnology Centre, Griffith University, Nathan, QLD 4111, Australia

School of Environment and Science, Griffith University, Nathan, QLD 4111, Australia

**Boron nitride (BN) is a material with outstanding technological promise because of 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, 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_{\rm c\to h}=335\pm 30$ K). We also reveal a low-symmetry monoclinic phase that is extremely competitive with the other low-energy polymorphs and which could explain the origins of the experimentally observed “compressed h–BN” phase. Our theoretical findings therefore should stimulate new experimental efforts in bulk BN as well as promote the use of high-level theories in modelling of technologically relevant van der Waals materials.

**

In spite of the tremendous technological interest of bulk boron nitride (BN) dean10 ; tang02 , fundamental knowledge of its phase diagram remains contentious to this day. The two most common BN polymorphs possess hexagonal (h–BN) and cubic (c–BN) symmetries and are structurally analogous to the graphite and diamond phases of carbon, respectively (Fig. 1). Based on empirical observations, and in analogy to the carbon phase diagram, h–BN generally is regarded as the most stable BN polymorph at ambient conditions pease50 ; corrigan75 ; narayan16 : c–BN does not exist in nature and its synthesis in laboratories requires high-temperature and high-pressure conditions. Strikingly, experimental phase diagrams based on thermodynamic and in situ x-ray diffraction measurements strongly suggest that c–BN is more stable than h–BN at normal conditions bundy62 ; solozhenko99 ; will00 . The reported c–BN$\leftrightarrow$h–BN transition temperatures as extrapolated to ambient pressure, however, vary by as much as $420<T^{\rm exp}_{\rm c\to h}<1500$ K will00 . The cause of this huge variation is the critical importance of kinetic effects on the c–BN$\leftrightarrow$h–BN transformation, which depends strongly on difficult-to-control parameters like grain size, defects concentration, and the purity of the starting material solozhenko99 ; will00 .

Calculations based on quantum mechanics are free of the abovementioned kinetic factors affecting experiments; however, weakly bound layered materials, like 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 – LDA – lda and generalized gradient – GGA – pbe approximations to the exchange-correlation energy). In fact, 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, kresse99 ; ahmed07 ), 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 like the D3 empirical correction d31 ; d32 ; d3bj and many-body dispersion methods mbd ; fi are able to reliably optimize complex layered polymorph structures, and predict energy differences between them with fair accuracy tawfik2018 . But 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 $\sim 1$ kJ/mol ($\sim 10$ meV per formula unit), that is, the characteristic scale of non-systematic 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, like the layered polymorphs of BN, exhibit collective dispersion interactions mbd ; ambrosetti2016 ; hermann2017 identified as “Type B” non-additivity by Dobson dobson2014 . 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 like 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 which has a full treatment of many-body interactions rpa0 ; rpa1 ; rpa2 ; rpa3 ; rpa4 . 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_{\rm c\to h}=335\pm 30$ K. By using a multi-stage approach (Methods, Supplementary Methods, and Supplementary Table 1) 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 which 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.

Figure 1 shows the zero-temperature energies of most popular BN polymorphs relative to that of h–BN$(AA^{\prime})$, calculated with several dispersion-corrected DFT methods and RPA (corresponding space groups and polymorph notation employed throughout the text are explained therein). These zero-temperature energies include quantum nuclear effects cazorla17 via zero-point energy (ZPE) corrections (Methods), although the ZPE effects do not lead to qualitative changes in the relative stability of BN polymorphs (Supplementary Table 2). Our results demonstrate two key points. Firstly, the RPA ordering of low-energy states confirms that c–BN is the ground state bundy62 ; solozhenko99 ; will00 , which is about $1$–$2$ kJ/mol lower in energy than the h–BN$(AA^{\prime})$, m–BN (discussed below), h–BN$(AB)$, and h–BN$(ABC)$ (also known as r–BN) polymorphs. And secondly, the many-body dispersion (MBD) method mbd agrees quite well with RPA, this being the only semi-empirical method that predicts correctly the energy ordering among all the most competitive phases. Although, the D3 method with Becke-Johnson damping [D3(BJ)] d3bj , and the fractionally ionic MBD (FI) fi method come close. Energy results obtained with the LDA and PBE methods turn out to be of unsatisfactory quality (Supplementary Table 2).

Our study reveals a low–symmetry monoclinic phase, denoted here as m–BN (space group $Cm$, different from a previously predicted monoclinic phase zhang13 hence the “2” in Fig. 1), that is energetically very competitive with respect to the h–BN polymorphs ($\Delta E^{\rm RPA}<1$ kJ/mol). This new phase is vibrationally stable and presents a reduced two–atoms unit cell with equilibrium parameters $a_{m}=4.34$, $b_{m}=2.51$, $c_{m}=3.56$ Å , and $\beta=113^{\circ}$, as obtained with MBD and FI methods (Figure 2, Supplementary Table 3, and Supplementary Figure 1). The predicted m–BN polymorph is structurally very similar to a previously reported monoclinic phase that was experimentally observed during the c–BN$\leftrightarrow$h–BN transformation occurring at high–$P$ high–$T$ conditions ($a^{\rm exp}_{m}=4.33$, $b^{\rm exp}_{m}=2.50$, $c^{\rm exp}_{m}=3.1$–$3.3$ Å , and $\beta^{\rm exp}=92$–$95^{\circ}$) and which was named as “compressed h–BN” horiuchi95 ; 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^{\prime})$, 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 horiuchi95 (Supplementary Table 4). 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–6% above their RPA counterparts. This outcome supports the use of LDA for assessing the vibrational properties of BN polymorphs accurately, a task which 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 (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 in order to provide precise phase transition data (Methods and Supplementary Figure 2). We find that the h–BN$(AA^{\prime})$ polymorph becomes entropically stabilized over c–BN at $T_{\rm c\to h}=335$ K, a temperature relatively close to ambient conditions that falls significantly below the corresponding experimental data bundy62 ; solozhenko99 ; will00 . 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_{\rm c\to h}=335\pm 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^{\prime})$ and c–BN at $T_{\rm c\to h}$ agree almost perfectly with the experimental measurements carried out at room temperature (Fig. 3) knittle89 ; solozhenko95 . The low-energy phonon excitations in the h–BN$(AA^{\prime})$ phase present much lower frequencies than in c–BN (Supplementary Figure 3), hence the vibrational entropy of the former polymorph becomes increasingly more favourable 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 follow closely that of h–BN$(AA^{\prime})$, falling all within an energy difference range of about $1$ kJ/mol at $300\leq T\leq 400$ K.

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$\leftrightarrow$h–BN transformation lies relatively close to ambient conditions. Consequently, the phase diagram of BN, if not analogous, certainly is more similar than previously thought to that of carbon. Actually, our $T_{\rm c\to 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, $420<T^{\rm exp}_{\rm c\to h}<1500$ K, may be the high-temperature high-pressure conditions and kinetic factors involved in the experimental synthesis and analysis of BN samples will00 . We hope that our theoretical results will motivate new experimental activity in bulk BN. On the theory side, we have shown: i) that LDA is a good method for calculating elastic properties of materials, and hence is probably good for estimating vibrational free energies, but that ii) for accurate prediction of energy ordering among van der Waals polymorphs one must employ methods which include many-body dispersion interactions, ideally at a high-level using RPA, but certainly approximately when RPA is infeasible.

Methods

Density functional theory and phonon calculations. First-principles calculations based on density functional theory (DFT) are performed to analyze the energy, structural, and vibrational properties of BN polymorphs. We perform these calculations with the VASP code vasp , using projector augmented-wave method potentials bloch94 . The electronic states $1s$-$2s$ of B and $2s$-$2p$ of N atoms are considered as valence. Wave functions are represented in a plane-wave basis truncated at $650$ eV. By using these parameters and dense ${\bf k}$-point grids for integration within the first Brillouin zone (IBZ), energies are converged to within $1$ meV per formula unit (0.1 kJ/mol, Supplementary Figure 4). In the geometry relaxations, a tolerance of $0.01$ eV$\cdot$Å*-1* is imposed in the atomic forces.

Ab initio phonon frequencies are calculated with the direct method in order 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 quasi-harmonic approach cazorla17 . In the direct method the force-constant matrix is calculated in real-space by considering the proportionality between atomic displacements and forces alfe09 . The quantities with respect to which our phonon calculations are converged include the size of the supercell, the size of the atomic displacements, and the numerical accuracy in the sampling of the IBZ. We find the following settings to provide quasi-harmonic free energies converged to within 0.1 kJ/mol: $4\times 4\times 4$ supercells (where the figures indicate the number of replicas of the unit cell along the corresponding lattice vectors, Supplementary Figure 4), atomic displacements of $0.02$ Å, and ${\bf q}$-point grids of $14\times 14\times 14$. The value of the phonon frequencies are obtained with the PHON code developed by Alfè alfe09 . In using this code we exploit the translational invariance of the system, to impose the three acoustic branches to be exactly zero at the center of the Brillouin zone, and apply central differences in the atomic forces.

Random phase approximation calculations. To overcome the accuracy barrier of semi-empirical 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 rpa0 ; rpa1 ; rpa2 ; rpa3 ; rpa4 . Due to its high numerical cost, our RPA calculations were performed using structures optimised at the FI level fi . Two sets of calculations were carried out. Firstly, exact exchange (EXX) and RPA correlation energy calculations,

[TABLE]

were carried out for all low-energy structures. These used an energy cutoff of 480 eV, a dense $12\times 12\times 12$ (or equivalent) ${\bf k}$-point grid for EXX, and a coarser $7\times 7\times 7$ one for RPA, evaluated on self-consistent PBE orbitals as per standard practice rpa0 . To further refine the energies of the c–BN, h–BN($AA^{\prime}$) 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 use the quasi-harmonic approach (QHA) cazorla17 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) kresse99 , hence we disregard them here. (We should note that in the unlikely case that anharmonicity played a role at $T\sim 300$ K, it would probably tend to further stabilise the hexagonal polymorph over the cubic kresse99 , thus additionally reducing $T_{\rm c\to h}$.) In the QHA approximation the vibrational free energy of a crystal $F_{\rm vib}$ with volume $V$ and at temperature $T$ is:

[TABLE]

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 $\bm{q}$ and phonon branches $s$, and $\omega_{\bm{q}s}$ are the vibrational frequencies of the crystal, which depend on volume. In the zero-temperature limit $F_{\rm vib}$ becomes:

[TABLE]

which usually is referred to as the “zero-point energy” (ZPE). The Gibbs free energy of a crystal then reads:

[TABLE]

where $E_{\rm 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:

[TABLE]

Finally, by using the thermodynamic constraint $P(V_{0},T)=0$ and performing $E_{\rm el}$ and $F_{\rm vib}$ calculations over a dense grid of volume points, it is possible to account precisely for $T$–induced volume expansion effects at zero pressure (Supplementary Figure 2). The zero-temperature energies reported in this study account for possible quantum nuclear effects by means of the expression:

[TABLE]

where “$m$” denotes the method of calculation, $V^{m}_{0}$ the resulting equilibrium volume, and $E^{m}_{\rm ZPE}$ the zero-point energy as given by Eq.(3). We have checked that the value of zero-point energy differences between BN polymorphs are practically independent of the employed method (Supplementary Figure 5), hence the reason for our fixed choice of $E_{\rm 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_{\rm 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 leconte17 , as we have explicitly corroborated in this study (Table 1). Based on numerical errors of 6% for the elastic constants of LDA, versus RPA, we assign a corresponding numerical error to the vibrational free energies of $\pm 0.06~{}$kJ/mol (calculated as 6% of the zero-point energy 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 (Supplementary Table 2).

Supplementary Materials

Supplementary material for this article is available at xxx.

table S1. Numerical tests performed for the RPA calculations.

table S2. Zero-temperature electronic energies of several BN polymorphs calculated with different first-principles methods.

table S3. Structural properties of the monoclinic phase m–BN reported in this study.

table S4. Elastic constants associated with compressive deformations calculated with the LDA method for different BN polymorphs.

fig. S1. Phonon spectrum of the new monoclinic phase m–BN calculated with the LDA method.

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 and h–BN polymorphs calculated with the LDA method.

fig. S4. Convergence tests of the electronic and vibrational free energies calculated with DFT methods for the h–BN polymorph.

fig. S5. Zero-point energy corrections for several BN polymorphs calculated with different DFT methods.

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 no competing interests. Data availability: All data needed to evaluate the conclusions in the study are present in the paper and/or the Supplementary Materials. The data that support the findings of this study are available from the corresponding authors (C.C. and T.G.) upon reasonable request.