Polymorphism of boron phosphide: Theoretical and experimental assessments

Stable crystal structures of wurtzite ( w -BP) and recently discovered rhombohedral ( rh -BP) polymorphic modifications of boron phosphide were obtained based on crystal chemistry rationale and unconstrained geometry optimization calculations within the density functional theory (DFT), and compared with known cubic polymorph ( c -BP). Both w -BP and rh -BP are mechanically (elastic constants) and dynamically (phonons) stable and exhibit thermodynamic and mechanical properties very close to those of c -BP. The electronic band structures depict semi-conducting behavior with band gap magnitudes close to 0.5 eV for cubic and rhombohedral polymorphs and 1 eV for w -BP.


Introduction
Boron phosphide BP is refractory wide bandgap A III B V semiconductor [1] characterized by a unique combination of mechanical [2,3], thermal and electrical properties [4][5][6], high thermoelectric power [7] and outstanding chemical and high-temperature stability that makes it a promising material for a wide range of engineering applications [8].
Cubic sphalerite (F-43m) BP (c-BP) [9] is the thermodynamically stable phase of boron phosphide at ambient conditions, and it has been found to be stable at pressures up to 110 GPa [10].
Hexagonal wurtzite (P6 3 mc) BP polymorph (w-BP) has also been reported [11], however, data on structure, stability region and properties of this phase are lacking in the literature. Very recently, the third polymorph, rhombohedral BP (rh-BP), has been discovered by transmission electron microscopy [12].
In the present work structure and properties of rhombohedral and wurtzite BP polymorphs have been studied using crystal chemistry rationale and unconstrained geometry optimization calculations within the density functional theory (DFT).

Computational framework
The search for the ground state structures and energies as well as energy dependent quantities within DFT was carried out using plane-wave Vienna Ab initio Simulation Package (VASP) code [13,14] with the projector augmented wave (PAW) method [14,15] for the atomic potentials. The exchange-correlation XC effects within DFT were considered using the generalized gradient approximation (GGA) [16]. For the methodology details, the reader is referred to our previous papers [17,18].
For each BP polymorph, the mechanical properties were calculated (i) from the set of elastic constants C ij , and (ii) using three contemporary theoretical models of hardness [19][20][21]. To verify the dynamic stability of BP phases, the calculations of the phonon bands have been performed. The approach consists of computing the phonon modes through finite displacements of the atoms off their equilibrium positions to subsequently obtain the forces from the summation over the different configurations. The phonon dispersion curves along the main lines of the respective Brillouin zone (BZ) were obtained using first-principles phonon calculation code "Phonopy" [22], and thermodynamic properties (Helmholtz free energy, heat capacity and entropy) were calculated as functions of temperature within the harmonic approximation. Finally, the electronic band structures showing small gap insulating behavior were assessed.

Energy and crystal structures
Starting from crystal structures of c-BP (ICSD No. 29050) and w-BP (ICSD No. 615155), the calculations were carried out for the ground state configurations. In the case of rh-BP, the structural data of rhombohedral ZnS [23] were used.
Energy is the prevailing criterion for a first assessment of the crystal structures. From the total electronic energies (E tot ) the cohesive energies were obtained by deducting the atomic contributions. The resulting atom cohesive energies of three BP polymorphs were then compared. The results provided in Table 1 show close atom average cohesive energies with a slightly more stable rhombohedral form. This validates the choice of rh-ZnS as the starting model for the structural refinement of rh-BP.
The resulting structures of three BP polymorphs shown in Figs. 1(a-c) exhibit tetrahedral stacking like that of diamond and dense polymorphs of BN. It should be noted that the structure of rhombohedral BP can be considered as three-dimensional (3D) form of a 2D analogue with 3 B-P layers shown in Fig. 1d. The similar 2D structure has been considered earlier for the study of rh-B 2 N 2 , a new rhombohedral polymorph of boron nitride [17]. However, the two-dimensional rh-BP is found less cohesive than the 3D rh-BP as shown in the last line of Table 1.
The crystal data provided in Table 2  Comparison of simulated X-ray diffraction patterns of all BP polymorphs ( Fig. 2) reveals that rhombohedral BP is rather difficult to distinguish from cubic BP by powder X-ray diffraction.

Charge density projections
Further illustration of electronic / crystal structure relationship focusing on the tetrahedral stacking can be provided by the charge density projections. Figs. 3(a-c) show for the three polymorphs the charge density yellow volume within the BP 3 and PB 3 tetrahedra, skewed toward P as expected from the higher electronegativity of P (χ =2.19) versus B (χ =2.03). Due to high electronegativity of N (χ = 3.04) in the case of polar covalent B-N bonding [17], the character of B-P bond is expected to be more covalent. For the sake of comparison, Fig. 3d presents the sp 2 -like planar charge density in 2D rh-BP (Fig. 1d).

Mechanical properties (i) Elastic constants
The elastic constants C ij needed to derive the mechanical properties were determined by performing finite distortions of the lattice and deriving C ij from the strain-stress relationship. Then fully the bulk B and shear G moduli are obtained by averaging the single-crystal elastic constants using, here, Voigt's method [25] based on a uniform strain. The calculated sets of elastic constants are given in Table 3. All values are positive. Their combinations obeying the rules pertaining to the mechanical stability of the phase, and the equations providing the bulk B V and shear G V moduli are as follows. • For the hexagonal (trigonal) system: ℎ . = 1/9 {2(C 11 + C 12 ) + 4C 13 + C 33 } G Voigt hex. = 1/30 {C 11 +C 12 + 2C 33 − 4C 13 + 12C 44 + 6(C 11 − C 12 )} As can be seen from Table 3, all three BP polymorphs exhibit the same bulk moduli while their shear moduli differ markedly, indicating different brittleness within the series with a clear tendency of decreasing brittleness from c-BP to rh-BP.

(ii) Hardness
Vickers hardness (H V ) was predicted using three contemporary theoretical models of hardness [19][20][21]. The thermodynamic model [19] is based on thermodynamic properties and crystal structure, Lyakhov-Oganov approach [20] Table 4. Table 5 summarizes hardness calculated using tree theoretical models and other mechanical properties i.e. shear moduli (G), Young's moduli (E), the Poisson's ratios (ν) and fracture toughness (K Ic ). Regardless of the model used, all three BP polymorphs are characterized by very close hardness values. As for the absolute H V values, the thermodynamic model seems to be the most reliable, because it shows the best agreement with experimental data for c-BP [2]. It is worth mentioning that all three polymorphs have virtually equal bulk moduli, while their shear and Young's moduli differ noticeably, decreasing in the c-BP -w-BP -rh-BP row (see Table 5). Fracture toughness of all three polymorphs is the same and almost twice lower than that of single-crystal cubic BN [26].

Dynamical stabilities from the phonons
Further stability criteria can be sought from the phonon dispersion relations in the Brillouin zone, i.e. the phonon band structures. Following the method presented in the 'Computational framework' section, the obtained phonon band structures for three BP polymorphs are shown in Fig. 4. In each panel the bands run along the main lines of the cubic BZ (Fig. 4a) and hexagonal BZ (Figs. 4b and   4c).
Along the vertical direction the frequency is given in units of Terahertz (THz). Since no negative energy magnitudes are observed, the three BP structures can be considered as dynamically stable.
There are 3N-3 optical modes at high energy and 3 acoustic modes. As mentioned in Table 3, the novel rhombohedral phase rh-BP is expressed in hexagonal coordinates, and BZ directions are the same as in w-BP (Fig. 4b).
The acoustic modes start from zero energy (ω = 0) at the Γ point, center of the Brillouin Zone, up to a few Terahertz. They correspond to the lattice rigid translation modes of the crystal (two transverse and one longitudinal). But the number of explicit bands is subjected to the lattice symmetry, i.e., the higher the symmetry, the more dispersion curves showing degeneracy are found for a given frequency, and therefore an apparent reduction in the number of dispersion curves is observed. In all three panels the energy range is the same for the three phases, i.e. from 0 to 25 THz, stressing furthermore their similitude but with a much lower magnitude than observed for diamond by Raman spectroscopy: ω ~40 THz [27].

Thermodynamic properties
Thermodynamic properties of BP polymorphs have been calculated from the phonon frequencies using the statistical thermodynamic expressions on a high precision sampling mesh in the Brillouin zone. The properties of all three polymorphs were found to be very similar over a wide temperature range. Fig. 5 shows the temperature changes of heat capacity at constant volume C V , entropy S and Helmholtz free energy F = U -TS, where U is the internal energy, S is the entropy and T is the absolute temperature.
For the heat capacity of cubic BP, an excellent agreement of calculated C V values with experimental data [28][29][30] is observed in the 5-700 K range. The deviation of experimental points from the calculated curve above 700 K is evidently due to the limitations of the method used (AC calorimetry [28]) at high temperatures.

Electronic structures
The electronic band structures are presented in Fig. 6 Table 4 Vickers hardness (H V ) and bulk moduli (B 0 ) of boron phosphide polymorphs calculated in the framework of the thermodynamic model of hardness [19] Space group