- Email: [email protected]
Phonon transport properties of bulk and monolayer GaN from first-principles calculations
Computational Materials Science 138 (2017) 419–425
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Phonon transport properties of bulk and monolayer GaN from first-principles calculations Yongqiang Jiang, Shuang Cai, Yi Tao, Zhiyong Wei ⇑, Kedong Bi, Yunfei Chen Jiangsu Key Laboratory for Design & Manufacture of Micro/Nano Biomedical Instruments and School of Mechanical Engineering, Southeast University, Nanjing 211189, People’s Republic of China
a r t i c l e
i n f o
Article history: Received 24 February 2017 Received in revised form 27 June 2017 Accepted 11 July 2017
Keywords: Phonon Thermal conductivity GaN First-principles Mean free path
a b s t r a c t The thermal conductivity of isotopically pure wurtzite GaN and the corresponding two-dimensional monolayer crystals are investigated based on first-principles calculations and phonon Boltzmann transport equation. It is found that the c-axis thermal conductivity of defect-free wurtzite GaN is 42–11% higher than that along the in-plane direction when temperature increases from 100 to 800 K. Our calculation also shows that the thermal conductivity of monolayer GaN with hexagonal lattice is higher than that with haeckelite lattice by a factor of four at room temperature. The phonon dispersion relation, phonon group velocity and the phonon scattering rate are extracted to uncover the underlying physical mechanism. The thermal conductivity accumulation function with respect to the phonon mean free path and the relative contribution of each phonon branch for the two kinds of monolayer GaN are further calculated as compared to that of graphene. These investigations would provide important guidance to develop GaN-based electronic and photonic device. Ó 2017 Elsevier B.V. All rights reserved.
1. Introduction Gallium nitride (GaN) is an important wide band-gap semiconductor that can be widely used in electronics and optoelectronics devices, such as microwave power transistors [1,2], lightemitting diodes [3] and laser diodes [4]. The most stable GaN in nature (wurtzite GaN) belongs to hexagonal crystal system, which is the most anisotropic among the seven kinds of crystal family [5]. However, few investigations are conducted to systematically study its anisotropic thermal conductivity. Most previous experiments results [6–8] showed that the thermal conductivity of wurtzite GaN is isotropic around room temperature due to its nearly isotropic elastic property. On the other hand, most theoretical works [9–11] often take the average value of the three directions or the value of one arbitrary direction as the thermal conductivity of GaN at 300 K, and omit its anisotropy. By performing nonequilibrium molecular dynamics simulation, recently Ju et al. [12] shows the thermal conductivity along the c-axis is 25% higher than that along the in-plane direction at room temperature. Although the phonon focusing effects [13] are used to explain the anisotropy, the anisotropic thermal conductivity in wurtzite GaN is still needed to be verified.
⇑ Corresponding author. E-mail address: [email protected] (Z. Wei). http://dx.doi.org/10.1016/j.commatsci.2017.07.012 0927-0256/Ó 2017 Elsevier B.V. All rights reserved.
In the past decade, the exfoliated single layer of graphene from bulk graphite by Novoselov and Geim [14,15] has attracted much attention for the study of atomic layer structure, such as BN, MoS2, black phosphorus, etc. [16–18]. Since there are no weak van der Waals interactions in bulk GaN, it is difficult to exfoliate single layer GaN from its corresponding bulk structure. Thus, little attention is focused on the physical properties of monolayer GaN. However, recent experiment has shown that the atomic layer GaN can be obtained by a migration-enhanced encapsulation growth (MEEG) technique [19,20]. This experimentally stable monolayer GaN is similar to graphene, which is verified by the scanning electron microscope (SEM) and scanning transmission electron microscope (STEM) [19,20]. Sarma et al. [21] shows that the thermal conductivity of this monolayer hexagonal GaN is about 25 W/m-K at room temperature by molecular dynamics simulations. Very recently a new allotrope of monolayer GaN consisting of square and octagonal rings (called Haeckelite 8–4) was proposed and Dulce et al. [22] demonstrated that both Hexagonal (Hexagonal-GaN) and Haeckelite (Haeckelite-GaN) structures are thermally stable. Understanding the phonon transport is crucial to predict the thermal performance of newly emerging materials. According to the lattice difference between the two kinds of monolayer GaN, it is expected that their thermal conductivity may be different. The electronic properties of single layer GaN has been extensively investigated recently, however, the thermal transport properties of them are not systematically investigated so far.
420
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
In this study, we systematically investigate the phonon transport properties of wurtzite GaN crystal and its lattice thermal conductivity by combining first-principles calculations [23] with phonon Boltzmann transport equation (BTE) [24]. The calculated results show that the lattice thermal conductivity of defect-free wurtzite GaN exhibits obvious anisotropic behavior at low temperature. At room temperature, the thermal conductivity of defectfree wurtzite GaN in the [0 0 1] direction is 13% higher than that in the [1 0 0] direction. We also study the thermal behavior of two monolayer GaN lattices by the same method. The phonon properties are also analyzed comprehensively based on phonon dispersion relation, phonon group velocity, phonon lifetimes and specific heat to reveal differences in the thermal transport properties of two kinds of atomic layer. 2. Computational details The first-principles calculations of GaN are performed using a plane-wave pseudopotential formulation within the framework of density functional theory (DFT) [25] as implemented in the Vienna Ab initio Simulation Package (VASP) [26]. The generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) parameterization is adopted for the exchange-correlation functional [27]. Projector-augmented wave (PAW) potentials are used for Ga and N atoms and the planewave cutoff energy is set to 520 eV, with the energy convergence threshold set as 106 eV. In order to avoid layer interactions, a vacuum layer of 15 Å is added along c-axis for two-dimensional GaN structures. Moreover, the Brillouin zone is sampled with 9 9 5, 13 13 3 and 7 7 3 electronic wave-vector grid for Wurtzite-GaN, Hexagonal-GaN and Haeckelite-GaN, respectively. The electronic wave-vector grid and planewave cutoff energy are two important factors that can influence the convergence of total energy, and can further influence the force constants and thermal conductivity [28]. So we checked the convergence of the total energy during the calculation. It is shown that the currently used parameters for electronic wave-vector grid and planewave cutoff energy are sufficient to calculate the total energy. The 4 4 3, 5 5 1 and 3 3 1 supercells are used to calculate the second-order and third-order interatomic force constants (IFCs) for Wurtzite-GaN, HexagonalGaN, and Haeckelite-GaN, respectively. In this paper, the secondorder IFCs that are required to calculate the phonon dispersions are generated in the framework of density functional perturbation theory (DFPT) [29–31]. The third-order IFCs can be used to calculate the three phonon scattering rate by solving the linearized BTE numerically as implemented in the ShengBTE software package [24,32]. It should be noted that the calculated thermal conductivity in this study does not consider the fourth-phonon interactions, although Feng and Ruan have shown it have important effects on lattice thermal conductivity based on experiential potential [33]. At thermal equilibrium, in the absence of temperature gradient or other thermodynamical forces, the phonons are distributed obeying the Bose-Einstein statistics f0(xk ) [32]. The k comprises both a phonon branch index p and a wave vector q. The xk is the angular frequency. With a temperature gradient rT, the phonon distribution function f k deviates from f0 and this deviation can be obtained from the BTE. In the steady state, the phonon distribution is affected by diffusion and scattering, and it can be expressed by the BTE [32]:
@f k @f k @f þ k ¼0 ¼ @t @t diffusion @t scattering
ð1Þ
@f k @f ¼ rT v k k @t diffusion @T
ð2Þ
where v k is the phonon group velocity. The diffusion term is caused by the temperature gradient rT. The scattering term depends on
the specific scattering processes, which can be analyzed in the framework of perturbation theory due to intrinsic phonon-phonon interactions and impurities such as isotopes [32]. In most practical situations, rT is small enough that f k can be expanded to the first order in rT [32]:
f k ¼ f 0 ðxk Þ Fk rT
@f k @T
ð3Þ
Fk ¼ s0k ðv k þ Dk Þ
ð4Þ
Here, s0k is the relaxation time of mode k obtained from perturbation theory and Dk is the measure of the number a specific phonon mode with a dimension of a velocity [32]. Then, the Eq. (4) is solved by iterating from a zeroth-order approximation F0k ¼ s0k v k with relevant formulas for Dk and s0k as:
Dk ¼
þ 1X 1X 1 Cþkk0 k00 ðnkk00 Fk00 nkk0 Fk0 Þ þ C 0 00 ðn 00 F 00 þ nkk0 Fk0 Þ N k0 k00 N k0 k00 2 kk k kk k
ð5Þ 1
s0k
þ X 1 X 1 ¼ Cþkk0 k00 þ Ckk0 k00 N k0 k00 0 00 2 kk
! ð6Þ
where N ¼ N1 N 2 N 3 is the number of discretizing q points in the Brillouin zone, and nkk0 = xk0 /xk . Cþ kk0 k00 and Ckk0 k00 are the threephonon scattering rates corresponding to absorption processes (xk + xk0 = xk00 ) and emission processes (xk = xk0 + xk00 ), respectively [32]. Conservation of quasimomentum requires that q00 ¼ q q0 þ Q , Q is a reciprocal lattice vector and the signs to the two types of three-phonon process. The iterative process has a large impact on the results when calculating the high thermal conductivity materials such as diamond and graphene where the normal three-phonon processes play a significant role [32,34]. The lattice thermal conductivity j can be obtained by solving the phonon BTE written as [32]
j¼
1 kB T XN 2
X 2 0 0 ðhxk Þ f k ðf k þ 1Þv ak F bk
ð7Þ
k
where kB is the Boltzmann constant, ⁄ is the reduced Planck constant, T is absolute temperature, X is the volume of the unit cell. a and b are the Cartesian directions. F bk is the b component of the mean free displacement Fk, which can be obtained from phonon BTE [32]. We define lk ¼ Fk v k =jv k j as the scaler mean free path (MFP) for mode k. It should be noted that the volume X is defined for crystalline bulk materials which are three-dimensional systems. As to two-dimensional systems, a layer thickness d is required to redefine the thermal conductivity [35,36]. This value is selected to be d = 3.7 Å in this investigation, equivalent to the sum of van der Waals Radii of Ga (2.1 Å) and N (1.6 Å) atoms [37]. The effects of phonon wave-vector grid on the convergence of thermal conductivity are also checked. The evaluated error is about 1% when the phonon wave-vector grid is 15 15 15, 15 15 15 and 12 12 12 for Wurtzite-GaN, Hexagonal-GaN and HaeckeliteGaN, respectively. 3. Results and discussion 3.1. Anisotropic thermal conductivity of wurtzite GaN We firstly optimized the lattice structure of wurtzite GaN. It is found the obtained lattice constants are a = 3.247 Å and c = 5.283 Å, which are about 1.7% higher than those from experiment of a = 3.190 Å and c = 5.189 Å [38]. Then, we extract the second-order IFCs and calculate the phonon dispersions of wurtzite
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
GaN. Fig. 1 displays there are totally 12 phonon branches along any direction in the first Bullirium zone. The three phonon branches originated from C point are named acoustic phonon, and the others are optical phonon branches. There is a clear band-gap between 10.2 and 15.8 THz in the phonon spectrum, which agrees well with the results calculated by Xu et al. [39]. Below the bandgap, the acoustic and optical modes are overlapped far from the C point. After extracting the third-order IFCs and obtaining the phonon lifetimes, we calculate the lattice thermal conductivity of wurtzite GaN. Fig. 2 shows the temperature dependence of the lattice thermal conductivity of defect-free wurtzite GaN along the [0 0 1] and [1 0 0] directions. At room temperature, the lattice thermal conductivity along the [0 0 1] and [1 0 0] directions are 409 and 363 W/m-K, respectively. We define the anisotropy factor R in thermal conductivity by R = j[0 0 1]/j[1 0 0], where j[0 0 1] and j[1 0 0] are the thermal conductivity along [0 0 1] and [1 0 0] directions, respectively. Our calculation shows the R of defect-free GaN is about 1.13 at room temperature. The anisotropy R decreases continuously from 1.42 to 1.09 as the temperature increases from 100 to 800 K, indicating that the anisotropy of wurtzite GaN cannot be neglected, especially at low temperature. The similar anisotropy in thermal conductivity of GaN are also found in a recent MD simula-
25
Ga
Frequency [THz]
20
N
15
10
5
0
Γ
K M Γ A
H L
A
Wavevector Fig. 1. Structure and calculated phonon dispersion of wurtzite GaN.
2500 [001] direction [100] direction Lindsay et al. [001]_pure(Wu et al.) [100]_pure(Wu et al.) Mion et al.
κ [W/m-K]
2000 1500
Wurtzite-GaN
1000 500 0 0
200
400
600
800
Temperature [K] Fig. 2. The thermal conductivity j of wurtzite GaN as a function of temperature. The results from Lindsay et al. [10], Wu et al. [38], Mion et al. [39] are also shown.
421
tion by Ju et al. [12] and first-principles calculations by Wu et al. [40]. Wu et al. find that the anisotropy of thermal conductivity in wurtzite GaN increases with decreasing temperature, and j[1 0 0] and j[0 0 1] for pure wurtzite GaN are 361 W/m-K and 404 W/mK at 300 K. The stronger anisotropy in low temperature may imply that the phonon focusing effects [13] are an important contribution. Phonon focusing means that at some regions of the first Brillouin zone (FBZ), the phonon group velocity vectors point in nearly the same direction. For bulk GaN, the phonon focusing can result in a higher phonon irradiation heat flux along c-axis than that of inplane direction [12]. However, at high temperature the strong scattering process among phonons counteracts the phonon focusing effects, leading to the decreasing of R with increasing temperature. By considering the isotope scattering and impurity scattering, Lindsay et al. [10] shows the anisotropy of GaN is less than 1% at room temperature. We also found that the experimental results by Mion et al. [41] shown in Fig. 2 are lower than most theoretically predicted thermal conductivity, which may be resulted from the isotope and defect scattering in the samples.
3.2. Lattice thermal conductivity and phonon properties of monolayer GaN The single atomic layer of GaN is a newly emerging material, and the related phonon transport properties are less investigated. Here, we also calculate the phonon dispersions for the Hexagonal-GaN monolayer and Haeckelite-GaN monolayer in Fig. 3. The unit cell of Hexagonal-GaN monolayer only contains two atoms, while that of Haeckelite-GaN monolayer contains eight atoms. So the number of phonon branches in Hexagonal-GaN monolayer is much less than that in Haeckelite-GaN monolayer. Fig. 3 also shows that the transverse acoustic (TA) and longitudinal acoustic (LA) phonon dispersions of both Hexagonal-GaN and Haeckelite-GaN single layer are linear when the wave vector approaches the C point, while the out-of-plane acoustic (ZA) phonon dispersions of both two kinds of monolayer structures are quadratic. This phenomenon is similar to the phonon dispersions of graphene [42], which reflects a general characteristic of twodimensional materials. Comparing Figs. 1 and 3, it is found that the band gaps in the three GaN phases are significantly different. This difference may be explained from the atomic interaction strength. This is because the bond strength is directly correlated to the phonon frequencies. For example, the atomic distance in wurtzite and hexagonal monolayer are more or less the same, while the shortest atomic distance in haeckelite monolayer is shorter than that in wurtzite and hexagonal monolayer. Shorter atomic distance can increase the atomic interaction strength, and thus increase the phonon frequency and reduce the bandgap. Fig. 4 shows the intrinsic lattice thermal conductivity of both monolayers with temperature increasing from 100 to 800 K. As can be seen, the thermal conductivity of both Hexagonal-GaN and Haeckelite-GaN decreases with increasing temperature. The fitted curves show an approximate inverse relationship of
j / T 0:90 and j / T 0:88 for Hexagonal-GaN and Haeckelite-GaN, respectively. Fig. 4 shows that the thermal conductivity of Haeckelite-GaN is much smaller than those of Hexagonal-GaN from 100 to 800 K. For example, the thermal conductivity of Haeckelite-GaN is about 9.3 W/m-K at room temperature, while the corresponding value for Hexagonal-GaN is about 37 W/m-K, by a factor of about four. The difference of the thermal conductivity between the two structures is gradually reduced as the temperature increases. According to the phonon kinetic theory, the lattice thermal conductivity can be completely determined by the specific heat, group velocity and mean free path (MFP) of every phonon mode. We first
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
(a) Hexagonal-GaN Frequency [THz]
Ga
N
(b) Haeckelite-GaN
30 25
Ga
20
N
30 25
Frequency [THz]
422
15 10 5 0
20 15 10 5
Γ
M
0
Γ
K
Wavevector
Γ
X
Γ
M
Wavevector
Fig. 3. Structures and phonon dispersion of GaN monolayers (a) Hexagonal-GaN and (b) Haeckelite-GaN.
4x106
120 100
[W/m-K]
80
T-0.90
60
Wurtzite-GaN Hexagonal-GaN Haeckelite-GaN
Fitted Fitted
Cv [J/m3-K]
Hexagonal-GaN Haeckelite-GaN
3x106
2x106
40
T-0.88
20 0
0
200
1x106
400
600
0
800
200
Temperature [K]
calculate the specific heat cv of these two kind of atomic layer structure by
2
hx f 0 ðf 0 þ 1Þ k BT k Z hx 2 kB 3 f 0 ðf 0 þ 1Þd q: ¼ 3 kB T ð2pÞ
cv ¼
kB
600
800
XN
ð8Þ
Fig. 5 shows that the specific heat of both monolayer structures is much lower than that of bulk GaN. Moreover, the specific heat of Hexagonal-GaN is higher than that of Haeckelite-GaN. For example, the difference at 300 K is about 9.3%. Hence, the specific heat can partially explain the difference in the thermal conductivity between two monolayer structures. It should be noted that the calculated cv depends on the definition of thickness, and increasing the thickness can reduce the cv. The definition of thickness has the same effects on both cv and j. As long as the same thickness is defined for both Hexagonal-GaN and Haeckelite-GaN, the cv of Hexagonal-GaN should be larger than that of Haeckelite-GaN. Next we checked the phonon group velocity, which is equal to the derivation of phonon dispersion with respect to the wavevector. Fig. 6 shows that the group velocity of Hexagonal-GaN is generally larger than that of Haeckelite-GaN. This is because the atomic number density of the former (0.223 atoms/Å2) is larger than that of the latter (0.189 atoms/Å2). The higher atomic number density results in a higher phonon group velocity in Hexagonal-GaN.
Fig. 5. Specific heat cv of both bulk and monolayers as a function of temperature from 100 to 800 K.
12
Wurtzite-GaN Hexagonal-GaN Haeckelite-GaN
10
Group velocity [km/s]
Fig. 4. Lattice thermal conductivity j of GaN monolayers as a function of temperature from 100 to 800 K. The inverse relationship between lattice thermal conductivity and temperature is also shown in black (for Hexagonal GaN) and red (Haeckelite GaN) lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
X
400
Temparature [K]
T=300 K
8 6 4 2 0 0
5
10
15
20
25
30
Frequency [THz] Fig. 6. Phonon group velocity as a function of frequency for Wurtzite-GaN, Hexagonal-GaN and Haeckelite-GaN.
Fig. 6 also show the phonon group velocity of bulk GaN. Although the maximum on the low frequency for bulk GaN is a little lower than that of monolayer, it is obvious that the phonon group velocity around 15–20 THz is higher than that of both monolayer structures. Another important characteristic of the phonon that affects the lattice thermal conductivity is the phonon lifetimes [42,43]. Fig.7 shows the phonon lifetimes as a function of frequency for Wurtzite-GaN, Hexagonal-GaN and Haeckelite-GaN. It can be seen
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
smaller than that of single layer GaN. So this is also an evident why the thermal conductivity of monolayer GaN is lower than that of wurtzite GaN.
Phonon lifetimes [ps]
104 Wurtzite-GaN Hexagonal-GaN Haeckelite-GaN
102
T=300 K
3.3. Phonon modes contribution to the thermal conductivity
100
10-2
0
5
10
15
20
25
30
Frequency [THz] Fig. 7. Phonon lifetimes as a function of frequency for Wurtzite-GaN, HexagonalGaN and Haeckelite-GaN.
clearly that the phonon lifetimes of Haeckelite-GaN are smaller than those of Hexagonal-GaN, and the phonon lifetimes of single layer GaN are smaller than that of bulk GaN. With the calculated phonon group velocity and phonon lifetimes, we further compared the phonon MFP between Hexagonal-GaN and Haeckelite-GaN. Fig.8 shows that the phonon MFP of Haeckelite-GaN is remarkably shorter than that of Hexagonal-GaN, especially for the lowfrequency phonons around 5 THz. In addition, it should be noted that an inverse relationship exist between the phase space for three-phonon scattering P3 and the lattice thermal conductivity of a material [44]. The space phase P3 is a measure of the scattering space available for three-phonon processes. It represents the relative number of three-phonon scattering processes that satisfy the energy and momentum conservation selection rules., The threephonon phase space is calculated as [32]
P3 ¼
X
X ZZ
ð6pnÞ2 k;k0 ;k00
423
3
3
d½xk ðqÞ þ xk0 ðq0 Þ xk00 ðq þ q0 qÞd q0 d q ð9Þ
We calculated phase space P3 for Hexagonal-GaN (P3 = 3.38e3) and Haeckelite-GaN (P3 = 3.61e3). Therefore, our calculation clear shows a lower specific heat capacity, lower phonon group velocity and lower phonon MFP result in a lower thermal conductivity in Haeckelite-GaN as compared to that of Hexagonal-GaN. On the other hand, the P3 for wurtzite GaN is P3 = 1.39e3, which is much
Based on the obtained phonon group velocity and phonon lifetimes of every mode, we decompose the thermal conductivity into different mode contribution. Fig. 9 shows there are a lot of similarities about the modes contribution to the thermal conductivity between Hexagonal-GaN and Haeckelite-GaN. First, the TA modes dominate the thermal conductivity on the investigated temperature range. Second, the optical phonon modes have negligible contribution to the thermal conductivity. Third, the mode contribution to thermal conductivity is not sensitive to temperature above room temperature. But we also found that the mode relative contribution is different quantitatively between Hexagonal-GaN and Haeckelite-GaN. For Hexagonal-GaN, the contribution from TA modes is about 63%, which is about 3.5 times larger than that from ZA modes and 4.5 times larger than that from LA modes. But for Haeckelite-GaN, the contribution from TA modes is about 45%, which is 2 times larger than that from ZA modes and 3 times larger than that from LA modes. There is a significant difference between monolayer GaN and graphene in the contribution of phonon branches. For example, the contributions from ZA modes of Hexagonal-GaN (18%) and Haeckelite-GaN (24%) is much smaller than that of graphene (80%) [45]. Previous investigation shows that the large contribution from ZA mode to the lattice thermal conductivity of graphene is due to the reflection symmetry and the resulted selection rule, which forbid the scattering of ZA modes and increase their lifetimes [46]. However, some new twodimensional nanomaterials with the buckled structure, such as silicene [45], stanene [47] and phosphorene [48], have unique physical and chemical properties different from graphene and break the selection rule. According to our calculations, Hexagonal-GaN possesses a planar honeycomb structure, which is similar to graphene and monolayer BN [49]. However, because of the large differences in atomic radius between Ga (0.126 nm) and N (0.075 nm) atoms as compared to the graphene and monolayer BN, the planar honeycomb structure of Hexagonal-GaN is not as perfectly smooth as that of graphene and monolayer BN. Due to the difference in atom radius, the cohesive energy and in-plane stiffness of a honeycomb structure show a significant fluctuation [50]. We also compared the force constants of monolayer hexagonal GaN and graphene. It is found that the variation of third-order force constants of monolayer hexagonal GaN is different from that of graphene, which may means the symmetry selection rule is slightly broken in single GaN.
Wurtzite-GaN
104
Hexagonal-GaN
3.4. Phonon accumulation function
Phonon MFP [nm]
Haeckelite-GaN
T=300 K
102
100
10-2
10-4
0
5
10
15
20
25
30
Frequency [THz] Fig. 8. Phonon mean free path (MFP) as a function of frequency for Wurtzite-GaN, Hexagonal-GaN and Haeckelite-GaN.
Yang and Dames has shown the thermal conductivity of nanostructure can be obtained by the characteristic size of the nanostructure and the thermal conductivity accumulation function of the corresponding bulk materials [51]. Until now, the thermal conductivity accumulation as a function of phonon mean free path has been predicted [52,53] and experimentally measured [54,55], which can be conventionally used to tune the thermal conductivity of nanostructures. In Fig. 10, we present the cumulative lattice thermal conductiv~ of both monolayers as a function of phonon MFP at 300 K. It ity j can be seen that a broad phonon MFP contributes to the lattice thermal conductivity ranging from tens of nanometers to a few microns. About 60% of the heat is conducted by phonons with MFP between 10 and 100 nm, which is beneficial for us to under-
424
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
1.0
ZA TA LA optical
(a) Hexagonal-GaN
Mode contribution
0.8
1.0 0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0
200
400
600
800
ZA TA LA optical
(b) Haeckelite-GaN
0.0
200
Temperature [K]
400
600
800
Temperature [K]
Fig. 9. Phonon modes contribution to the lattice thermal conductivity of (a) Hexagonal-GaN and (b) Haeckelite-GaN.
50
Hexagonal-GaN Haeckelite-GaN
[W/mK]
40
1.0
Fitted Fitted
(a) Hexagonal-GaN
0.8
30
T=300K
300K 500K 800K
0.6
0.6
0.4
0.4
0.2
0 10-2
100
102
MFP [nm] ~ of both monolayers as a function Fig. 10. Cumulative lattice thermal conductivity j of the phonon MFP at 300 K.
stand how size affects thermal conductivity. Fitting curves of the ~ , according to a single cumulative lattice thermal conductivity j parametric function [32,56]
j~ ðl < lmax Þ ¼
jmax
1 þ l0 =lmax
;
0.0 100
102
MFP [nm]
104
ð10Þ
where jmax is the ultimate thermal conductivity, lmax is the maximal MFP concerned, and l0 is the parameter that can be interpreted as representative of the MFP of relevant heat-carrying phonons in monolayer GaN [32]. Fig. 10 shows that the calculated values l0 for monolayer GaN are 33.75 nm (Hexagonal-GaN) and 15.50 nm (Haeckelite-GaN), which is an estimate of the size at which the phonon boundary scattering becomes dominant over the three-phonon scattering. This is critical to the thermal design for modulating the thermal conductivity when the nanostructure size is below the characteristic size in the practical design of thermoelectric nanodevices. Furthermore, we investigate the normalized cumulative lattice ~ =jmax as an accumulation function of the thermal conductivity j phonon MFP at different temperatures. Fig. 11 shows that the normalized lattice thermal conductivity of Hexagonal-GaN and Haeckelite-GaN keeps increases as MFP increases, until reaching to the cutoff MFP of heat-carrying [57]. Their lattice thermal conductivity distribution curves look very similar, except for the different cutoff MFP of heat-carrying. At 300, 500 and 800 K, the
300K 500K 800K
0.2
0.0 10-2
(b) Haeckelite-GaN
0.8
20 10
1.0
104
10-2
100
102
104
MFP [nm]
~ =jmax as an accumulation Fig. 11. Normalized lattice thermal conductivity j function of the phonon MFP at 100 K, 300 K, 500 K and 800 K for (a) HexagonalGaN and (b) Haeckelite-GaN.
longest phonon MFP of heat-carrying are 0.34, 0.23 and 0.13 lm for Hexagonal-GaN as shown as in Fig. 11(a), and 1.79, 1.02 and 0. 59 lm for Haeckelite-GaN as shown in Fig. 11(b). Fig. 11 also shows that the cutoff MFP of heat-carrying is longer at lower temperature. This is because higher temperature increases the anharmonic interactions, and thus reduces the phonon MFP. By examining Fig. 11(a), it yields the parameter l0 of 33.75 nm, 18.80 nm and 10.48 nm for Hexagonal-GaN at 300, 500 and 800 K, respectively. The corresponding values for l0 of HaeckeliteGaN in Fig. 11(b) is 12.73, 7.09 and 3.70 nm lm from 300 to 800 K. It can be found that phonons with the MFP shorter than l0 contribute around 50% of the total thermal conductivity. The parameter l0 is helpful in studying size effect on ballistic or diffusive phonon transport. 4. Conclusions In summary, we systematically investigate the phonon properties of wurtzite GaN and the corresponding two-dimensional single layer structure from first-principles calculations and phonon BTE. It is found that the thermal conductivity along [0 0 1] direction is 13% higher than that along [1 0 0] direction at room temperature for defect-free wurtzite GaN. The thermal anisotropy is significantly reduced from 42% to 11% with increasing temperature from 100 to 800 K. We also found that the thermal conductivity of monolayer GaN with hexagonal structure is higher than that with
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
haeckelite structure by a factor of four. By comparing the specific heat, phonon group velocity and the phonon MFP between the two kinds of single layer structures, the underlying physics is uncovered to explain the thermal conductivity difference. Our calculation shows that the complex unit cell in Haeckelite-GaN single layer structure leads to a lower phonon group velocity and shorter phonon lifetimes. We further extract the thermal conductivity accumulation function with phonon MFP. It is found that the dominant phonon mean free path in Hexagonal monolayer is higher than that in Haeckelite. It is also found that the TA mode dominants the thermal conductivity for the GaN-single layer, which is not the same with that in graphene. Acknowledgements The authors thank Prof. H. Bao at Shanghai Jiao Tong University for useful discussion. The authors acknowledge the financial support from the Natural Science Foundation of China (51406034), and National Basic Research Program of China (2011CB707605). References [1] W. Chikhaoui, J.M. Bluet, M.A. Poisson, N. Sarazin, C. Dua, C. Bru-Chevallier, Appl. Phys. Lett. 96 (2010) 072107. [2] S. Ganguly, A. Konar, Z.Y. Hu, H.L. Xing, D. Jena, Appl. Phys. Lett. 101 (2012) 253519. [3] X.H. Wang, L.W. Guo, H.Q. Jia, Z.G. Xing, Y. Wang, X.J. Pei, J.M. Zhou, H. Chen, Appl. Phys. Lett. 94 (2009) 111913. [4] X. Li, Z.S. Liu, D.G. Zhao, D.S. Jiang, P. Chen, J.J. Zhu, J. Yang, L.C. Le, W. Liu, X.G. He, X.J. Li, F. Liang, L.Q. Zhang, J.P. Liu, H. Yang, Y.T. Zhang, G.T. Du, J. Vac. Sci. Technol. B 34 (2016) 041211. [5] B. Auld, Acoustic Fields and Waves in Solids, John Wiley & Sons, 1973. [6] A. Jezowski, B.A. Danilchenko, M. Bockowski, I. Grzegory, S. Krukowski, T. Suski, T. Paszkiewicz, Solid State Commun. 128 (2003) 69–73. [7] G.A. Slack, L.J. Schowalter, D. Morelli, J.A. Freitas, J. Cryst. Growth 246 (2002) 287–298. [8] M. Mukherjee, N. Mazumder, S.K. Roy, K. Goswami, Semicond. Sci. Technol. 22 (2007) 1258. [9] J. Zou, X. Lange, C. Richardson, J. Appl. Phys. 100 (2006) 104309. [10] L. Lindsay, D.A. Broido, T.L. Reinecke, Phys. Rev. Lett. 109 (2012) 095901. [11] X.F. Wu, J. Lee, V. Varshney, J.L. Wohlwend, A.K. Roy, T.F. Luo, Sci. Rep-Uk 6 (2016) 22504. [12] W.J. Ju, Z.Y. Zhou, Z.Y. Wei, Aip. Adv. 6 (2016) 065328. [13] Z.Y. Wei, Y.F. Chen, C. Dames, Appl. Phys. Lett. 102 (2013) 011901. [14] K.S. Novoselov, D. Jiang, F. Schedin, T.J. Booth, V.V. Khotkevich, S.V. Morozov, A. K. Geim, Proc. Natl. Acad. Sci. U.S.A. 102 (2005) 10451–10453. [15] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Science 306 (2004) 666–669. [16] M.X. Liu, Y.C. Li, P.C. Chen, J.Y. Sun, D.L. Ma, Q.C. Li, T. Gao, Y.B. Gao, Z.H. Cheng, X.H. Qiu, Y. Fang, Y.F. Zhang, Z.F. Liu, Nano Lett. 14 (2014) 6342–6347.
425
[17] A. Splendiani, L. Sun, Y.B. Zhang, T.S. Li, J. Kim, C.Y. Chim, G. Galli, F. Wang, Nano Lett. 10 (2010) 1271–1275. [18] L.K. Li, Y.J. Yu, G.J. Ye, Q.Q. Ge, X.D. Ou, H. Wu, D.L. Feng, X.H. Chen, Y.B. Zhang, Nat. Nanotechnol. 9 (2014) 372–377. [19] Z.Y. Al Balushi, K. Wang, R.K. Ghosh, R.A. Vila, S.M. Eichfeld, J.D. Caldwell, X.Y. Qin, Y.C. Lin, P.A. DeSario, G. Stone, S. Subramanian, D.F. Paul, R.M. Wallace, S. Datta, J.M. Redwing, J.A. Robinson, Nat. Mater. 15 (2016) 1166–1171. [20] N.A. Koratkar, Nat. Mater. 15 (2016) 1153–1154. [21] J.V.N. Sarma, R. Chowdhury, R. Jayaganthan, J. Appl. Phys. 113 (2013) 243504. [22] D.C. Camacho-Mojica, F. Lopez-Urias, Sci. Rep-Uk 5 (2015) 17902. [23] D.A. Broido, M. Malorny, G. Birner, N. Mingo, D.A. Stewart, Appl. Phys. Lett. 91 (2007) 231922. [24] M. Omini, A. Sparavigna, Physica B 212 (1995) 101–112. [25] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [26] G. Kresse, J. Furthmuller, Phys. Rev. B 54 (1996) 11169–11186. [27] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [28] A. Jain, A.J.H. McGaughey, Comput. Mater. Sci 110 (2015) 115–120. [29] S. Baroni, S. de Gironcoli, A. Dal Corso, P. Giannozzi, Rev. Mod. Phys. 73 (2001) 515–562. [30] G. Deinzer, G. Birner, D. Strauch, Phys. Rev. B 67 (2003) 144304. [31] L. Lindsay, D.A. Broido, T.L. Reinecke, Phys. Rev. B 87 (2013) 165201. [32] W. Li, J. Carrete, N.A. Katcho, N. Mingo, Comput. Phys. Commun. 185 (2014) 1747–1758. [33] T.L. Feng, X.L. Ruan, Phys. Rev. B 93 (2016) 045202. [34] A. Ward, D.A. Broido, D.A. Stewart, G. Deinzer, Phys. Rev. B 80 (2009). [35] X.J. Tan, H.Z. Shao, T.Q. Hu, G.Q. Liu, J. Jiang, H.C. Jiang, PCCP 17 (2015) 22872– 22881. [36] L.L. Zhao, S. Xu, M.C. Wang, S.C. Lin, J. Phys. Chem. C 120 (2016) 27675–27681. [37] S.S. Batsanov, Inorg. Mater. 37 (2001) 871–885. [38] H. Schulz, K. Thiemann, Solid State Commun. 23 (1977) 815–819. [39] L.C. Xu, R.Z. Wang, X.D. Yang, H. Yan, J. Appl. Phys. 110 (2011) 043528. [40] R.K. Wu, R. Hu, X.B. Luo, J. Appl. Phys. 119 (2016) 145706. [41] C. Mion, J.F. Muth, E.A. Preble, D. Hanser, Appl. Phys. Lett. 89 (2006) 092123. [42] N. Bonini, J. Garg, N. Marzari, Nano Lett. 12 (2012) 2673–2678. [43] A. Togo, L. Chaput, I. Tanaka, Phys. Rev. B 91 (2015) 094306. [44] L. Lindsay, D.A. Broido, J. Phys. – Condens. Matter 20 (2008) 165209. [45] X.K. Gu, R.G. Yang, J. Appl. Phys. 117 (2015) 025102. [46] L. Lindsay, D.A. Broido, N. Mingo, Phys. Rev. B 82 (2010) 115427. [47] B. Peng, H. Zhang, H.Z. Shao, Y.C. Xu, X.C. Zhang, H.Y. Zhu, Sci. Rep-Uk 6 (2016) 20225. [48] G.Z. Qin, Q.B. Yan, Z.Z. Qin, S.Y. Yue, M. Hu, G. Su, PCCP 17 (2015) 4854–4858. [49] L. Lindsay, D.A. Broido, Phys. Rev. B 84 (2011) 155421. [50] H. Sahin, S. Cahangirov, M. Topsakal, E. Bekaroglu, E. Akturk, R.T. Senger, S. Ciraci, Phys. Rev. B 80 (2009) 155453. [51] F. Yang, C. Dames, Phys. Rev. B 87 (2013) 035437. [52] Z.Y. Wei, J.K. Yang, W.Y. Chen, K.D. Bi, D.Y. Li, Y.F. Chen, Appl. Phys. Lett. 104 (2014) 081903. [53] A.S. Henry, G. Chen, J. Comput. Theor. Nanos. 5 (2008) 141–152. [54] A.J. Minnich, J.A. Johnson, A.J. Schmidt, K. Esfarjani, M.S. Dresselhaus, K.A. Nelson, G. Chen, Phys. Rev. Lett. 107 (2011) 095901. [55] K.T. Regner, D.P. Sellan, Z.H. Su, C.H. Amon, A.J.H. McGaughey, J.A. Malen, Nat. Commun. 4 (2013) 1640. [56] L. Lindsay, D.A. Broido, N. Mingo, Phys. Rev. B 83 (2011) 235428. [57] G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, F. Mauri, Nano Lett. 14 (2014) 6109–6114.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Phonon transport properties of bulk and monolayer GaN from first-principles calculations Yongqiang Jiang, Shuang Cai, Yi Tao, Zhiyong Wei ⇑, Kedong Bi, Yunfei Chen Jiangsu Key Laboratory for Design & Manufacture of Micro/Nano Biomedical Instruments and School of Mechanical Engineering, Southeast University, Nanjing 211189, People’s Republic of China
a r t i c l e
i n f o
Article history: Received 24 February 2017 Received in revised form 27 June 2017 Accepted 11 July 2017
Keywords: Phonon Thermal conductivity GaN First-principles Mean free path
a b s t r a c t The thermal conductivity of isotopically pure wurtzite GaN and the corresponding two-dimensional monolayer crystals are investigated based on first-principles calculations and phonon Boltzmann transport equation. It is found that the c-axis thermal conductivity of defect-free wurtzite GaN is 42–11% higher than that along the in-plane direction when temperature increases from 100 to 800 K. Our calculation also shows that the thermal conductivity of monolayer GaN with hexagonal lattice is higher than that with haeckelite lattice by a factor of four at room temperature. The phonon dispersion relation, phonon group velocity and the phonon scattering rate are extracted to uncover the underlying physical mechanism. The thermal conductivity accumulation function with respect to the phonon mean free path and the relative contribution of each phonon branch for the two kinds of monolayer GaN are further calculated as compared to that of graphene. These investigations would provide important guidance to develop GaN-based electronic and photonic device. Ó 2017 Elsevier B.V. All rights reserved.
1. Introduction Gallium nitride (GaN) is an important wide band-gap semiconductor that can be widely used in electronics and optoelectronics devices, such as microwave power transistors [1,2], lightemitting diodes [3] and laser diodes [4]. The most stable GaN in nature (wurtzite GaN) belongs to hexagonal crystal system, which is the most anisotropic among the seven kinds of crystal family [5]. However, few investigations are conducted to systematically study its anisotropic thermal conductivity. Most previous experiments results [6–8] showed that the thermal conductivity of wurtzite GaN is isotropic around room temperature due to its nearly isotropic elastic property. On the other hand, most theoretical works [9–11] often take the average value of the three directions or the value of one arbitrary direction as the thermal conductivity of GaN at 300 K, and omit its anisotropy. By performing nonequilibrium molecular dynamics simulation, recently Ju et al. [12] shows the thermal conductivity along the c-axis is 25% higher than that along the in-plane direction at room temperature. Although the phonon focusing effects [13] are used to explain the anisotropy, the anisotropic thermal conductivity in wurtzite GaN is still needed to be verified.
⇑ Corresponding author. E-mail address: [email protected] (Z. Wei). http://dx.doi.org/10.1016/j.commatsci.2017.07.012 0927-0256/Ó 2017 Elsevier B.V. All rights reserved.
In the past decade, the exfoliated single layer of graphene from bulk graphite by Novoselov and Geim [14,15] has attracted much attention for the study of atomic layer structure, such as BN, MoS2, black phosphorus, etc. [16–18]. Since there are no weak van der Waals interactions in bulk GaN, it is difficult to exfoliate single layer GaN from its corresponding bulk structure. Thus, little attention is focused on the physical properties of monolayer GaN. However, recent experiment has shown that the atomic layer GaN can be obtained by a migration-enhanced encapsulation growth (MEEG) technique [19,20]. This experimentally stable monolayer GaN is similar to graphene, which is verified by the scanning electron microscope (SEM) and scanning transmission electron microscope (STEM) [19,20]. Sarma et al. [21] shows that the thermal conductivity of this monolayer hexagonal GaN is about 25 W/m-K at room temperature by molecular dynamics simulations. Very recently a new allotrope of monolayer GaN consisting of square and octagonal rings (called Haeckelite 8–4) was proposed and Dulce et al. [22] demonstrated that both Hexagonal (Hexagonal-GaN) and Haeckelite (Haeckelite-GaN) structures are thermally stable. Understanding the phonon transport is crucial to predict the thermal performance of newly emerging materials. According to the lattice difference between the two kinds of monolayer GaN, it is expected that their thermal conductivity may be different. The electronic properties of single layer GaN has been extensively investigated recently, however, the thermal transport properties of them are not systematically investigated so far.
420
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
In this study, we systematically investigate the phonon transport properties of wurtzite GaN crystal and its lattice thermal conductivity by combining first-principles calculations [23] with phonon Boltzmann transport equation (BTE) [24]. The calculated results show that the lattice thermal conductivity of defect-free wurtzite GaN exhibits obvious anisotropic behavior at low temperature. At room temperature, the thermal conductivity of defectfree wurtzite GaN in the [0 0 1] direction is 13% higher than that in the [1 0 0] direction. We also study the thermal behavior of two monolayer GaN lattices by the same method. The phonon properties are also analyzed comprehensively based on phonon dispersion relation, phonon group velocity, phonon lifetimes and specific heat to reveal differences in the thermal transport properties of two kinds of atomic layer. 2. Computational details The first-principles calculations of GaN are performed using a plane-wave pseudopotential formulation within the framework of density functional theory (DFT) [25] as implemented in the Vienna Ab initio Simulation Package (VASP) [26]. The generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) parameterization is adopted for the exchange-correlation functional [27]. Projector-augmented wave (PAW) potentials are used for Ga and N atoms and the planewave cutoff energy is set to 520 eV, with the energy convergence threshold set as 106 eV. In order to avoid layer interactions, a vacuum layer of 15 Å is added along c-axis for two-dimensional GaN structures. Moreover, the Brillouin zone is sampled with 9 9 5, 13 13 3 and 7 7 3 electronic wave-vector grid for Wurtzite-GaN, Hexagonal-GaN and Haeckelite-GaN, respectively. The electronic wave-vector grid and planewave cutoff energy are two important factors that can influence the convergence of total energy, and can further influence the force constants and thermal conductivity [28]. So we checked the convergence of the total energy during the calculation. It is shown that the currently used parameters for electronic wave-vector grid and planewave cutoff energy are sufficient to calculate the total energy. The 4 4 3, 5 5 1 and 3 3 1 supercells are used to calculate the second-order and third-order interatomic force constants (IFCs) for Wurtzite-GaN, HexagonalGaN, and Haeckelite-GaN, respectively. In this paper, the secondorder IFCs that are required to calculate the phonon dispersions are generated in the framework of density functional perturbation theory (DFPT) [29–31]. The third-order IFCs can be used to calculate the three phonon scattering rate by solving the linearized BTE numerically as implemented in the ShengBTE software package [24,32]. It should be noted that the calculated thermal conductivity in this study does not consider the fourth-phonon interactions, although Feng and Ruan have shown it have important effects on lattice thermal conductivity based on experiential potential [33]. At thermal equilibrium, in the absence of temperature gradient or other thermodynamical forces, the phonons are distributed obeying the Bose-Einstein statistics f0(xk ) [32]. The k comprises both a phonon branch index p and a wave vector q. The xk is the angular frequency. With a temperature gradient rT, the phonon distribution function f k deviates from f0 and this deviation can be obtained from the BTE. In the steady state, the phonon distribution is affected by diffusion and scattering, and it can be expressed by the BTE [32]:
@f k @f k @f þ k ¼0 ¼ @t @t diffusion @t scattering
ð1Þ
@f k @f ¼ rT v k k @t diffusion @T
ð2Þ
where v k is the phonon group velocity. The diffusion term is caused by the temperature gradient rT. The scattering term depends on
the specific scattering processes, which can be analyzed in the framework of perturbation theory due to intrinsic phonon-phonon interactions and impurities such as isotopes [32]. In most practical situations, rT is small enough that f k can be expanded to the first order in rT [32]:
f k ¼ f 0 ðxk Þ Fk rT
@f k @T
ð3Þ
Fk ¼ s0k ðv k þ Dk Þ
ð4Þ
Here, s0k is the relaxation time of mode k obtained from perturbation theory and Dk is the measure of the number a specific phonon mode with a dimension of a velocity [32]. Then, the Eq. (4) is solved by iterating from a zeroth-order approximation F0k ¼ s0k v k with relevant formulas for Dk and s0k as:
Dk ¼
þ 1X 1X 1 Cþkk0 k00 ðnkk00 Fk00 nkk0 Fk0 Þ þ C 0 00 ðn 00 F 00 þ nkk0 Fk0 Þ N k0 k00 N k0 k00 2 kk k kk k
ð5Þ 1
s0k
þ X 1 X 1 ¼ Cþkk0 k00 þ Ckk0 k00 N k0 k00 0 00 2 kk
! ð6Þ
where N ¼ N1 N 2 N 3 is the number of discretizing q points in the Brillouin zone, and nkk0 = xk0 /xk . Cþ kk0 k00 and Ckk0 k00 are the threephonon scattering rates corresponding to absorption processes (xk + xk0 = xk00 ) and emission processes (xk = xk0 + xk00 ), respectively [32]. Conservation of quasimomentum requires that q00 ¼ q q0 þ Q , Q is a reciprocal lattice vector and the signs to the two types of three-phonon process. The iterative process has a large impact on the results when calculating the high thermal conductivity materials such as diamond and graphene where the normal three-phonon processes play a significant role [32,34]. The lattice thermal conductivity j can be obtained by solving the phonon BTE written as [32]
j¼
1 kB T XN 2
X 2 0 0 ðhxk Þ f k ðf k þ 1Þv ak F bk
ð7Þ
k
where kB is the Boltzmann constant, ⁄ is the reduced Planck constant, T is absolute temperature, X is the volume of the unit cell. a and b are the Cartesian directions. F bk is the b component of the mean free displacement Fk, which can be obtained from phonon BTE [32]. We define lk ¼ Fk v k =jv k j as the scaler mean free path (MFP) for mode k. It should be noted that the volume X is defined for crystalline bulk materials which are three-dimensional systems. As to two-dimensional systems, a layer thickness d is required to redefine the thermal conductivity [35,36]. This value is selected to be d = 3.7 Å in this investigation, equivalent to the sum of van der Waals Radii of Ga (2.1 Å) and N (1.6 Å) atoms [37]. The effects of phonon wave-vector grid on the convergence of thermal conductivity are also checked. The evaluated error is about 1% when the phonon wave-vector grid is 15 15 15, 15 15 15 and 12 12 12 for Wurtzite-GaN, Hexagonal-GaN and HaeckeliteGaN, respectively. 3. Results and discussion 3.1. Anisotropic thermal conductivity of wurtzite GaN We firstly optimized the lattice structure of wurtzite GaN. It is found the obtained lattice constants are a = 3.247 Å and c = 5.283 Å, which are about 1.7% higher than those from experiment of a = 3.190 Å and c = 5.189 Å [38]. Then, we extract the second-order IFCs and calculate the phonon dispersions of wurtzite
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
GaN. Fig. 1 displays there are totally 12 phonon branches along any direction in the first Bullirium zone. The three phonon branches originated from C point are named acoustic phonon, and the others are optical phonon branches. There is a clear band-gap between 10.2 and 15.8 THz in the phonon spectrum, which agrees well with the results calculated by Xu et al. [39]. Below the bandgap, the acoustic and optical modes are overlapped far from the C point. After extracting the third-order IFCs and obtaining the phonon lifetimes, we calculate the lattice thermal conductivity of wurtzite GaN. Fig. 2 shows the temperature dependence of the lattice thermal conductivity of defect-free wurtzite GaN along the [0 0 1] and [1 0 0] directions. At room temperature, the lattice thermal conductivity along the [0 0 1] and [1 0 0] directions are 409 and 363 W/m-K, respectively. We define the anisotropy factor R in thermal conductivity by R = j[0 0 1]/j[1 0 0], where j[0 0 1] and j[1 0 0] are the thermal conductivity along [0 0 1] and [1 0 0] directions, respectively. Our calculation shows the R of defect-free GaN is about 1.13 at room temperature. The anisotropy R decreases continuously from 1.42 to 1.09 as the temperature increases from 100 to 800 K, indicating that the anisotropy of wurtzite GaN cannot be neglected, especially at low temperature. The similar anisotropy in thermal conductivity of GaN are also found in a recent MD simula-
25
Ga
Frequency [THz]
20
N
15
10
5
0
Γ
K M Γ A
H L
A
Wavevector Fig. 1. Structure and calculated phonon dispersion of wurtzite GaN.
2500 [001] direction [100] direction Lindsay et al. [001]_pure(Wu et al.) [100]_pure(Wu et al.) Mion et al.
κ [W/m-K]
2000 1500
Wurtzite-GaN
1000 500 0 0
200
400
600
800
Temperature [K] Fig. 2. The thermal conductivity j of wurtzite GaN as a function of temperature. The results from Lindsay et al. [10], Wu et al. [38], Mion et al. [39] are also shown.
421
tion by Ju et al. [12] and first-principles calculations by Wu et al. [40]. Wu et al. find that the anisotropy of thermal conductivity in wurtzite GaN increases with decreasing temperature, and j[1 0 0] and j[0 0 1] for pure wurtzite GaN are 361 W/m-K and 404 W/mK at 300 K. The stronger anisotropy in low temperature may imply that the phonon focusing effects [13] are an important contribution. Phonon focusing means that at some regions of the first Brillouin zone (FBZ), the phonon group velocity vectors point in nearly the same direction. For bulk GaN, the phonon focusing can result in a higher phonon irradiation heat flux along c-axis than that of inplane direction [12]. However, at high temperature the strong scattering process among phonons counteracts the phonon focusing effects, leading to the decreasing of R with increasing temperature. By considering the isotope scattering and impurity scattering, Lindsay et al. [10] shows the anisotropy of GaN is less than 1% at room temperature. We also found that the experimental results by Mion et al. [41] shown in Fig. 2 are lower than most theoretically predicted thermal conductivity, which may be resulted from the isotope and defect scattering in the samples.
3.2. Lattice thermal conductivity and phonon properties of monolayer GaN The single atomic layer of GaN is a newly emerging material, and the related phonon transport properties are less investigated. Here, we also calculate the phonon dispersions for the Hexagonal-GaN monolayer and Haeckelite-GaN monolayer in Fig. 3. The unit cell of Hexagonal-GaN monolayer only contains two atoms, while that of Haeckelite-GaN monolayer contains eight atoms. So the number of phonon branches in Hexagonal-GaN monolayer is much less than that in Haeckelite-GaN monolayer. Fig. 3 also shows that the transverse acoustic (TA) and longitudinal acoustic (LA) phonon dispersions of both Hexagonal-GaN and Haeckelite-GaN single layer are linear when the wave vector approaches the C point, while the out-of-plane acoustic (ZA) phonon dispersions of both two kinds of monolayer structures are quadratic. This phenomenon is similar to the phonon dispersions of graphene [42], which reflects a general characteristic of twodimensional materials. Comparing Figs. 1 and 3, it is found that the band gaps in the three GaN phases are significantly different. This difference may be explained from the atomic interaction strength. This is because the bond strength is directly correlated to the phonon frequencies. For example, the atomic distance in wurtzite and hexagonal monolayer are more or less the same, while the shortest atomic distance in haeckelite monolayer is shorter than that in wurtzite and hexagonal monolayer. Shorter atomic distance can increase the atomic interaction strength, and thus increase the phonon frequency and reduce the bandgap. Fig. 4 shows the intrinsic lattice thermal conductivity of both monolayers with temperature increasing from 100 to 800 K. As can be seen, the thermal conductivity of both Hexagonal-GaN and Haeckelite-GaN decreases with increasing temperature. The fitted curves show an approximate inverse relationship of
j / T 0:90 and j / T 0:88 for Hexagonal-GaN and Haeckelite-GaN, respectively. Fig. 4 shows that the thermal conductivity of Haeckelite-GaN is much smaller than those of Hexagonal-GaN from 100 to 800 K. For example, the thermal conductivity of Haeckelite-GaN is about 9.3 W/m-K at room temperature, while the corresponding value for Hexagonal-GaN is about 37 W/m-K, by a factor of about four. The difference of the thermal conductivity between the two structures is gradually reduced as the temperature increases. According to the phonon kinetic theory, the lattice thermal conductivity can be completely determined by the specific heat, group velocity and mean free path (MFP) of every phonon mode. We first
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
(a) Hexagonal-GaN Frequency [THz]
Ga
N
(b) Haeckelite-GaN
30 25
Ga
20
N
30 25
Frequency [THz]
422
15 10 5 0
20 15 10 5
Γ
M
0
Γ
K
Wavevector
Γ
X
Γ
M
Wavevector
Fig. 3. Structures and phonon dispersion of GaN monolayers (a) Hexagonal-GaN and (b) Haeckelite-GaN.
4x106
120 100
[W/m-K]
80
T-0.90
60
Wurtzite-GaN Hexagonal-GaN Haeckelite-GaN
Fitted Fitted
Cv [J/m3-K]
Hexagonal-GaN Haeckelite-GaN
3x106
2x106
40
T-0.88
20 0
0
200
1x106
400
600
0
800
200
Temperature [K]
calculate the specific heat cv of these two kind of atomic layer structure by
2
hx f 0 ðf 0 þ 1Þ k BT k Z hx 2 kB 3 f 0 ðf 0 þ 1Þd q: ¼ 3 kB T ð2pÞ
cv ¼
kB
600
800
XN
ð8Þ
Fig. 5 shows that the specific heat of both monolayer structures is much lower than that of bulk GaN. Moreover, the specific heat of Hexagonal-GaN is higher than that of Haeckelite-GaN. For example, the difference at 300 K is about 9.3%. Hence, the specific heat can partially explain the difference in the thermal conductivity between two monolayer structures. It should be noted that the calculated cv depends on the definition of thickness, and increasing the thickness can reduce the cv. The definition of thickness has the same effects on both cv and j. As long as the same thickness is defined for both Hexagonal-GaN and Haeckelite-GaN, the cv of Hexagonal-GaN should be larger than that of Haeckelite-GaN. Next we checked the phonon group velocity, which is equal to the derivation of phonon dispersion with respect to the wavevector. Fig. 6 shows that the group velocity of Hexagonal-GaN is generally larger than that of Haeckelite-GaN. This is because the atomic number density of the former (0.223 atoms/Å2) is larger than that of the latter (0.189 atoms/Å2). The higher atomic number density results in a higher phonon group velocity in Hexagonal-GaN.
Fig. 5. Specific heat cv of both bulk and monolayers as a function of temperature from 100 to 800 K.
12
Wurtzite-GaN Hexagonal-GaN Haeckelite-GaN
10
Group velocity [km/s]
Fig. 4. Lattice thermal conductivity j of GaN monolayers as a function of temperature from 100 to 800 K. The inverse relationship between lattice thermal conductivity and temperature is also shown in black (for Hexagonal GaN) and red (Haeckelite GaN) lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
X
400
Temparature [K]
T=300 K
8 6 4 2 0 0
5
10
15
20
25
30
Frequency [THz] Fig. 6. Phonon group velocity as a function of frequency for Wurtzite-GaN, Hexagonal-GaN and Haeckelite-GaN.
Fig. 6 also show the phonon group velocity of bulk GaN. Although the maximum on the low frequency for bulk GaN is a little lower than that of monolayer, it is obvious that the phonon group velocity around 15–20 THz is higher than that of both monolayer structures. Another important characteristic of the phonon that affects the lattice thermal conductivity is the phonon lifetimes [42,43]. Fig.7 shows the phonon lifetimes as a function of frequency for Wurtzite-GaN, Hexagonal-GaN and Haeckelite-GaN. It can be seen
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
smaller than that of single layer GaN. So this is also an evident why the thermal conductivity of monolayer GaN is lower than that of wurtzite GaN.
Phonon lifetimes [ps]
104 Wurtzite-GaN Hexagonal-GaN Haeckelite-GaN
102
T=300 K
3.3. Phonon modes contribution to the thermal conductivity
100
10-2
0
5
10
15
20
25
30
Frequency [THz] Fig. 7. Phonon lifetimes as a function of frequency for Wurtzite-GaN, HexagonalGaN and Haeckelite-GaN.
clearly that the phonon lifetimes of Haeckelite-GaN are smaller than those of Hexagonal-GaN, and the phonon lifetimes of single layer GaN are smaller than that of bulk GaN. With the calculated phonon group velocity and phonon lifetimes, we further compared the phonon MFP between Hexagonal-GaN and Haeckelite-GaN. Fig.8 shows that the phonon MFP of Haeckelite-GaN is remarkably shorter than that of Hexagonal-GaN, especially for the lowfrequency phonons around 5 THz. In addition, it should be noted that an inverse relationship exist between the phase space for three-phonon scattering P3 and the lattice thermal conductivity of a material [44]. The space phase P3 is a measure of the scattering space available for three-phonon processes. It represents the relative number of three-phonon scattering processes that satisfy the energy and momentum conservation selection rules., The threephonon phase space is calculated as [32]
P3 ¼
X
X ZZ
ð6pnÞ2 k;k0 ;k00
423
3
3
d½xk ðqÞ þ xk0 ðq0 Þ xk00 ðq þ q0 qÞd q0 d q ð9Þ
We calculated phase space P3 for Hexagonal-GaN (P3 = 3.38e3) and Haeckelite-GaN (P3 = 3.61e3). Therefore, our calculation clear shows a lower specific heat capacity, lower phonon group velocity and lower phonon MFP result in a lower thermal conductivity in Haeckelite-GaN as compared to that of Hexagonal-GaN. On the other hand, the P3 for wurtzite GaN is P3 = 1.39e3, which is much
Based on the obtained phonon group velocity and phonon lifetimes of every mode, we decompose the thermal conductivity into different mode contribution. Fig. 9 shows there are a lot of similarities about the modes contribution to the thermal conductivity between Hexagonal-GaN and Haeckelite-GaN. First, the TA modes dominate the thermal conductivity on the investigated temperature range. Second, the optical phonon modes have negligible contribution to the thermal conductivity. Third, the mode contribution to thermal conductivity is not sensitive to temperature above room temperature. But we also found that the mode relative contribution is different quantitatively between Hexagonal-GaN and Haeckelite-GaN. For Hexagonal-GaN, the contribution from TA modes is about 63%, which is about 3.5 times larger than that from ZA modes and 4.5 times larger than that from LA modes. But for Haeckelite-GaN, the contribution from TA modes is about 45%, which is 2 times larger than that from ZA modes and 3 times larger than that from LA modes. There is a significant difference between monolayer GaN and graphene in the contribution of phonon branches. For example, the contributions from ZA modes of Hexagonal-GaN (18%) and Haeckelite-GaN (24%) is much smaller than that of graphene (80%) [45]. Previous investigation shows that the large contribution from ZA mode to the lattice thermal conductivity of graphene is due to the reflection symmetry and the resulted selection rule, which forbid the scattering of ZA modes and increase their lifetimes [46]. However, some new twodimensional nanomaterials with the buckled structure, such as silicene [45], stanene [47] and phosphorene [48], have unique physical and chemical properties different from graphene and break the selection rule. According to our calculations, Hexagonal-GaN possesses a planar honeycomb structure, which is similar to graphene and monolayer BN [49]. However, because of the large differences in atomic radius between Ga (0.126 nm) and N (0.075 nm) atoms as compared to the graphene and monolayer BN, the planar honeycomb structure of Hexagonal-GaN is not as perfectly smooth as that of graphene and monolayer BN. Due to the difference in atom radius, the cohesive energy and in-plane stiffness of a honeycomb structure show a significant fluctuation [50]. We also compared the force constants of monolayer hexagonal GaN and graphene. It is found that the variation of third-order force constants of monolayer hexagonal GaN is different from that of graphene, which may means the symmetry selection rule is slightly broken in single GaN.
Wurtzite-GaN
104
Hexagonal-GaN
3.4. Phonon accumulation function
Phonon MFP [nm]
Haeckelite-GaN
T=300 K
102
100
10-2
10-4
0
5
10
15
20
25
30
Frequency [THz] Fig. 8. Phonon mean free path (MFP) as a function of frequency for Wurtzite-GaN, Hexagonal-GaN and Haeckelite-GaN.
Yang and Dames has shown the thermal conductivity of nanostructure can be obtained by the characteristic size of the nanostructure and the thermal conductivity accumulation function of the corresponding bulk materials [51]. Until now, the thermal conductivity accumulation as a function of phonon mean free path has been predicted [52,53] and experimentally measured [54,55], which can be conventionally used to tune the thermal conductivity of nanostructures. In Fig. 10, we present the cumulative lattice thermal conductiv~ of both monolayers as a function of phonon MFP at 300 K. It ity j can be seen that a broad phonon MFP contributes to the lattice thermal conductivity ranging from tens of nanometers to a few microns. About 60% of the heat is conducted by phonons with MFP between 10 and 100 nm, which is beneficial for us to under-
424
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
1.0
ZA TA LA optical
(a) Hexagonal-GaN
Mode contribution
0.8
1.0 0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0
200
400
600
800
ZA TA LA optical
(b) Haeckelite-GaN
0.0
200
Temperature [K]
400
600
800
Temperature [K]
Fig. 9. Phonon modes contribution to the lattice thermal conductivity of (a) Hexagonal-GaN and (b) Haeckelite-GaN.
50
Hexagonal-GaN Haeckelite-GaN
[W/mK]
40
1.0
Fitted Fitted
(a) Hexagonal-GaN
0.8
30
T=300K
300K 500K 800K
0.6
0.6
0.4
0.4
0.2
0 10-2
100
102
MFP [nm] ~ of both monolayers as a function Fig. 10. Cumulative lattice thermal conductivity j of the phonon MFP at 300 K.
stand how size affects thermal conductivity. Fitting curves of the ~ , according to a single cumulative lattice thermal conductivity j parametric function [32,56]
j~ ðl < lmax Þ ¼
jmax
1 þ l0 =lmax
;
0.0 100
102
MFP [nm]
104
ð10Þ
where jmax is the ultimate thermal conductivity, lmax is the maximal MFP concerned, and l0 is the parameter that can be interpreted as representative of the MFP of relevant heat-carrying phonons in monolayer GaN [32]. Fig. 10 shows that the calculated values l0 for monolayer GaN are 33.75 nm (Hexagonal-GaN) and 15.50 nm (Haeckelite-GaN), which is an estimate of the size at which the phonon boundary scattering becomes dominant over the three-phonon scattering. This is critical to the thermal design for modulating the thermal conductivity when the nanostructure size is below the characteristic size in the practical design of thermoelectric nanodevices. Furthermore, we investigate the normalized cumulative lattice ~ =jmax as an accumulation function of the thermal conductivity j phonon MFP at different temperatures. Fig. 11 shows that the normalized lattice thermal conductivity of Hexagonal-GaN and Haeckelite-GaN keeps increases as MFP increases, until reaching to the cutoff MFP of heat-carrying [57]. Their lattice thermal conductivity distribution curves look very similar, except for the different cutoff MFP of heat-carrying. At 300, 500 and 800 K, the
300K 500K 800K
0.2
0.0 10-2
(b) Haeckelite-GaN
0.8
20 10
1.0
104
10-2
100
102
104
MFP [nm]
~ =jmax as an accumulation Fig. 11. Normalized lattice thermal conductivity j function of the phonon MFP at 100 K, 300 K, 500 K and 800 K for (a) HexagonalGaN and (b) Haeckelite-GaN.
longest phonon MFP of heat-carrying are 0.34, 0.23 and 0.13 lm for Hexagonal-GaN as shown as in Fig. 11(a), and 1.79, 1.02 and 0. 59 lm for Haeckelite-GaN as shown in Fig. 11(b). Fig. 11 also shows that the cutoff MFP of heat-carrying is longer at lower temperature. This is because higher temperature increases the anharmonic interactions, and thus reduces the phonon MFP. By examining Fig. 11(a), it yields the parameter l0 of 33.75 nm, 18.80 nm and 10.48 nm for Hexagonal-GaN at 300, 500 and 800 K, respectively. The corresponding values for l0 of HaeckeliteGaN in Fig. 11(b) is 12.73, 7.09 and 3.70 nm lm from 300 to 800 K. It can be found that phonons with the MFP shorter than l0 contribute around 50% of the total thermal conductivity. The parameter l0 is helpful in studying size effect on ballistic or diffusive phonon transport. 4. Conclusions In summary, we systematically investigate the phonon properties of wurtzite GaN and the corresponding two-dimensional single layer structure from first-principles calculations and phonon BTE. It is found that the thermal conductivity along [0 0 1] direction is 13% higher than that along [1 0 0] direction at room temperature for defect-free wurtzite GaN. The thermal anisotropy is significantly reduced from 42% to 11% with increasing temperature from 100 to 800 K. We also found that the thermal conductivity of monolayer GaN with hexagonal structure is higher than that with
Y. Jiang et al. / Computational Materials Science 138 (2017) 419–425
haeckelite structure by a factor of four. By comparing the specific heat, phonon group velocity and the phonon MFP between the two kinds of single layer structures, the underlying physics is uncovered to explain the thermal conductivity difference. Our calculation shows that the complex unit cell in Haeckelite-GaN single layer structure leads to a lower phonon group velocity and shorter phonon lifetimes. We further extract the thermal conductivity accumulation function with phonon MFP. It is found that the dominant phonon mean free path in Hexagonal monolayer is higher than that in Haeckelite. It is also found that the TA mode dominants the thermal conductivity for the GaN-single layer, which is not the same with that in graphene. Acknowledgements The authors thank Prof. H. Bao at Shanghai Jiao Tong University for useful discussion. The authors acknowledge the financial support from the Natural Science Foundation of China (51406034), and National Basic Research Program of China (2011CB707605). References [1] W. Chikhaoui, J.M. Bluet, M.A. Poisson, N. Sarazin, C. Dua, C. Bru-Chevallier, Appl. Phys. Lett. 96 (2010) 072107. [2] S. Ganguly, A. Konar, Z.Y. Hu, H.L. Xing, D. Jena, Appl. Phys. Lett. 101 (2012) 253519. [3] X.H. Wang, L.W. Guo, H.Q. Jia, Z.G. Xing, Y. Wang, X.J. Pei, J.M. Zhou, H. Chen, Appl. Phys. Lett. 94 (2009) 111913. [4] X. Li, Z.S. Liu, D.G. Zhao, D.S. Jiang, P. Chen, J.J. Zhu, J. Yang, L.C. Le, W. Liu, X.G. He, X.J. Li, F. Liang, L.Q. Zhang, J.P. Liu, H. Yang, Y.T. Zhang, G.T. Du, J. Vac. Sci. Technol. B 34 (2016) 041211. [5] B. Auld, Acoustic Fields and Waves in Solids, John Wiley & Sons, 1973. [6] A. Jezowski, B.A. Danilchenko, M. Bockowski, I. Grzegory, S. Krukowski, T. Suski, T. Paszkiewicz, Solid State Commun. 128 (2003) 69–73. [7] G.A. Slack, L.J. Schowalter, D. Morelli, J.A. Freitas, J. Cryst. Growth 246 (2002) 287–298. [8] M. Mukherjee, N. Mazumder, S.K. Roy, K. Goswami, Semicond. Sci. Technol. 22 (2007) 1258. [9] J. Zou, X. Lange, C. Richardson, J. Appl. Phys. 100 (2006) 104309. [10] L. Lindsay, D.A. Broido, T.L. Reinecke, Phys. Rev. Lett. 109 (2012) 095901. [11] X.F. Wu, J. Lee, V. Varshney, J.L. Wohlwend, A.K. Roy, T.F. Luo, Sci. Rep-Uk 6 (2016) 22504. [12] W.J. Ju, Z.Y. Zhou, Z.Y. Wei, Aip. Adv. 6 (2016) 065328. [13] Z.Y. Wei, Y.F. Chen, C. Dames, Appl. Phys. Lett. 102 (2013) 011901. [14] K.S. Novoselov, D. Jiang, F. Schedin, T.J. Booth, V.V. Khotkevich, S.V. Morozov, A. K. Geim, Proc. Natl. Acad. Sci. U.S.A. 102 (2005) 10451–10453. [15] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Science 306 (2004) 666–669. [16] M.X. Liu, Y.C. Li, P.C. Chen, J.Y. Sun, D.L. Ma, Q.C. Li, T. Gao, Y.B. Gao, Z.H. Cheng, X.H. Qiu, Y. Fang, Y.F. Zhang, Z.F. Liu, Nano Lett. 14 (2014) 6342–6347.
425
[17] A. Splendiani, L. Sun, Y.B. Zhang, T.S. Li, J. Kim, C.Y. Chim, G. Galli, F. Wang, Nano Lett. 10 (2010) 1271–1275. [18] L.K. Li, Y.J. Yu, G.J. Ye, Q.Q. Ge, X.D. Ou, H. Wu, D.L. Feng, X.H. Chen, Y.B. Zhang, Nat. Nanotechnol. 9 (2014) 372–377. [19] Z.Y. Al Balushi, K. Wang, R.K. Ghosh, R.A. Vila, S.M. Eichfeld, J.D. Caldwell, X.Y. Qin, Y.C. Lin, P.A. DeSario, G. Stone, S. Subramanian, D.F. Paul, R.M. Wallace, S. Datta, J.M. Redwing, J.A. Robinson, Nat. Mater. 15 (2016) 1166–1171. [20] N.A. Koratkar, Nat. Mater. 15 (2016) 1153–1154. [21] J.V.N. Sarma, R. Chowdhury, R. Jayaganthan, J. Appl. Phys. 113 (2013) 243504. [22] D.C. Camacho-Mojica, F. Lopez-Urias, Sci. Rep-Uk 5 (2015) 17902. [23] D.A. Broido, M. Malorny, G. Birner, N. Mingo, D.A. Stewart, Appl. Phys. Lett. 91 (2007) 231922. [24] M. Omini, A. Sparavigna, Physica B 212 (1995) 101–112. [25] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [26] G. Kresse, J. Furthmuller, Phys. Rev. B 54 (1996) 11169–11186. [27] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [28] A. Jain, A.J.H. McGaughey, Comput. Mater. Sci 110 (2015) 115–120. [29] S. Baroni, S. de Gironcoli, A. Dal Corso, P. Giannozzi, Rev. Mod. Phys. 73 (2001) 515–562. [30] G. Deinzer, G. Birner, D. Strauch, Phys. Rev. B 67 (2003) 144304. [31] L. Lindsay, D.A. Broido, T.L. Reinecke, Phys. Rev. B 87 (2013) 165201. [32] W. Li, J. Carrete, N.A. Katcho, N. Mingo, Comput. Phys. Commun. 185 (2014) 1747–1758. [33] T.L. Feng, X.L. Ruan, Phys. Rev. B 93 (2016) 045202. [34] A. Ward, D.A. Broido, D.A. Stewart, G. Deinzer, Phys. Rev. B 80 (2009). [35] X.J. Tan, H.Z. Shao, T.Q. Hu, G.Q. Liu, J. Jiang, H.C. Jiang, PCCP 17 (2015) 22872– 22881. [36] L.L. Zhao, S. Xu, M.C. Wang, S.C. Lin, J. Phys. Chem. C 120 (2016) 27675–27681. [37] S.S. Batsanov, Inorg. Mater. 37 (2001) 871–885. [38] H. Schulz, K. Thiemann, Solid State Commun. 23 (1977) 815–819. [39] L.C. Xu, R.Z. Wang, X.D. Yang, H. Yan, J. Appl. Phys. 110 (2011) 043528. [40] R.K. Wu, R. Hu, X.B. Luo, J. Appl. Phys. 119 (2016) 145706. [41] C. Mion, J.F. Muth, E.A. Preble, D. Hanser, Appl. Phys. Lett. 89 (2006) 092123. [42] N. Bonini, J. Garg, N. Marzari, Nano Lett. 12 (2012) 2673–2678. [43] A. Togo, L. Chaput, I. Tanaka, Phys. Rev. B 91 (2015) 094306. [44] L. Lindsay, D.A. Broido, J. Phys. – Condens. Matter 20 (2008) 165209. [45] X.K. Gu, R.G. Yang, J. Appl. Phys. 117 (2015) 025102. [46] L. Lindsay, D.A. Broido, N. Mingo, Phys. Rev. B 82 (2010) 115427. [47] B. Peng, H. Zhang, H.Z. Shao, Y.C. Xu, X.C. Zhang, H.Y. Zhu, Sci. Rep-Uk 6 (2016) 20225. [48] G.Z. Qin, Q.B. Yan, Z.Z. Qin, S.Y. Yue, M. Hu, G. Su, PCCP 17 (2015) 4854–4858. [49] L. Lindsay, D.A. Broido, Phys. Rev. B 84 (2011) 155421. [50] H. Sahin, S. Cahangirov, M. Topsakal, E. Bekaroglu, E. Akturk, R.T. Senger, S. Ciraci, Phys. Rev. B 80 (2009) 155453. [51] F. Yang, C. Dames, Phys. Rev. B 87 (2013) 035437. [52] Z.Y. Wei, J.K. Yang, W.Y. Chen, K.D. Bi, D.Y. Li, Y.F. Chen, Appl. Phys. Lett. 104 (2014) 081903. [53] A.S. Henry, G. Chen, J. Comput. Theor. Nanos. 5 (2008) 141–152. [54] A.J. Minnich, J.A. Johnson, A.J. Schmidt, K. Esfarjani, M.S. Dresselhaus, K.A. Nelson, G. Chen, Phys. Rev. Lett. 107 (2011) 095901. [55] K.T. Regner, D.P. Sellan, Z.H. Su, C.H. Amon, A.J.H. McGaughey, J.A. Malen, Nat. Commun. 4 (2013) 1640. [56] L. Lindsay, D.A. Broido, N. Mingo, Phys. Rev. B 83 (2011) 235428. [57] G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, F. Mauri, Nano Lett. 14 (2014) 6109–6114.