On the thermal conductivity anisotropy in wurtzite GaN

Gallium nitride (GaN), a wide bandgap semiconductor, has recently emerged as a very promising material for high-power (HP) and high-frequency (HF) electronics.^1,2 Due to its outstanding material properties, e.g., high critical electric field, high electron mobility, and high saturation velocity, GaN-based electronic devices have potential to overperform their Si- and SiC-based counterparts. For electronic devices operating at high currents, high voltages, and high frequencies, the dissipation of the Joule heat generated in the active region becomes an important issue. Self-heating has been shown to cause current degradation and reduction of the breakdown voltage.^3–6 Hence, thermal management is essential for the device’s performance and reliability. The thermal conductivity (ThC) of GaN is a key parameter for the design and thermal management of the device structures.

The ThC of GaN has been widely studied both experimentally and theoretically. In bulk GaN, ThC in the range of 200–295 W/m K has been measured at room temperature.^7–17 Such a variation can be attributed to the difference in the dislocation density and the impurity concentration in the studied samples grown by hydride phase vapor epitaxy (HVPE), high-pressure high-temperature (HPHT) growth, or ammonothermal (AT) growth. Theoretical values for room-temperature ThC of intrinsic GaN obtained from first principles calculations are in the range of 240–340 W/m K.^18–21 

GaN crystallizes in the wurtzite crystal structure, where an anisotropy in the thermal conductivity along the two main crystallographic directions can be expected. In wurtzite crystals, the ThC anisotropy factor (ρ) can be defined as ρ = ∣2(k_‖ − k_⊥)/(k_‖ + k_⊥)∣, where k_‖ is the ThC along of c-axis and k_⊥ is the ThC in a direction perpendicular to c-axis. First-principles calculations have predicted a quite different ρ for intrinsic bulk GaN. Lindsay et al.¹⁸ have reported a very small anisotropy (ρ = 1.2%). In contrast, Ma et al.¹⁹ and Wu et al.²⁰ have calculated ρ = 13%. In another study, ρ = 5.5% has been determined.²¹ It is worth noting that the calculations of Ma et al.¹⁹ and Wu et al.²⁰ have revealed k_‖ > k_⊥, while an opposite relation (k_‖ < k_⊥) has been calculated by Lindsay et al.¹⁸ and Yuan et al.²¹ 

The experimental studies of the ThC anisotropy in bulk GaN are very limited. Zheng et al.¹⁵ have measured ThC of HVPE grown c-plane and m-plane thick GaN layers by time-domain transient thermoreflectance (TDTR). A slightly lower k_‖ than k_⊥ has been found, but the exact anisotropy could not be determined because the studied c-plane and m-plane layers were with very different impurity concentrations, which strongly affects the ThC. Li et al.²² have measured ThC of 400 μm thick GaN c-plane layers grown by HVPE and by variation of the modulation frequency and the excitation laser spot extracted k_‖ and k_⊥. They have obtained k_‖ > k_⊥ with an anisotropy factor of ρ = 4%. Despite these theoretical and experimental studies, a comprehensive understanding of the thermal conductivity anisotropy of GaN is still lacking.

In this study, we provide insights into the thermal conductivity anisotropy of GaN via combined experimental and theoretical investigations. We measure the thermal conductivity in low-defect-density and isotopically natural bulk c-plane and m-plane GaN samples in a temperature range of 80–400 K. Theoretical calculations based on a solution of the Boltzmann transport equation (BTE) within a relaxation time approximation (RTA) are also performed and the factors that determine the ThC anisotropy are analyzed.

GaN samples with a thickness of 1 mm sliced perpendicular to the c-axis (with c-plane surface orientation) and along the c-axis (with m-plane surface orientation) from boules grown by HVPE were examined. Secondary ion mass spectroscopy (SIMS) revealed a concentration of background impurities (Si, O, and C) below 5 × 10¹⁶ cm⁻³. The dislocation density in the GaN samples was estimated from x-ray diffraction (XRD) and cathodoluminescence (CL) measurements to be in the order of 1 × 10⁷ cm⁻².

Thermal conductivity measurements were performed using the transient thermoreflectance (TTR) technique. More details of the method can be found in our previous reports.^17,23,24 A gold (Au) transducer layer having a thickness of L = 200 ± 5 nm was deposited on the sample surface. The size of the probing laser is controlled to be more than one order of magnitude smaller than that of the heating laser. Therefore, the heat transport probed during the measurements can be assumed to be along the direction perpendicular to the surface. In the processing of data, a least-square minimization fitting using one-dimensional heat transport equations was performed, and the thermal conductivity of the GaN and thermal resistance at the Au/GaN interface were simultaneously obtained. In the fitting procedure, the thickness and specific heat capacity (C_p) of Au and GaN layers and the thermal conductivity of the Au were used as input parameters. The ThC and C_p of Au and the C_p of GaN were taken from the available literature data.^25–27 

The lattice thermal conductivity of semiconductors is commonly calculated by solving the BTE with the RTA.²⁸ In a simplified version of this approach, a Debye-like phonon dispersion, i.e., a linear dependence of the phonon frequency on the wavevector, is assumed.^{13,17,23,24,28–32} In this work, the GaN ThC is calculated using a non-Debye RTA model, where the exact dispersion of all phonon branches is taken into account.

Only two resistive phonon scattering processes are considered in our calculations, namely, the Umklapp three-phonon scattering (U) and the isotope-phonon scattering (I). The reason is that our GaN samples are thick bulk-like layers with low dislocation densities and low background doping, so the phonon-dislocation, phonon-point-defect, and phonon-boundary scattering have a negligible contribution to the ThC. The scattering rate of the Umklapp process is given by^31, $τU−1=BUω2(q)Te−θD/3T$ with $BU=ℏγ2(q)M−1vg−2(q)θD−1$⁠. Here, γ(q) is the Grüneisen parameter of the corresponding phonon mode, θ_D = ℏω_max(q_max)/k_B is the Debye temperature, and M is the averaged atomic mass. The scattering rate of the I-process is^31,35 $τI−1=BIΓIω4(q)$ with $BI=V/4πvg3(q)$⁠, where Γ_I is the isotope scattering parameter (Γ_I = 2.74 × 10⁻⁴)^17,31 and V is the volume per atom. The total scattering rate is obtained by using Matthiessen’s rule, i.e., $τ−1=τU−1+τI−1$⁠.

The phonon dispersion in GaN is calculated using the ab initio approach within the framework of density functional theory (DFT). A supercell approach was used with a cell dimension of (4 × 4 × 4) and (6 × 6 × 6) Monkhorst–Pack q-point mesh. The force constants are obtained from Hellmann–Feynman forces induced by atomic displacements in the supercells and calculated with Quantum Espresso package³⁶ under the local density approximation (LDA) with the use of Vanderbilt ultrasoft pseudopotentials. The energy cutoffs of 70 Ry are used to guarantee a good convergence. The phonon dispersions are then obtained from the diagonalization of dynamical matrices with the use of Phonopy code.³⁷ An analytical correction due to the dielectric constant and the Born effective charge are included.

Figure 1 shows the calculated phonon dispersion along the [0001], $[101̄0]$⁠, and $[112̄0]$ wurtzite crystallographic directions. A good agreement with the experimental data measured by inelastic x-ray scattering³⁸ is obtained. From the calculated phonon dispersion, the q_max, the Debye temperature, and the Grüneisen parameter $γ(q)=−Vω(q)δω(q)δV$ for all phonon modes were extracted. These parameters are listed in Table I. Due to the large frequency gap between the acoustic and optical phonons, the contribution of the phonons above the frequency gap in the lattice ThC is assumed to be negligible.³⁹ Then, in our calculations, only the phonons below the energy gap, namely, three longitudinal and transverse acoustic [LA, TA(z) TA(x, y)], the two optical [$E21$(low) and $E22$(low)], and the silent B₁(low) phonon modes are considered. Since our ThC measurements were done along the [0001] and $[101̄0]$ directions, the parameters only for these two directions are presented in Table I. Note that the TA(z) and TA(x, y), as well as $E21$(low) and $E22$(low) modes along the [0001] direction, are degenerate. In Table I, γ_cal is the Grüneisen parameter averaged over the phonon wavevector, i.e., γ_cal = $⟨γ2(q)⟩$⁠.³¹ 

Figure 2 shows the measured ThC along the [0001] and $[101̄0]$ crystallographic directions in the temperature range 80–400 K together with the ThC temperature dependence calculated by Eq. (3). At room temperature, we have measured k_‖ = (225 ± 11) W/m K and k_⊥ = (215 ± 11) W/m K resulting in an anisotropy factor of ρ = 4.5%. This value is very close to that obtained by Li et al.²² It is worth noting that thermal conductance at Au/GaN interface is ∼100 MW/m⁻² K⁻¹ at a measurement temperature of 300 K, and this value is comparable with 56.7–256.6 MW/m⁻² K⁻¹ as reported in existing literature.⁴⁰ In addition, the ThC of GaN and thermal resistance at Au/GaN interface are individually analyzed; the ThC obtained from this analysis, therefore, remains unaffected by the properties of the interface, thus ensuring its reliability and accuracy. There is a very good agreement between the calculated temperature dependence of ThC and the experimental data. Our calculations reveal k_‖ < k_⊥ at room temperature and an increasing anisotropy at low temperature in accordance with the first principle calculations of Lindsay et al.¹⁸ and Yuan et al.²¹ The estimated relative contributions of different phonon modes to the ThC are listed in Table I. It is seen that the thermal conductivity of GaN is entirely dominated by the scattering of acoustic phonons. A comparison with previously published results on the GaN ThC along the [0001] and $[101̄0]$ directions^15,22 reveal that our experimental data is generally higher (Fig. 2). At lower temperatures, these differences become more significant. Nevertheless, the obtained anisotropy factor ρ at room temperature is consistent with our findings. The lower ThC measured in Refs. 15 and 22 can be attributed to the higher concentration of background impurities in the studied samples, which leads to an enhanced phonon-point-defect scattering. In Fig. 2, we also show the results from the previous first principle calculations of GaN ThC in the directions parallel and perpendicular to the c-axis.^18,19 The calculations of Lindsay et al.¹⁸ reveal a negligible anisotropy factor ρ in the entire temperature range studied. In contrast, the results of Ma et al.¹⁹ show a significant ThC anisotropy (ρ = 13% at room temperature). Note that the calculations in Ref. 19 were done for [0001] and $[112̄0]$ crystallographic directions. Despite the similar phonon dispersion along $[101̄0]$ and $[112̄0]$ directions, a difference between the ThC in both directions may be expected.

The small ThC anisotropy observed experimentally is quite surprising having in mind the high crystal anisotropy of wurtzite GaN (c-to a-axis lattice constant ratio of 1.63).⁴¹ To get further insight into the ThC anisotropy in GaN, we have made an analysis based on Eq. (3). The larger extension of the first Brillouin zone (in q-space) in directions perpendicular to the c-axis compared to that along the c-axis implies higher q_max and higher mode Debye temperatures (see Table I). As a consequence, k_‖ < k_⊥ is expected. On the other hand, the higher phonon group velocity along [0001] direction (averaged over the six modes in consideration) suggests k_‖ > k_⊥. The strength of Umklapp three-phonon scattering rate is weighted by the Grüneisen parameters, which measure the crystal anharmonicity. The higher the mode Grüneisen parameter, the shorter the scattering time. We have found that the averaged Grüneisen parameter over TA, LA, E₂(low), and B₁(low) phonon modes along the $[101̄0]$ direction is by about 10% higher than that for [0001] direction. Thus, k_‖ > k_⊥ can be predicted from Eq. (3). Obviously, the ThC anisotropy is determined by the interplay between the crystal anisotropy (i.e., the anisotropy of the first Brillouin zone) and the anisotropy of the group velocity and Grüneisen parameter of different phonon modes. As a result, the contributions of these effects cancel out, which can explain the small ThC anisotropy in GaN.

To get insight into the thermal transport in undoped bulk GaN, we have calculated the cumulative thermal conductivity as a function of the phonon mean free path (MFP). The normalized cumulative ThC along the [0001] and $[101̄0]$ crystallographic directions is shown in Fig. 3. We have found an MFP value at 50% of the cumulative ThC (MFP₅₀) of 1.2 and 2.0 μm for the $[101̄0]$ and [0001] directions, respectively, implying a higher averaged MFP for the phonons along the [0001] direction. The cumulative ThC along the [0001] calculated by the Debye–Callaway model in which only the three acoustic phonons are considered²⁴ starts at a rather high MFP and has an abrupt increase with increasing the MFP. This result is explained by the overestimated contribution of the phonons at high frequencies (i.e., with a low MFP), which is due to the assumed linear phonon dispersion and the constant sound velocity in the Debye-Callaway model. The available experimental data for the cumulative k_‖⁴⁵ are also presented in Fig. 3. Note that the slope of the curve matches much better the calculations with our non-Debye RTA model than those by the Debye–Callaway model. Surprisingly, the first principle calculations⁴⁶ reveal a much lower MFP₅₀ (about 400 nm for both the [0001] and the $[101̄0]$ directions) and a lower layer thickness beyond which the size effect becomes negligible, thus predicting a weaker size effect of ThC.

In conclusion, we have revisited the anisotropy of the thermal conductivity of GaN. We have determined the ThC of c- and m-plane oriented GaN bulk samples at temperatures from 80 to 400 K. The experimental results reveal a small ThC anisotropy in the entire temperature range studied. This finding is in very good agreement with the calculations based on our non-Debye RTA model for lattice thermal conductivity. The small ThC anisotropy is explained by the interplay between the group velocity, Debye temperature, and Grüneisen parameter anisotropy.

The authors are indebted to Dr. J. H. Leach from Kyma Technologies, Inc., for providing the HVPE samples for this study. This work was performed within the competence centers for III-Nitride Technology (C3NiT-Janzén) supported by the Swedish Governmental Agency for innovation systems (VINOVA) under the Competence Center Program Grant No. 2022-03139. We further acknowledge support from the Swedish Research Council VR (Grant Nos. 2016-00889, 2017-03714, and 2022-04812), the Swedish Foundation for Strategic Research (Grant Nos. RIF14-055 and EM16-0024), and the Swedish Government Strategic Research Area in Materials Science on Functional Materials at Linköping University, Faculty Grant SFO Mat LiU Nos. CBET-1336464 and DMR-1506159. The work at NCSU was supported by NSF (Grant Nos. CBET-1336464 and DMR-1506159). The first-principles computations were enabled by resources provided by the Swedish National Infrastructure of Computing (SNIC).

The authors have no conflicts to disclose.

Dat Q. Tran: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Writing – original draft (lead); Writing – review & editing (equal). Tania Paskova: Resources (equal); Validation (equal); Writing – review & editing (equal). Vanya Darakchieva: Funding acquisition (lead); Project administration (lead); Resources (equal); Validation (equal); Writing – review & editing (equal). Plamen P. Paskov: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Supervision (lead); Validation (equal); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.