Phonons of hexagonal BN under pressure: Effects of isotopic composition

Raman scattering experiments on isotopically enriched hexagonal boron nitride have been performed under pressure up to 11 GPa at room temperature. The sublinear increase of the Raman-active ${E}_{2g}$ mode frequencies has been characterized. The pressure behavior has been analyzed by means of a bond-stiffness--bond-length scaling parameter $\ensuremath{\gamma}$ which takes into consideration the vast differences in $a$- and $c$-axis compressibilities. The interlayer shear mode exhibits a $\ensuremath{\gamma}$ parameter similar to that of graphite, and the mode frequency in isotopically pure samples separates faster at low pressures as a result of van der Waals interactions. Because of the extremely low $a$-axis compressibility, the intralayer mode exhibits a rather large scaling parameter. Within experimental uncertainties, no systematic departures related to isotopic mass have been observed.


I. INTRODUCTION
Hexagonal boron nitride (h-BN) is a layered material with a prominent role in present-day electronic and photonic applications based on graphene and related two-dimensional (2D) materials [1] owing to its unique properties as a dielectric insulator, strong ultraviolet emitter [2,3], and natural optical hyperbolic material [4]. Bulk cubic boron nitride (c-BN) and h-BN are good thermal conductors and therefore promising candidates for thermal management of electronic and optoelectronic applications involving high power, high temperature, or high photon energy. h-BN has the potential for providing an insulating substrate with high thermal conductivity for flexible electronic devices and novel two-dimensional van der Waals structures [5]. In fact, thermal conductivity has been reported to increase in high-quality h-BN monolayers due to the reduction in phonon states available for umklapp scattering [6].
Another route to increase thermal conductivity is isotope engineering, which allows the reduction of phonon scattering by isotopic disorder to be achieved in isotopically purified crystals [7,8]. Very recently, an ∼90% enhancement of the room-temperature thermal conductivity to ultrahigh values over 1600 W m −1 K −1 was reported in isotopically enriched c-10 BN [9]. In-plane thermal conductivity as high as 585 W m −1 K −1 was measured at room temperature in monoisotopic 10 BN [10]. A substantial increase of the thermal conductivity up to 1009 W m −1 K −1 was also observed in atomically thin monoisotopic h-BN owing to the suppression of the out-of-plane optical phonon scattering as well as of the interlayer interactions leading to umklapp scattering of the out-of-plane acoustic phonons [11]. Besides improving the thermal properties, isotope engineering has profound implications for other application areas of h-BN. 10 B is one of the best neutron absorbers, and devices based on 10 B-enriched h-BN have demonstrated the highest detection efficiency for solid-state neutron detectors to date [12]. On the other hand, the suppression of phonon scattering by isotopic disorder and hence the substantial increase of phonon lifetime in isotopically purified samples [13] have opened up new avenues for phonon-polariton applications of h-BN [14]. Hyperlensing in the infrared [15] and high-quality-factor resonators in subdiffractional cavities [4] have been demonstrated in low-loss isotopically purified h-BN [16].
The availability of high-quality single-crystal h-BN samples [17] and the precise control of the isotopic composition [18] have enabled us to extensively address the optical and vibrational properties of bulk h-BN in a number of recent works [13,[19][20][21][22][23][24][25][26][27]. Raman scattering has provided detailed information about the intralayer and interlayer phonon modes and their lifetimes [20,21], and recently, phonon behavior under pressure was studied in natural-composition h-BN up to the phase transition to the wurtzite phase [26]. Raman scattering experiments on isotopically enriched crystals under pressure can yield detailed information about the interplay between mechanisms determining phonon lifetime such as anharmonicity and isotopic disorder. Indeed, large variations of phonon lifetimes as a function of pressure and isotopic composition were reported in zinc-blende and wurtzite semiconductors, namely, CuI [28], ZnS [29], and ZnO [30]. Similar studies on layered materials are lacking in the literature.
In this work, we present Raman scattering measurements of the intralayer and interlayer Raman-active modes under pressure of isotopically enriched samples over the whole range of isotopic compositions. The phonon linewidth of the high-frequency intralayer mode does not show significant variations across the pressure range studied, and the main differences among the different isotopic compositions are related to the degree of isotopic disorder. The low-frequency interlayer shear mode is most sensitive to the weak interlayer van der Waals interactions. Differences in the pressure dependence of the E low 2g frequency have been revealed in isotopically pure samples. The change in phonon frequency with pressure is modeled by using density-functional perturbation theory (DFPT) calculations, which provide insight into the sublinear pressure dependence of the phonon frequencies.

II. EXPERIMENT
The isotopically enriched samples were grown from highpurity elemental 10 B (99.2 at %) or 11 B (99.4 at %) powders. A Ni-Cr-B powder mixture at 48, 48, and 4 wt %, respectively, was loaded into an alumina crucible and placed in a furnace. h-BN with the natural distribution of boron isotopes was made with 50 wt % Ni and 50 wt % Cr mixtures with a hot-pressed ceramic h-BN crucible, which also served as the source material. The furnace was evacuated and then filled with N 2 and forming gas (5% hydrogen in balance argon) to a constant pressure of 850 Torr. During the reaction process, the N 2 and forming gases continuously flowed through the system with rates of 125 and 25 sccm, respectively. With the pure boron sources, all the nitrogen in the h-BN crystal originated from the flowing N 2 gas. The forming gas was used to minimize oxygen and carbon impurities in the h-BN crystal. After a dwell time of 24 h at 1550 • C, the h-BN crystals were precipitated by cooling at a rate of 1 • C/h to 1500 • C, and then the system was quickly quenched to room temperature [18]. Samples with thickness ranging from 30 to 100 μm were cleaved from h-BN single-crystalline platelets.
The high quality of these h-BN crystals has been verified by several techniques. The crystals are highly oriented with the [0001] axis perpendicular to the sample surface, as determined by x-ray diffraction (XRD) [31]. The rms surface roughness over a 2 × 2 μm 2 area was 0.4 nm, as measured by atomic force microscopy. X-ray photoelectron spectroscopy (XPS) detected only boron and nitrogen, and these peaks shift to slightly higher energies as the 11 B isotope concentration increases from 0% to 100%. [18] The low-temperature photoluminescence spectra have multiple narrow phonon replicas above 5.75 eV, which are present only in high-quality h-BN [23]. The etch pit density was 5 × 10 7 cm −2 on h-BN crystals grown by this technique, the pits forming where dislocations (identified by transmission electron microscopy) intersect the surface upon defect-sensitive etching [32].
The high-pressure Raman scattering spectra were measured up to the phase-transition pressure (∼11 GPa) [25,26] at room temperature in a membrane diamond anvil cell (DAC) [33] equipped with type-IIac diamond anvils. Ne was used as a (quasi)hydrostatic pressure-transmitting medium. Pressure was determined using the ruby scale [34]. The Raman spectra were recorded in the backscattering configuration from the c face using the 514.5-nm line of an Ar + laser and a confocal microscope setup with a 50× objective. The scattered light was analyzed using a Jobin-Yvon T64000 triple spectrometer equipped with a LN 2 -cooled charge-coupled device detector. The spectral bandwidth of the experimental setup was ∼2 cm −1 .
DFPT calculations were performed using the ABINIT code [35]. The calculations were carried out in the generalized gradient approximation (GGA). Nonlocal van der Waals interactions between layers [36] were taken into account by considering the density-functional dispersion correction DFT-D3(BJ), as proposed by Grimme et al. [37] and Becke and Johnson [38]. Following the recommendations in Ref. [36], Zhang and Wang's revPBE pseudopotential [39] was used. Contrary to plain GGA, which does not bind the hexagonal layers, the incorporation of nonlocal corrections to the functional yields a good description of the layered h-BN crystal [36]. Figure 1 shows representative spectra of isotopically enriched samples ( 10 B x 11 B 1−x N, with x = 0, 0.5, 1) at three different pressures. The Raman spectra exhibit two peaks, corresponding to the Raman-active modes of E 2g symmetry. The low-frequency mode is due to an interlayer shear mode which involves a gliding motion of the rigid hexagonal layers, whereas the high-frequency mode involves in-plane boronnitrogen vibrations. Both modes display similar blueshifts with pressure, despite their vast difference in frequency. This is a consequence of the high anisotropy of the h-BN layered crystal, where the compressibility along the c direction is about 60 times higher than along the in-plane directions [40]. Owing to the low frequency of the E low 2g mode, there is a lack of efficient channels for anharmonic decay into lowerenergy modes, and this mode exhibits a long lifetime [21,22]. Therefore, the linewidth of its Raman spectrum corresponds to the slit width of the Raman spectrometer, and it is essentially the same for all the samples and pressures studied. The peak intensity increases with pressure for all samples, suggesting a stronger interlayer coupling and an increase of mode polarizability with increasing pressure [26]. The linewidth of the high-frequency mode E high 2g changes very little with pressure. This is in contrast to the case of ZnO, where the E high 2g frequency lies on a ridgelike structure of the two-phonon density of states, and the effects of anharmonic decay are strongly modulated by isotopic composition and pressure effects [30]. which does not significantly modulate the anharmonic decay rate [19,20]. As a result, the E high 2g line shape is symmetric in isotopically pure h-BN crystals and its linewidth does not change significantly with pressure, in contrast to the large variations of the E high 2g linewidth observed in ZnO [30]. In the isotopically mixed sample (50% 10 BN /50% 11 BN), the E high 2g linewidth is noticeably larger but does not show significant variations with pressure either. Again, the region of the phonon density of states around the E high 2g frequency is slowly varying for all isotopic compositions [19,20], and therefore, pressure does not produce a significant modulation of the isotopic disorder scattering rate. Figure 2 shows the Raman shifts of the two E 2g modes of h-BN as a function of pressure for different isotopic compositions up to the phase-transition pressure to the wurtzite phase [25,26]. For completeness, we have also included the data on natural-composition h-BN reported in a previous publication [26]. Because of the high structural anisotropy of h-BN, the compressibility along the c direction is large and nonlinear with pressure, whereas the in-plane lattice parameter shows only a very small linear decrease with pressure. The nonlinearity in the c lattice constant results in a sublinear increase of the phonon frequencies that can be described using a Murnaghantype dependence [26]

III. RESULTS AND DISCUSSION
Here, δ 0 is the normalized linear pressure coefficient at ambient pressure, and δ is the sublinearity parameter. The dashed lines in Fig. 2  guide to the eye since the frequencies of the different isotopes lie very close together. A more detailed analysis of the E low 2g frequency variation among the isotopically purified samples will be given later (see discussion of Fig. 4). The values of the fitted δ 0 and δ parameters are listed in Table I composition, and the small variations can be considered to be within experimental uncertainty.
The results of ab initio calculations are displayed in Fig. 2 as solid lines. A correction for the frequency shift of the E high 2g mode due to isotopic disorder was introduced for the 10 B0.5 11 B0.5 N ( E high 2g = 4.6 cm −1 ) and natural-composition samples ( E high 2g = 3.7 cm −1 ) according to the results reported in Ref. [13]. After that correction, the curves for all compositions were slightly above the experimental points. Density-functional theory calculations in the local-density approximation are known to overbind and thus overestimate phonon frequencies [41]. In our case, the GGA approach with the DFT-D3(BJ) correction ameliorates this limitation, but some overbinding still remains, which leads to a marginal overestimation of the phonon frequencies of just E high 2g ∼ 4 cm −1 . By correcting all the calculated frequencies by this amount, we obtain very good agreement with the experimental pressure dependence of the E high 2g frequency for all isotopic compositions. Note that the overestimation of the ab initio calculations cancels the isotopic disorder shift for the case of the natural-composition sample. This explains why the GGA calculations with the DFT-D3(BJ) correction provided a remarkably accurate description of the pressure dependence of the E high 2g mode in a previous study on a natural-composition sample [26].
Vibrational frequencies increase under compression because of anharmonicity effects that stiffen the bonds as they shorten. In the h-BN layered crystal under pressure, the compression of the interlayer distance is much higher than the shortening of the in-plane boron-nitrogen distances. Thus, the interlayer force constant governing the E low 2g mode is much more affected by pressure than the intralayer force constant. To describe this situation, following Zallen [42], we assume that each bond type i individually obeys a scaling law where r i is the corresponding bond length (r 0 = a and r 1 = c) and γ is a scaling parameter which should be a constant of order unity for each bond type. The scaling law (2) implies a pressure dependence of the mode frequency ( For an isotropic crystal, the scaling parameter γ coincides with the macroscopic Grüneisen parameter γ (i) G = −d (ln ω i )/d (ln V ). However, in the case of the anisotropic h-BN system, the macroscopic Grüneisen parameters of the intralayer and interlayer modes vastly differ, as γ intralayer G is affected by the ratio of linear compressibilities along the a and c directions [42]. According to Solozhenko and coworkers' data [40], the a-axis compressibility of h-BN is extremely small. A linear fit a = a 0 (1 + β −1 0 P) to their experimental data yields β 0 = 2130 GPa, which is nearly twice the value reported for graphite (β 0 = 1250 GPa) by Hanfland et al. [43]. We have corroborated the small a-axis compressibility of h-BN by carrying out XRD measurements under pressure on our samples. Within experimental error, our results are fully consistent with those reported in Ref. [40]. Remarkably, the pressure coefficient of the E high 2g mode is higher in h-BN than in graphite (5.3 vs 4.7 GPa, respectively) [26,43], which suggests a more pronounced stiffening of the boron-nitrogen bond with pressure. This is probably due to the sizable ionicity of the boron-nitrogen bond [44] Table I.
Within experimental uncertainties, the scaling parameter γ for the E low 2g interlayer mode is about 1.56 for all samples. This value is only slightly higher than the value of 1.4(1) reported for the corresponding mode of graphite [43], indicating that the interlayer interactions in h-BN and graphite are very similar in nature [45]. As in the case of graphite, the larger value of γ in h-BN suggests a strong anharmonicity of the E low 2g mode. Indeed, in h-BN, anharmonic effects were found to induce significant relative shifts of the E low 2g mode frequency with temperature [19,21].
The scaling parameter for the E high 2g intralayer mode is, within experimental uncertainties, about 2.3-2.4 and does not show any clear trend with isotopic composition either. Contrary to the case of graphite, the scaling parameter of the intralayer mode is sizably higher than that of the interlayer mode. In cubic and zinc-blende semiconductors, the Grüneisen parameter γ G spans values from about 1 to about 2, and a correlation between γ G and ionicity was established [46]. The ionicity of the boron-nitrogen bond could partly explain the higher γ value of the E high 2g mode. Hanfland et al. [43] compared the scaling factor γ found for the intralayer mode of graphite with the Grüneisen parameter γ G of diamond. Despite the different natures of the carbon-carbon bonds in graphite (sp 2 ) and diamond (sp 3 ), they found that the scaling parameter in graphite was comparable to the Grüneisen parameter of tetrahedral diamond. In the case of h-BN, we find from our analysis a γ value that is significantly higher than accepted values of γ G for the TO mode of c-BN [47,48]. The differences between sp 2 and sp 3 bonding could partly explain this observation. The γ values we find for bulk h-BN crystals are slightly higher than the Grüneisen parameters determined in h-BN flakes under uniaxial or biaxial strain [49,50]. However, a direct comparison of the three-dimensional Grüneisen parameter determined by measurements under hydrostatic compression (defined by relative volume changes) with the parameters obtained in 2D-like flakes under uniaxial or biaxial traction (defined by relative length or surface changes) is not straightforward. On the other hand, the high values of γ that we find for the E high 2g mode are linked to the exceedingly small variation of the a lattice parameter with pressure and thus could be affected by experimental uncertainties in the a parameter. Figure 3 shows the pressure coefficients at ambient pressure of the two E 2g modes of h-BN for different isotopic compositions. No definite trend is observed in the experimental points, which fluctuate around 8 cm −1 GPa −1 for the low-frequency mode and around 5 cm −1 GPa for the highfrequency mode. Assuming that, as expected, the isotopic mass does not have any significant effect on the change in the force constants under pressure, the mass effect on the pressure coefficient for the intralayer mode would be on the order of μ −3/2 δμ ∼ 2%, which is too small to observe. The slightly decreasing trend with isotopic mass in the ab initio calculations (dashed lines) corresponds to this effect. The calculated pressure coefficients are in reasonable agreement with the experiment. As previously discussed [26], the nonlocal correction introduced in the calculations tends to overestimate the van der Waals interactions and leads to a smaller interlayer spacing than found experimentally. As the interlayer distances are further compressed under pressure, the contribution of Pauli repulsion becomes stronger, and further compression is more difficult, resulting in a calculated E low 2g pressure coefficient that is smaller than measured in the experiments.
The calculated pressure dependence of the E low 2g frequency plotted in Fig. 2 shows some dispersion with the isotopic mass. To demonstrate more clearly this dispersion, in Fig. 4(a) we plot the E low 2g frequencies after subtracting a common linear trend obtained as a straight line between the average of the frequencies at ambient pressure and the average of the frequencies at 10.4 GPa. Note the nonlinear bowing of the frequencies as a function of pressure, which rapidly separates the E low 2g frequencies of the different isotopic compositions. Since the interlayer shear mode is most sensitive to the weak van der Waals interactions, it is interesting to study the evolution with pressure of the difference in E low 2g frequencies between the two isotopically pure samples. In Fig. 4(b) we plot the difference between the E low 2g frequencies of 10 BN and 11 BN samples as a function of pressure. Note that both samples were measured in the same DAC run so that pressure conditions were the same. At low pressures, the E low 2g frequencies tend to separate faster, whereas at higher pressures the frequency difference tends to saturate to a value around 2.2 cm −1 . This behavior is well described by the ab initio calculations that include the densityfunctional dispersion correction DFT-D3(BJ) to account for van der Waals interactions, which are plotted in Fig. 4(b) as a solid line. In contrast, when the nonlocal correction is not included, the calculations predict a uniform increase of the frequency difference [dotted line in Fig. 4(b)]. This suggests that, initially, the applied pressure modifies the van der Waals forces differently in the isotopically pure samples. The separation of the E low 2g frequencies between the 10 BN and 11 BN samples increases by up to ∼1.5 cm −1 under pressure, which 085204-5 is a large relative variation given the low frequency of the mode. The higher pressure coefficient of the 10 BN sample could be related to a more diffuse electronic distribution between adjacent layers reported in monoisotopic 10 BN [23]. At higher pressures, the Pauli repulsion takes over, and the frequency behavior assumes a more covalent character. While the experimental points of the frequency vs pressure curves tend to separate at low pressures, at higher pressures the ω(P) curves for both isotopes run parallel.

IV. CONCLUSIONS
We have carried out a detailed study of both E 2g Ramanactive modes of h-BN under high pressure on a set of isotopically purified samples spanning the whole isotopic composition range. Frequency shifts due to isotopic disorder in mixed composition samples are relevant to account for the observed E high 2g frequencies. There is good agreement between DFPT calculations in the GGA approximation with nonlocal corrections and the observed pressure dependence of the E 2g modes. The use of the bond-stiffness-bond-length scaling parameter γ allows the analysis of the pressure dependence of the intralayer and interlayer modes in a more transparent way than conventional Grüneisen parameters, as it avoids the distorting effects of the vast differences between a-and c-axis compressibilities.
Within experimental uncertainties, the E low 2g interlayer shear mode exhibits a scaling parameter γ ∼ 1.56, without significant differences among the isotopic compositions studied. The γ value we find for h-BN is comparable to the value reported for graphite, indicating similar interlayer van der Waals interactions in both compounds. On the other hand, the E high 2g intralayer mode exhibits a scaling parameter around 2.3-2.4. This rather high γ value stems from the extremely low a-axis compressibility of h-BN, as reported from x-ray measurements with synchrotron radiation. This result indicates that the E high 2g mode exhibits strong anharmonicity owing to the sp 2 bonding and the substantial ionic character of the in-plane boron-nitrogen bond.
A signature of the van der Waals interactions was observed by comparing the pressure dependence of the E low 2g interlayer shear mode in isotopically pure samples. At low pressures, where van der Waals forces dominate the interlayer interactions, the mode frequencies diverge faster than at high pressures, where Pauli repulsion sets in. This behavior was confirmed by DFPT calculations with nonlocal corrections.
On the whole, the pressure behavior of the Raman modes in isotopically enriched h-BN samples is consistent and well understood, and no significant departures of the bond-stiffness-bond-length scaling parameters associated with isotopic mass have been detected within the experimental uncertainties.