Thermal conductivity for III-V and II-VI semiconductor wurtzite and zinc-blende polytypes: the role of anharmonicity and phase space

TL;DR

This study uses ab initio calculations to analyze and predict the thermal conductivity differences between cubic and hexagonal phases of semiconductors, considering anharmonicity and phase space effects, and explores tuning these properties in nanowires.

Contribution

It introduces a predictive method for determining the more conductive phase in semiconductors based on anharmonicity and phase space, advancing understanding of thermal transport in different crystal structures.

Findings

Cubic phase exhibits higher anharmonicity than hexagonal.

Accessible phase space favors hexagonal phase conductivity.

Nanowire diameter can tune the ratio of thermal conductivities.

Abstract

We calculate the lattice thermal conductivity ($\kappa$) for cubic (zinc-blende) and hexagonal (wurtzite) phases for 8 semiconductors using $\textit{ab initio}$ calculations and solving the Phonon Boltzmann Transport Equation, explaining the different behavior of the ratio $\kappa_{\rm hex}/\kappa_{\rm cub}$ between the two phases. We show that this behavior depends on the relative importance of two antagonistic factors: anharmonicity, which we find to be always higher in the cubic phase; and the accessible phase space, which is higher for the less symmetric hexagonal phase. Based on that, we develop a method that predicts the most conducting phase---cubic or hexagonal---where other more heuristic approaches fail. We also present results for nanowires made of the same materials, showing the possibility to tune $\kappa_{\rm hex}/\kappa_{\rm cub}$ over a wide range by modifying their diameter, thus making them attractive materials for complex phononic and thermoelectric applications/systems.

Tables

Table 1. Table 1 : Calculated κ 𝜅 \kappa and κ h ​ e ​ x / κ c ​ u ​ b subscript 𝜅 ℎ 𝑒 𝑥 subscript 𝜅 𝑐 𝑢 𝑏 \kappa_{hex}/\kappa_{cub} ratios at 300 K. The experimental κ 𝜅 \kappa at room temperature is also presented at normal conditions. The values for all κ 𝜅 \kappa ’s are given in W/m ⋅ ⋅ \cdot K.

	${\kappa }_{cub}^{calc}$ Ref. 19	$\kappa _{hex}^{calc}{}_{†}{}^{}$ Ref. 19	${\kappa }_{cub}^{calc}$ Ref. 20	$\kappa _{hex}^{calc}{}_{}{}^{†}$ Ref. 20	${\kappa }_{cub}^{calc}$ This work	$\kappa _{hex}^{calc}{}_{}{}^{†}$ This work	${}_{}{}^{{\kappa }_{hex}}/_{{\kappa }_{cub}}^{}$	$\kappa _{}^{exp}{}_{}{}^{‡}$
BN	726	592 (602/573)§	940††	–	1071	887.7 (906.6/849.9)§	0.845	760a (cub) 1200b (cub)
AlAs	86.8	72.9 (73.9/71.0)	105	–	100.0	65.84 (65.54/66.45)	0.655	98a (cub) 91c (cub)
GaN	181	171 (171/172)	215	241 (242/239)	290.7	304.3 (293.8/325.3)	1.047	253d (hex) 269e (hex) 294f (hex) 280g (hex) 300g (hex) 330g (hex) 380g (hex)
GaP	104	92.8 (96.5/85.4)	131	–	157.1	144.7 (148.8/136.4)	0.915	77h (cub) 100a (cub) 110i (cub)
GaAs	32.1	27.2 (27.8/25.9)	54	–	47.23	39.52 (38.97/40.61)	0.837	45a (cub) 45.5j (cub) 46k (cub)
InP	85.2	68.9 69.3/68.2	89	–	106.2	87.82 (85.46/92.55)	0.827	93a (cub) 67k (cub) 68h (cub)
InAs	25.2	18.3 (18.5/18.0)	36	–	36.63	33.29 (33.13/33.61)	0.909	30a (cub) 27.3j (cub) 26.5h (cub)
ZnSe	15.6	14.0 (13.8/14.5)	–	–	26.06	22.35 (21.65/23.76)	0.858	19a,l (cub) 33m (cub)

Table 2. Table 2 : Seven standard conditions of high κ 𝜅 \kappa for all materials in both phases. Boldface indicates the characteristic more favorable to a higher κ 𝜅 \kappa . M a ​ v ​ g subscript 𝑀 𝑎 𝑣 𝑔 M_{avg} stands for the average atomic mass of the unit cell, f D subscript 𝑓 𝐷 f_{D} for the Debye frequency and γ 300 ​ K subscript 𝛾 300 𝐾 \gamma_{300K} for the Grüneisen parameter at room temperature.

material	(I) crystal structure	(II) $M_{\mathrm{avg}}$ $\left(amu\right)$	(III) $f_D$ $\left(\mathrm{THz}\right)$	(IV) ${\gamma }_{300K}$	(V) “a-o” gap $\left(\mathrm{THz}\right)$	(VI) “acoustic” bunching $\left(\mathrm{THz}\right)$	(VII) optical bandwidth $\left(\mathrm{THz}\right)$
BN	w zb	12.41 12.41	39.15 39.15	0.67 0.73	0.00 0.00	30.11 15.55	11.41 9.37
AlAs	w zb	50.95 50.95	5.44 5.41	0.29 0.40	3.54 3.55	6.32 3.56	1.93 1.93
GaN	w zb	41.86 41.86	9.75 9.71	0.64 0.70	6.18 6.50	10.27 6.43	5.87 5.36
GaP	w zb	50.35 50.35	6.66 5.57	0.51 0.57	2.67 2.49	7.01 4.53	2.10 1.99
GaAs	w zb	72.32 72.32	4.53 4.67	0.55 0.60	0.00 0.00	6.21 4.35	2.89 2.53
InP	w zb	72.90 72.90	3.51 3.48	0.38 0.48	3.88 3.69	4.92 3.41	1.42 1.35
InAs	w zb	94.87 94.87	2.92 3.07	0.36 0.38	0.82 0.51	4.56 3.62	1.55 1.39
ZnSe	w zb	72.17 72.17	4.70 4.47	0.67 0.71	0.00 0.00	5.29 4.00	2.03 2.06

Table 3. Table 3 : Mean of effective anharmonicity and accessible phase space for zinc-blende and wurtzite phases at 77 and 300K, together with κ 𝜅 \kappa ratios. Boldface indicates the characteristic more favorable to a higher κ 𝜅 \kappa . Mean of effective anharmonicities and accessible phase spaces are given in eV 2 /(amu ⋅ 3 Å 6 {}^{3}\cdot\mathrm{\AA}^{6} ) and ps, respectively.

material	$\overline{{\left|V_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}\right|}_{occ,77K}^2}$	${\delta }_{occ,77K}$	${\kappa }_{hex}^{77K}/{\kappa }_{cub}^{77K}$	$\overline{{\left|V_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}\right|}_{occ,300K}^2}$	${\delta }_{occ,300K}$	${\kappa }_{hex}^{300K}/{\kappa }_{cub}^{300K}$
BN ZB WZ WZ/ZB	$3.387\times {10}^{-10}$ $7.904\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟏}}$ 0.233	$2.324\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$ $1.231\times {10}^{-5}$ 5.297	1.284	$1.350\times {10}^{-4}$ $3.422\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}$ 0.253	$8.757\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$ $7.284\times {10}^{-2}$ 8.318	0.829
AlAs ZB WZ WZ/ZB	$1.538\times {10}^{-7}$ $4.457\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$ 0.290	$1.100\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏}}$ $9.662\times {10}^{-1}$ 8.784	0.617	$4.378\times {10}^{-5}$ $1.051\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}$ 0.240	$1.186\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏}}$ $1.018\times {10}^2$ 8.583	0.658
GaN ZB WZ WZ/ZB	$7.227\times {10}^{-8}$ $1.522\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$ 0.210	$5.877\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$ $4.088\times {10}^{-2}$ 6.956	1.359	$1.157\times {10}^{-4}$ $2.424\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}$ 0.210	$1.217\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟎}}$ $9.302\times {10}^0$ 7.643	1.047
GaP ZB WZ WZ/ZB	$1.631\times {10}^{-7}$ $3.797\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$ 0.233	$6.190\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$ $5.566\times {10}^{-1}$ 8.992	0.882	$7.073\times {10}^{-5}$ $1.471\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}$ 0.208	$8.147\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟎}}$ $7.103\times {10}^1$ 8.719	0.921
GaAs ZB WZ WZ/ZB	$1.914\times {10}^{-7}$ $4.817\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$ 0.252	$2.587\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏}}$ $2.230\times {10}^0$ 8.620	0.916	$3.777\times {10}^{-5}$ $9.459\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$ 0.250	$3.019\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏}}$ $2.559\times {10}^2$ 8.476	0.837
InP ZB WZ WZ/ZB	$1.029\times {10}^{-7}$ $2.501\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$ 0.243	$3.581\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏}}$ $3.030\times {10}^0$ 8.461	0.741	$2.115\times {10}^{-5}$ $4.692\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$ 0.222	$3.114\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏}}$ $2.610\times {10}^2$ 8.382	0.827
InAs ZB WZ WZ/ZB	$2.226\times {10}^{-7}$ $5.220\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$ 0.235	$8.224\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏}}$ $6.651\times {10}^0$ 8.087	0.832	$2.966\times {10}^{-5}$ $6.776\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟔}}$ 0.228	$6.909\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏}}$ $5.548\times {10}^2$ 8.030	0.909
ZnSe ZB WZ WZ/ZB	$3.025\times {10}^{-7}$ $7.207\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$ 0.238	$4.817\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏}}$ $3.548\times {10}^0$ 7.366	0.971	$4.648\times {10}^{-5}$ $1.105\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟓}}$ 0.238	$4.902\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏}}$ $3.689\times {10}^2$ 7.25	0.857

Table 4. Table 4 : Anharmonicity ratio | Φ 4 / Φ 3 | 2 / | Φ 2 | superscript subscript Φ 4 subscript Φ 3 2 subscript Φ 2 \left|\Phi_{4}/\Phi_{3}\right|^{2}/\left|\Phi_{2}\right| to get an insight of the relative importance of four-phonon processes in GaN-ZB and GaN-WZ. IFCs have been calculated along the stacking direction ([111] or [0001]), and they are given in eV/Å n .

phase	${\Phi }_2$	${\Phi }_3$	${\Phi }_4$	${\left|{\Phi }_4/{\Phi }_3\right|}^2/\left|{\Phi }_2\right|$
GaN-ZB	20.1	70.0	244	0.604
GaN-WZ	19.7	69.3	184	0.358

Equations

$${\kappa }^{\alpha \beta }=\frac{1}{N{k}_B\Omega T^2}\sum _{\lambda }{n}_0(n_0+1)(\hslash {\omega }_{\lambda }{)}^2{\upsilon }_{\lambda }^{\alpha }{F}_{\lambda }^{\beta },$$

$${\Gamma }_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{\pm }=\frac{\hslash \pi |V_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{\pm }{|}^2}{4{\omega }_{\lambda }{\omega }_{{\lambda }^{\prime }}{\omega }_{{\lambda }^{\prime \prime }}}\left\{\begin{array}{c}{n}_0^{\prime }-n_0^{\prime \prime }\\ n_0^{\prime }+n_0^{\prime \prime }+1\end{array}\right\}\delta ({\omega }_{\lambda }\pm {\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{\prime \prime }}),$$

$${\delta }_{occ,T}=\frac{2}{3}({\delta }_{occ,T}^{+}+\frac{1}{2}{\delta }_{occ,T}^{-})$$

$$\begin{gathered} {\delta }_{occ,T}^{+}=2\pi {\left[\frac{V_{BZ}}{(2\pi {)}^3}\right]}^3\sum _{n,n^{\prime },n^{\prime \prime }}\underset{BZ}{\iiint }\delta ({\omega }_{\lambda }+{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{\prime \prime }})\times \\ {n}_0{n}_0^{\prime }(1+n_0^{\prime \prime }){\delta }_{q+q^{\prime },q^{\prime \prime }+G}{d}^3{q}^{\prime \prime }{d}^3{q}^{\prime }{d}^{3}q, \end{gathered}$$

$$\begin{gathered} {\delta }_{occ,T}^{-}=2\pi {\left[\frac{V_{BZ}}{(2\pi {)}^3}\right]}^3\sum _{n,n^{\prime },n^{\prime \prime }}\underset{BZ}{\iiint }\delta ({\omega }_{\lambda }-{\omega }_{{\lambda }^{\prime }}-{\omega }_{{\lambda }^{\prime \prime }})\times \\ {n}_0(1+n_0^{\prime })(1+n_0^{\prime \prime }){\delta }_{q-q^{\prime },q^{\prime \prime }+G}{d}^3{q}^{\prime \prime }{d}^3{q}^{\prime }{d}^{3}q, \end{gathered}$$

$$\overline{|V_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}{|}_{occ,T}^2}=\frac{2}{3}(\overline{|V_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{+}{|}_{occ,T}^2}+\frac{1}{2}\overline{|V_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{-}{|}_{occ,T}^2})$$

$$|V_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{+}{|}_{occ,T}^2$$

$$|V_{\lambda {\lambda }^{\prime }{\lambda }^{\prime \prime }}^{-}{|}_{occ,T}^2$$

Full text

I Introduction

Heat management stands out as one of the key problems for several technologies. The reduction in size Moore (1965) and the required increment in operating speeds of devices in electronics, or the necessity to enhance energy conversion ratio for thermoelectrics have increased the importance of phonon engineering or phononics Cahill et al. (2014). For instance, for thermoelectric applications a material with low thermal conductivity and high electrical conductivity, the so called “phonon-glass and electron-crystal”, is desired in order to have a good efficiency Slack (1995). Historically, one of the preferred approaches to thermoelectric materials was alloying, as it introduces point-mass defects that strongly scatter phonons, significantly reducing lattice thermal conductivity Kanatzidis (2010). A more recent approach to solve this problem is the nanostructuration as an effective way of reducing the lattice thermal conductivity Maldovan (2013)—hereafter just referred to as thermal conductivity, $\kappa$—or a combination of both Kanatzidis (2010).

At the same time, recent advances in semiconductor synthesis at the nanoscale have granted access to different phases that, in bulk, are only observed under extreme conditions. Namely, the wurtzite phase has theoretically been proved to be the most stable for several bulk-cubic Yeh et al. (1992) semiconductors when grown as nanowires (NWs) Akiyama et al. (2006) provided their diameter is smaller than a given critical value, but it can also be obtained at larger sizes, though thermodynamically it is only metastable. Indeed, the wurtzite phase in NWs has been reported for a great range of bulk-cubic semiconductors: InP Mattila et al. (2006); Paiman et al. (2009), Si Lopez et al. (2011); Hauge et al. (2015); Tang et al. (2017), InAs Hiruma et al. (1995); Caroff et al. (2009), GaAs Tchernycheva et al. (2006); Zardo et al. (2009) or GaP Assali et al. (2013); Berg et al. (2014) to name but a few.

It is in that context where these new semiconductor phases can become a cornerstone for thermoelectricity and phononics in general, as they normally show a lower thermal conductivity than their more stable counterparts while maintaining similar electronic properties. For example, Togo et al. calculated thermal conductivity of wurtzite and zinc-blende phases for 33 different materials Togo et al. (2015), Lindsay et al. Lindsay et al. (2013a) studied the zinc-blende/wurtzite thermal conductivity of GaN, and Li et al. addressed the thermal conductivity of bulk and nanowire InAs, AlN, and BeO polymorphs Li and Mingo (2013). However, none of these works explained in detail the origin of the reduction in the thermal conductivity. In a more recent example, we have reported and discussed a reduction of $40\%$ in the thermal conductivity between the more stable 3C cubic-Si and 2H hexagonal-Si Raya-Moreno et al. (2017).

With regard to the qualitative behavior of the thermal conductivity of materials, Mukhopadhyay et al. Mukhopadhyay et al. (2016), building on earlier work by Slack Slack (1973) and Lindsay et al. Lindsay et al. (2013b), have provided seven criteria to interpret the relative magnitude of the thermal conductivity between different materials.

In this work we study the thermal conductivity ($\kappa$) of several bulk and nanowire semiconductors for their cubic (zinc-blende, ZB) and the hexagonal (wurtzite, WZ) phases from first-principles. After seeing that the seven criteria provided in Ref. Mukhopadhyay et al., 2016 are unsuitable for rationalizing the $\kappa_{\rm hex}/\kappa_{\rm cub}$ value for the different materials, we propose an approach that successfully addresses this issue, providing insight into the factors determining the thermal conductivity of the materials we have studied.

The paper is structured as follows: after discussing the applied methodology for the thermal conductivity calculation in Sec. II, we present the results in Sec. III. Results are provided for bulk systems (Sec. III.2) and nanowires (Sec. III.3). Summary and conclusions are given in Sec. IV.

II Methodology

Harmonic and anharmonic interatomic force constants (IFCs), needed to calculate the thermal conductivity, were obtained using the supercell method, as implemented in Phonopy Togo and Tanaka (2015) for the harmonic IFCs and thirdorder.py Li et al. (2014) for anharmonic IFCs. They were calculated in cubic polytypes in a 5x5x5 supercell for harmonic IFCs and 4x4x4 for the anharmonic IFCs. For hexagonal polytypes, the used supercell was 4x4x3 for both types of IFC. In order to minimize the computational burden, anharmonic IFCs were computed from interactions up to fourth nearest neighbors, while it has been previously reported that including up to third nearest neighbors was sufficient to give a satisfactorily converged value of $\kappa$ Li et al. (2014). The unit cells used to span the supercells were optimized until strict limits for stress ($3\cdot 10^{-3}$GPa) and forces ($5\cdot 10^{-4}$ eV/Å) were attained. These optimizations were conducted within density-functional theory (DFT), using the VASP Kresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b) code with projector augmented-wave (PAW) potentials Blöchl (1994); Kresse and Joubert (1999). Local density approximation (LDA) for the exchange-correlation as parametrized by Perdew and Zunger Perdew and Zunger (1981) to Ceperley-Alder Ceperley and Alder (1980) data was used. For each system, the $k$-point mesh size had been previously optimized, taking into account that the supercell used to calculate 3rd-order IFCs should be commensurate with the mesh in order not to introduce spurious forces. Therefore, the selected primitive cell $k$-mesh for cubic materials was a 16$\times$16$\times$16 shifted mesh, and a 16$\times$16$\times$12 $\Gamma$-centered mesh for hexagonal polytypes, which are converged meshes for all systems. After optimization, a density functional perturbation theory (DFPT) run using VASP, with a doubled $k$-mesh was performed to obtain the Born charges ($Z^{*}$) and dielectric constant at high frequency ($\varepsilon^{\infty}$), needed to calculate the non-analytic term correction for the dynamical matrix near $\Gamma$.

After obtaining the IFCs, we solved the full linearized Boltzmann Transport Equation (LBTE) iteratively, as implemented in the almaBTE Carrete et al. (2017) code for bulk and the ShengBTE code Li et al. (2014) for nanowires. Consequently, the thermal conductivity tensor ($\kappa^{\alpha\beta}$) is obtained as:

[TABLE]

where $\alpha$ and $\beta$ label the Cartesian axes $x$, $y$ and $z$; ${N,k_{B},\Omega,T}$ are the number of q-points, the Boltzmann constant, the cell volume and the temperature, respectively. The summation is done over all phonon modes $\lambda$, which are characterized by the band label $n$ and their $q$-vector. $n_{0}$ is the equilibrium Bose-Einstein distribution function, $\hbar$ is the reduced Planck constant, $\omega_{\lambda}$ is the phonon mode frequency and $\upsilon_{\lambda}^{\alpha}$ is the group velocity of phonon mode $\lambda$ along the $\alpha$ direction. $F_{\lambda}^{\beta}$ is the generalized mean free path of the phonon mode along the $\beta$ direction, and it is calculated as $\tau_{\lambda}(\upsilon_{\lambda}^{\beta}+\Delta_{\lambda}^{\beta})$. $\tau_{\lambda}$ is the lifetime of the phonon in the relaxation time approximation (RTA) and $\Delta_{\lambda}^{\beta}$ is a measure of how much the population of a specific phonon mode and its associated heat current deviates from the RTA prediction Li et al. (2014). This correction is obtained iteratively starting from the RTA $(\Delta_{\lambda}^{\beta}=0)$. In addition to anharmonic three-phonon scattering processes, mass-disorder scattering coming from isotope presence is treated using the Tamura model Tamura (1983). Boundary scattering, whose effect is negligible in the temperature range considered, is not explicitly accounted for in the bulk systems (see Sec. III.3 for a short description on the implementation of boundary scattering in NWs).

Regarding convergence with $q$-mesh, we found (see Figs. LABEL:conv-WZ and LABEL:conv-ZB in the Supplemental Material SI ) that for all materials but GaN and BN a $q$-mesh of 30$\times$30$\times$30 (ZB) and 30$\times$30$\times$19 (WZ) is enough to obtain converged values—less than 5% change with respect to a higher accuracy 34$\times$34$\times$34/21 (ZB/WZ) mesh—of $\kappa$ and $\kappa_{pure}$ (i.e. without isotopes) for the 50–1000 K range. For GaN the converged $q$-mesh is found to be 34$\times$34$\times$34 (ZB) and 34x34x21 (WZ) for $\kappa$ and $\kappa_{pure}$, with respect to a 38$\times$38$\times$38/24 (ZB/WZ) mesh. Finally, the harder phonon modes of BN translate into a more demanding convergence behavior. On one hand, $q$-mesh convergence was only achieved for a 38$\times$38$\times$38 and 38$\times$38$\times$24 mesh, but only for the non-isotopically pure material, while the iterative calculation did not reach convergence for the isotopically pure material at low temperatures. Note that all these $q$-mesh values are significantly denser than the 24$\times$24$\times$24 mesh needed to achieve convergence in Si Li et al. (2014).

III Results

III.1 Interatomic Force Constants test

Owing to the computational cost required to obtain the anharmonic IFCs (168 and 208 DFT runs for cubic and hexagonal polytypes, respectively), prior to their calculation an accuracy test to harmonic IFCs was done using Phonopy Togo and Tanaka (2015); Wang et al. (2010); Togo et al. (2015), by checking the optical phonon frequencies at the $\Gamma$ point. These results, comparing calculated and experimental values, are shown in the phonon dispersion relations (see Figs. LABEL:disp1-LABEL:disp3 SI ), exhibiting, despite small differences, a good agreement between our calculations and the experimental values for the TO modes. The non-analytical correction (NAC), needed in polar materials to get the LO-TO splitting, is underestimated because of it being inversely proportional to $\epsilon^{\infty}$ Pick et al. (1970), which is overestimated due to LDA inability to take into account the polarization dependence on the exchange-correlation functional under a field Ortiz et al. (1998). This leads to an underestimation of the LO frequency with respect to the experiments. It should be noted, nonetheless, that the effect of LDA and GGA shortcomings have been proven to have a smaller effect in the anharmonic properties than in the harmonic ones Narasimhan and de Gironcoli (2002); Arrigoni and Madsen (2019).

III.2 Bulk

In this section we study the thermal conductivities of cubic (ZB) and hexagonal (WZ) phases for different compound semiconductors, using the methodology discussed in Sec. II. These two phases have the same coordination (4-fold, tetrahedral) and very similar first neighbor distances, the main difference being the stacking of atomic layers.

We further substantiate our results by comparing our calculated thermal conductivity to available experimental values (see Table 1). It is reasonable to expect that, for some materials, our results overestimate the experimental values, because samples used in such experiments might contain defects (impurities, vacancies, dislocations, etc), which can strongly suppress thermal conductivity, and they have not been considered in our simulations. Moreover, there might be a dependence of the measured value on the experimental technique, and it is quite challenging to obtain experimental thermal conductivity values with less than 5% error Zhao et al. (2016). We note that our values are in excellent agreement with experimental results for GaAs and AlAs, while keeping a good agreement for InP, InAs, ZnSe and GaP. In the nitrides, BN and GaN, we obtain values within the dispersion of the experimentally reported magnitudes. Regarding comparison with other first-principles calculations, our results are in reasonable agreement with those of Lindsay et al. Lindsay et al. (2013a), while there is a stronger disagreement with the values of Togo et al. Togo et al. (2015), which are obtained by a different approach to the LBTE. These two approaches are known to provide different values for the thermal conductivity in transition metal dichalcogenides as well Torres et al. (2019).

We found (see Table 1 and Figs. 1-8) that most of materials under study (GaAs, GaP, InP, InAs, ZnSe and AlAs) follow the silicon behavior of reducing their $\kappa$ with symmetry Raya-Moreno et al. (2017). However, we observe that is possible for some materials (GaN) to have the opposite behavior, namely $\kappa$ increases when symmetry is reduced, i.e. going from ZB to WZ. Moreover, as previously observed in other materials Li and Mingo (2013), BN can show these two opposite behaviors at different temperatures.

III.2.1 Standard analysis of $\mathrm{\kappa_{hex}/\kappa_{cub}}$

Four conditions that a crystal must fulfill to have a high $\kappa$ have been long established  Slack (1973): (I) to be structurally simple, (II) to be composed of light elements, (III) to have strong covalent bonds—represented normally by a high Debye temperature/frequency ($f_{D}$)—and (IV) to be harmonic, which is normally associated to a low Grüneisen parameter ($\gamma$).

Recently, three additional conditions for a high $\kappa$ have been introduced by Lindsay et al. Lindsay et al. (2013b) and Mukhopadhyay et al. Mukhopadhyay et al. (2016). They demonstrated that, in addition to the four previous conditions, a higher $\kappa$ is obtained if the material also has (V) high “acoustic-optical” (a-o) gap, (VI) a high “acoustic bunching” and (VII) low optical bandwidth.

We have collected those seven conditions for all materials under study, with values from our first-principles calculations, in Table 2. Interestingly, by observing those conditions together with the $\kappa$ values, we see that they do not clarify which of the phases—cubic or hexagonal—is the most conductive: for instance, if one looks at GaN four out of seven criteria suggest that zinc-blende phase should be the most conductive one, while one of them is neutral and only two of them predict a larger $\kappa$ for the wurtzite. Yet, the calculations predict the latter to be more conductive. Additionally, those conditions cannot either explain the reason for some materials, like BN, to change their $\kappa_{hex}/\kappa_{cub}$ ratio behavior with temperature. Such a shortcoming is due to being based on arguments of a qualitative character that do not quantify the relative importance of each condition, thus the necessity of a more quantitative viewpoint to predict/understand which phase is the most conductive at a given temperature.

III.2.2 Effective anharmonicity and accessible phase space

To gain insight on the $\kappa_{hex}/\kappa_{cub}$ ratio at different temperatures, we focus our analysis on two quantities that together contain all conditions: the three-phonon scattering matrix elements or anharmonicity (high-$\kappa$ conditions I,II, III and IV) and phase space, i.e. all the energy conserving three-phonon combinations (high-$\kappa$ conditions I,II,V,VI and VII). In fact, these two quantities directly contribute to $\kappa$ via the three-phonon scattering rate ($\Gamma_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}^{\pm}$) Li et al. (2014):

[TABLE]

where $V^{\pm}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}$ stands for the three-phonon scattering matrix element (+ for absorption processes and - for emission) and $n_{0}^{\prime}$ is a shorthand for $n_{0}(\omega_{\lambda^{\prime}})$ and similarly for $n_{0}^{\prime\prime}$.

As it can be seen from Eq. (2), an increment in anharmonicity ($V^{\pm}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}$) produces an increment in the scattering rate, therefore a reduction in $\kappa$. In the same way, an increment in phase space, represented in Eq. (2) by the energy conservation delta $\delta(\omega_{\lambda}\pm\omega_{\lambda^{\prime}}-\omega_{\lambda^{\prime\prime}})$, also reduces $\kappa$.

Comparing the phonon dispersion of both phases for all materials (see Figs. LABEL:disp1-LABEL:disp3 SI ) it becomes obvious that symmetry reduction, which causes the appearance of new low/medium optical eigenmodes with non-vanishing scattering matrix elements, also increases the phase space for transitions at a given temperature. In order to obtain an actual measurement of accessible phase space while taking into account the effect of temperature, we calculate it ($\delta_{occ,T}$) defined as follows:

[TABLE]

[TABLE]

[TABLE]

where $\delta_{q\pm q^{\prime},q^{\prime\prime}+G}$ is the momentum conservation condition; $\delta^{+}_{occ,T}$ and $\delta^{-}_{occ,T}$ are the contribution to the accessible phase space of absorption and emission processes, and the $2/3$ and $1/2$ weighting factors ensure the normalization and non double counting of processes Lindsay and Broido (2008). Therefore, $\delta_{occ,T}$ gives an idea for a given temperature of the accessible part of the available phase space. Certainly, the attainable phase space for the different materials (see Table  3 and Figs. LABEL:DELTAS1-LABEL:DELTAS3 SI ) confirms that, as inferred from dispersion relations, the hexagonal phase has a greater accessible phase space which contributes to reducing the $\kappa_{hex}/\kappa_{cub}$ ratio. This is a general feature of all the studied materials.

On the other hand, from the scattering matrix elements ($|V^{\pm}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}|^{2}$) of energy-conserving three-phonon processes one obtains direct information of the material anharmonicity. Notwithstanding that this bare anharmonicity is an interesting quantity by itself, it is not useful for us because it gives the same importance to processes in which the involved modes are occupied to those in which they are not. Therefore, analogously to what we did with the accessible phase space, we define the temperature-dependent mean effective anharmonicity as:

[TABLE]

where $\overline{|V^{+}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}|^{2}_{occ,T}}$ and $\overline{|V^{-}_{\lambda\lambda^{\prime}\lambda^{\prime\prime}}|^{2}_{occ,T}}$ are the arithmetic means of the population weighted three-phonon matrix elements squared modulus for absorption and emission processes:

[TABLE]

From the mean effective anharmonicity, it can be seen that, for all materials under study, the cubic phase is more anharmonic than the hexagonal phase both at 77 K (a representative value of the low temperature regime) and at 300 K (see Table 3 and Figs. LABEL:vpp_fulls1-LABEL:vpp_fulls3 SI ).

Therefore, we have two antagonistic processes occurring when reducing the symmetry from cubic to hexagonal: an increment in the available phase space for phonon-phonon scattering events and a lowering of the anharmonicity, which makes the strength of those events weaker when compared to the cubic ones, with no indication of their relative importance as regards to $\kappa$.

III.2.3 Discussion

As we already mentioned, despite being widely used to predict the relationship between thermal conductivity for different materials Mukhopadhyay et al. (2016), the criteria listed in Sec. III.2.1 lack the capacity to discern the relative importance of opposed processes.

To overcome such a limitation and owing to the fact that scattering rates are a product of the anharmonicity with phase space [see Eq. (2)], we represent in Fig. 9 $\kappa_{hex}/\kappa_{cub}$ versus the hexagonal-cubic ratio of the mean effective anharmonicity, Eq. (6), and the accessible phase space product, Eq. (3) at 77K and 300K. We call this ratio, which is a central magnitude in our discussion, REAAPS, standing for Ratio of Effective Anharmonicity and Accessible Phase Space product. Isotopic scattering, despite being the clear mechanism determining the ratio at low temperature for some materials like GaP, is not accounted in such analysis. To avoid the influence of the isotope effect, the same procedure, but with the ratios of isotopically pure materials, can be repeated (see Fig. 10).

By this analysis of the relative importance of the accessible phase space and effective anharmonicity at a given temperature, it can be observed from Fig. 10 that the $\kappa_{hex}^{pure}/\kappa_{cub}^{pure}$ ratio shows a strong correlation with the REAAPS, and that $\kappa_{hex}^{pure}>\kappa_{cub}^{pure}$ when REAAPS is lower than $\sim 1.66$, indicating that in these cases the effective anharmonicity increment in the zinc-blende phase compensates the wurtzite accessible phase space increment, thus giving a higher conductivity in wurtzite phase. For higher values of REAAPS, the dominant factor is the increment in the accessible phase space of the wurtzite phase, hence the higher conductivity of the zinc-blende phase. Interestingly, the plots with (Fig. 9) and without (Fig. 10) isotopic scattering are quite similar at both temperatures, indicating that three-phonon processes are the ones dominating the ratio behavior at the studied temperatures. We note that both the threshold of $\sim 1.66$ and the slopes of the linear regressions are independent, to a large extent, of temperature (see Fig. LABEL:reaaps700 SI for a plot at a high temperature, 700 K).

Our analysis helps to explain some of the ratio behaviors such as in GaN, which has $\kappa_{hex}>\kappa_{cub}$ over all the temperature range considered. As observed in Table 3, the anharmonicity of the hexagonal phase is particularly small when compared to that of the cubic phase. Additionally, the ratio of the accessible phase space is somewhat lower than in other materials. Both factors favor a higher $\kappa$ for the WZ phase, as it is indeed the case. This is in agreement with previous works by Lindsay et al. Lindsay et al. (2013a) but in opposition to the results by Togo et al. Togo et al. (2015), which, as mentioned before, are obtained by a different approach to the LBTE that yields different values for the thermal conductivity in other materials as well Torres et al. (2019).

The case of BN is also interesting. Both Table 3 and Fig. 9 show that there is nothing particular to BN at 300K, having $\kappa_{hex}^{300K}<\kappa_{cub}^{300K}$. However, when the temperature is lowered to 77K, the accessible phase space decreases much more than in the other materials (a consequence of the hard phonon modes) and the ratio of accessible phase space takes a low value of 5.297, significantly different from the rest of the materials. From this it can be concluded that the change from $\kappa_{hex}^{300K}<\kappa_{cub}^{300K}$ to $\kappa_{hex}^{77K}>\kappa_{cub}^{77K}$ is due to the abnormally large increase of $\kappa_{hex}$ as temperature is decreased because of the quicker decrease of the accessible phase space in the WZ phase.

Moreover, one can also see from Table 3 that for the vast majority of materials the behavior of the ratio with temperature is the opposite to the behavior of REAAPS with temperature. Therefore, when REAAPS decreases (increases) with the temperature, the ratio increases (decreases) for almost all materials, explaining the temperature dependence of $\kappa_{hex}^{pure}/\kappa_{cub}^{pure}$. However, this does not occur in the case of GaAs, where contrarily to the rest of the materials, the evolution of the REAAPS seems to indicate that the ratio should increase with temperature. This disagreement may show a limitation of our analysis, where we have not taken into account the correlations between anharmonicity and accessible phase space, which can play an important role in borderline cases. When isotopic scattering is added into consideration, it can significantly alter the monotonic dependence of the ratio vs. the temperature at very low temperatures, see the case of GaP for an example.

III.2.4 Four-phonon scattering for GaN

Feng et al. Feng et al. (2017) have conducted a rigorous study of four-phonon scattering on three representative materials, showing that a significant reduction of $\kappa$ is caused at high temperatures. This effect was shown to be particularly strong in BAs, a material with a large a-o gap, where the inclusion of four-phonon processes opened up scattering channels that were forbidden in a three-phonon event. GaN also presents a large a-o gap, and thus it is interesting to estimate how four-phonon scattering can affect the predicted values for $\kappa$.

Although a full study is out of scope of the present work, Ref. 61 provides some guidelines on how to estimate the effect of four-phonon scattering on the thermal conductivity. We have followed the procedure there reported for the estimation, computing the anharmonicity ratio $\left|\Phi_{4}/\Phi_{3}\right|^{2}/\left|\Phi_{2}\right|$, where $\Phi_{n}$ is the $n$-th order on-site force constant for Ga along the stacking direction ([111] or [0001]), e.g. $\Phi_{3}=\Phi^{\rm WZ}_{0,{\rm Ga},z;0,{\rm Ga},z;0,{\rm Ga},z}$, and comparing the cubic and hexagonal phases. These directions are chosen because (a) they correspond to the directions of the cation-anion bond, (b) they are the directions maximum structural difference between ZB and WZ, and (c) the IFCs do not take zero values. The results, shown in Table 4, indicate that the inclusion of four-order processes maintains that the anharmonicity of GaN-ZB is stronger than that of GaN-WZ. Also, those anharmonicity ratios, although the IFCs were computed along different directions than in Ref. 61 and we do not know up to what point they can be directly compared, have a higher numerical value than those provided for diamond, BAs and Si Feng et al. (2017). Furthermore, given that GaN-ZB has a larger a-o gap than its WZ counterpart, we expect that at high temperatures ZB will be more affected by the inclusion of four-phonon processes than the hexagonal phase. From all these considerations, we can predict that $\kappa_{hex}/\kappa_{cub}$ will increase at high temperature once four-phonon processes are included in the analysis.

III.3 Nanowires

As previously mentioned, phases that are not thermodynamically stable at room temperature and atmospheric pressure in bulk form can naturally occur when the materials are grown as nanowires, and thus it is in nanowires that the zinc-blende and wurtzite phase can both be easily accessed. Therefore, in this section we discuss the conductivity along the stacking direction, which is also the common growth direction ([111] for cubic and [0001] for hexagonal) for nanowires of different materials and diameters.

Despite being one dimensional structures, nanowire phonon dispersions can be approximated to the bulk ones for nanowires with diameters $\gtrsim 60-70$ nm Li et al. (2014). However, the reduction in the symmetry caused by the nanowire boundaries make it unrealistic to use bulk-like scattering rates, as they become position dependent in the direction perpendicular to the nanowire. To solve this problem, an iterative solution under the diffusive regime as proposed by Li et al. Li et al. (2012) is used as implemented in the ShengBTE code Li et al. (2014).

$\kappa$ along [111] (cubic) and [0001] (hexagonal) is plotted as function of the nanowire diameter for several temperatures in Fig. 11 together with $\kappa_{hex}/\kappa_{cub}$. In order to understand the ratio behavior in NWs, it becomes essential to obtain an insight of the size effects (boundary scattering). To do so, the cumulative thermal conductivity with the phonon mean free path (MFP) was plotted together with the nanowire $\kappa$ for both phases at different temperatures (see Fig. 12). Thus, from Fig. 12 we can confirm that the ratio in NWs is mostly controlled by size effects. Although the cumulative thermal conductivity with MFP do not hold accurate predictive power Li et al. (2012), they overall reproduce the behavior of the ratios as a function of diameter (Fig. LABEL:MFP-ratios SI ), especially at high temperatures 111The lack of agreement at 77K might be due to the larger number of phonons (lower temperatures mean longer MFPs) affected by the different cutoff behaviors for a given MFP/diameter.. For instance, InAs having a ratio larger than 1 at small diameters is due to the behavior of the ratio of the cumulative functions. Notwithstanding its utility for determining size effects (i.e: boundary scattering) trends, it is also clear that MFP plot cannot be used to perfectly predict the detailed ratio behavior as the cumulative and NW function differ, due to the cumulative function not accounting for phonon propagation axis or the expression of boundary scattering in NWs having an $\exp(-d/\lambda_{\rm MFP})$ behavior—where $d$ is the distance traveled by the phonon hitting the NW surface and $\lambda_{\rm MFP}$ is its MFP—as opposed to an abrupt cutoff Li et al. (2012).

Moreover, from the nanowire $\kappa$s we can also observe that BN and GaN show a higher $\kappa$ at 300K than at 77K, thus indicating a displacement in the $\kappa(T)$ peak to higher temperatures when compared to the bulk. Such results are coherent with experimental observations for silicon NWs Li et al. (2003), and they are associated to the domination of boundary scattering over Umklapp scattering in the NW geometry when compared to the bulk one. Finally, it is worth mentioning that the available tuning of $\kappa$ ratios by modifying their diameter, especially for AlAs with a range between 0.7 and 1.1, makes these NWs interesting blockpieces for complex thermoelectric and/or phononic systems.

IV Summary and conclusions

In this work we have presented ab initio calculations of the lattice thermal conductivity for the cubic and hexagonal phases of GaAs, GaN, GaP, InAs, InP, AlAs, BN and ZnSe using density functional theory and iteratively solving the phonon Boltzmann Transport Equation. We find that which phase is the most conductive one depends on the relative strength between effective anharmonicity and accessible phase space. Such factors are shown to be antagonistic for all materials due to the higher effective anharmonicity of the cubic phase when compared to the hexagonal one, which, on the other hand, has a higher accessible phase space. Furthermore, we carry out an analysis of which factor is dominant when three phonon processes are the only ones present in each material at 77 and 300K, showing that, when anharmonicity (phase space) dominates, it leads to a higher (lower) conductivity in the hexagonal phase compared to the cubic one. Moreover, we have observed that when the hexagonal-cubic ratio of temperature-weighted anharmonicity and accessible phase space product is less than $\sim 1.66$, the dominating factor determining $\kappa$ is anharmonicity ($\kappa_{hex}>\kappa_{cub}$). On the contrary, when that product is higher than $\sim 1.66$, the dominating factor determining $\kappa$ is the accessible phase space($\kappa_{hex}<\kappa_{cub}$), thereby making such quantity an excellent tool to predict which is the most conductive phase at a given temperature when other more qualitative analysis fail. We also present results for NWs, showing the effect of boundary scattering on ($\kappa_{hex}/\kappa_{cub}$) and its relation to the phonon mean free paths. Moreover, we find that NWs have the ability to have their $\kappa$ ratio tuned over a wide range with diameter hence making them appealing materials for phononic and thermoelectric applications. The case of AlAs, with a $\kappa_{hex}/\kappa_{cub}$ range between 0.7 and 1.1, is of special interest.

Acknowledgements.

We acknowledge financial support by the Ministerio de Economía, Industria y Competitividad (MINECO) under grants TEC2015-67462-C2-1-R (MINECO/FEDER) and RTI2018-097876-B-C21 (MICINN/FEDER) and FEDER-MAT2017-90024-P, the Severo Ochoa Centres of Excellence Program under Grant SEV-2015-0496, the Generalitat de Catalunya under grant no. 2017 SGR 1506 and the Ministerio de Educación, Cultura y Deporte programme of Formación de Profesorado Universitario under Grant Nos. FPU2016/02565. The authors thankfully acknowledge the computer resources, technical expertise and assistance provided by the Supercomputing and Visualization Center of Madrid (CeSViMa) and the Spanish Supercomputing Network (RES) under project FI-2018-1-0027.