Computational Simulations of Thermoelectric Transport Properties

Article information

J. Korean Ceram. Soc.. 2016;53(3):273-281
Publication date (electronic) : 2016 May 31
doi :
Thermoelectric Conversion Research Center, Creative and Fundamental Research Division, Korea Electrotechnology Research Institute (KERI), Changwon 51543, Korea
*Department of Advanced Materials Engineering, Hanbat National University, Daejeon 34158, Korea
Corresponding author: Min-Wook Oh, E-mail:, Tel: +82-42-821-1247 Fax: +82-42-821-1592
Received 2016 April 23; Revised 2016 May 10; Accepted 2016 May 12.


This review examines computational simulations of thermoelectric properties, such as electrical conductivity, Seebeck coefficient, and thermal conductivity. With increasing computing power and the development of several efficient simulation codes for electronic structure and transport properties calculations, we can evaluate all the thermoelectric properties within the first-principles calculations with the relaxation time approximation. This review presents the basic principles of electrical and thermal transport equations and how they evaluate properties from the first-principles calculations. As a model case, this review presents results on Bi2Te3 and Si. Even though there is still an unsolved parameter such as the relaxation time, the effectiveness of the computational simulations on the transport properties will provide much help to experimental scientist researching novel thermoelectric materials.

1. Introduction

Thermoelectric materials that recover waste heat and transform it into electricity or vice versa have attracted much attention because of their potential to help solve the current global energy crisis. Even though, up to now, thermoelectric materials have mainly been used in niche markets such as small refrigerators, a cooling device for laser diodes and cooling seats in automobiles, there have been many efforts to widen the applications of thermoelectric materials by enhancing their properties. The thermoelectric performance of the materials is characterized by the so-called dimensionless figure of merit, ZT=S2σT/κ, where S is the Seebeck coefficient; σ is the electrical conductivity; κ is the thermal conductivity, and T is the absolute temperature. 1,2) Conventional thermoelectric materials show ZT~1 at room temperature.1,2) Recently, the enhancement of ZT with ZT > 1 has been reported by many researchers. Heremans et al. reported that the value of ZT was 1.5 at 773K in Thallium doped PbTe due to the local distortion of the density of states and the resulting enhancement of the Seebeck coefficient.3) Pei et al. reported a ZT value of 1.8 at about 850K in Na-doped PbTe1−xSex alloys due to convergence of many valleys in the valence bands.4) Other improvements of ZT in the PbTe alloys have been reported in many instances.5) Nevertheless, the developed materials have limited use and commercialization due to environmental, economic and technical limitations including toxic elements in the materials, high cost elements in the materials, and a lack value for the ZT, respectively.

To design and discovery novel thermoelectric materials with high thermoelectric performance, computational simulations of the thermoelectric properties could have an important role in discovering novel materials as well as in guiding experiments and understanding the results. The computational simulations of three properties, the Seebeck coefficient, the electrical conductivity, and the thermal conductivity, have been a long-standing interest to many researchers. The core properties such as the band structure and the density of states, in calculating the transport properties, can be obtained with well-established computational tools. On the other hand, the computational tools for estimating the transport properties were recently established and released to the public. This review covers the latest advancements in computational simulations of thermoelectric properties mainly focusing on the first-principles calculations with the aim to increase attention on thermoelectric technologies in other researchers in computational materials science and highlight potential routes for the discovery of new thermoelectric materials. As a model case, this review presents the results of the electrical and thermal properties of Bi2Te3 and Si in detail.

2. Electrical Transport Properties: Seebeck Coefficient and Electrical Conductivity

The computational simulations for electrical conductivity and the Seebeck coefficient can be done with the first-principles calculations and the Boltzmann transport equation. The code BoltzTrap is the most widely used simulation program for electrical properties.6) This code was first written to use the eigen energy from the WIEN2k code.7) Furthermore, other codes for the first-principles calculations such as VASP and Quantum Espresso can also be used to produce the eigen energy for which the correspondent interface program was also developed.810) The BoltzWann code is also available to evaluate electrical transport properties with a maximally-localized Wannier function basis set.11)

Charge transport occurs when an electric field, and/or thermal gradient is present, and this phenomenon can be described as shown below:12)

(1) J=efν=σE

where J is the charge flux; e is the electronic charge; f is the charge distribution; σ is the electrical conductivity; E is the electric field and v is the charge velocity. Namely, if we know the distribution of the charge with respect to time and space, the flux of the charge can be determined. The charge distribution as a function of time can be described as follows: 12)

(2) ft+drdt·rf+dpdt·pf=(ft)c

where r is the position of the electrons; p is the momentum and the subscript c means collision. Thus, the above equation indicates the change in charge distribution after collision. And then, the relaxation time approximation is assumed here.

(3) ft=f-f0τ
(4) f-f0=Ce-t/τ

From equation (2)~(4), we can obtain

(5) f=f0+e(-f0ɛ)τv·E

and when we replace (1) with equation (5), the electrical conductivity can be described as

(6) σ=e2(-f0ɛ)ν2τ.

This equation can be rewritten as the tensor form within electronic structure calculations given as

(7) σ(ɛ)=e2NΩdɛ(-f0ɛ)n,kτn,kνn,kνn,kδ(ɛ-ɛn,k)

where Ω is the unit cell volume; e is the carrier charge; ɛ is the band energy; N is the number of k-points used in the calculation; f0 is the Fermi-Dirac distribution function; τ is the relaxation time; ν is the group velocity of charges and δ is the delta function.6,13) The subscripts k and n mean the crystal momentum and the band index, and the velocity, ν, can be estimated from the band structure with the following relationship:

(8) νn,k=1ɛn,kk

where ℏ is the reduced Planck constant.

When a temperature difference is present, the total electric field is different from equation (1) due to the field originating from the Seebeck effect, and then, equation (1) is changed to

(9) J=σE-σST.

And the heat flux (q) produced from the temperature difference can be described as

(10) q=TJs=STJ-κT

where Js is the entropy flux. Both equations (9) and (10) are called the Onsager relationship.14) Using equations (9) and (10), we can obtain the Seebeck coefficient and electronic thermal conductivity (ke) obtained from the band structure calculation, and the results are shown below:6,13)

(11) S=ekBNΩσ-1dɛ(-f0ɛ)(ɛ-μkBT)n,kτn,kνn,kνn,kδ(ɛ-ɛn,k)
(12) κe=kB2TNΩdɛ(-f0ɛ)(ɛ-μkBT)2n,kτn,kνn,kνn,kδ(ɛ-ɛn,k)-TσS2

where μ is the chemical potential, and kB is the Boltzmann’s constant. From the equations above, we can identify that the electrical conductivity, Seebeck coefficient, and electronic thermal conductivity can be directly estimated from the band structure. Only one parameter we do not know yet is the value of the relaxation time. If we assume that the relaxation time varies less within the energy window of interest, namely, assuming that the relaxation time is constant, the Seebeck coefficient can be derived from the band structure without any fitting parameter because the constant relaxation time in the numerator and denominator disappears. Of course, the relaxation time is dependent on the band energy in real materials. Most codes for electrical transport properties such as BoltzTraP and BoltzWann evaluate the properties with a constant relaxation time approximation and with energies of bands obtained from very fine meshes in the reciprocal lattices.

2.1 Case: Bi2Te3

Bi2Te3 is a well-known thermoelectric material for room temperature and is used in commercially available thermoelectric modules. Therefore, there are many reports on the electronic structures of Bi2Te3, and the results of transport properties simulations have been reported since the early 2000’s.1529) The crystal structure of Bi2Te3 is shown in Fig. 1, and the space group of the crystal structure is R m.30) The quintuple sequences of Te(1)-Bi-Te(2)-Bi-Te(1) are stacked along the c-axis in the hexagonal cell, where Te(1)-Te(1) has weak van der Waals bonding, and the others have a mixture of covalent and ionic bonding. It is well known from experiments that the transport properties show anisotropic features. The electrical and thermal conductivity are larger along the basal plane than along the trigonal axis which originates from a layered-like crystal structure of Bi2Te3. The major defects in Bi2Te3 are antisite defects which are responsible for carrier concentration.31)

Fig. 1

Crystal structure of Bi2Te3 in hexagonal (left) and rhombohedral (right) cell.

The characteristic electronic structure of Bi2Te3 and its related compounds such Sb2Te3 and Bi2Se3 are listed in Table 1, and the band structure of Bi2Te3 is also shown in Fig. 2. The band structure was obtained with the WIEN2k code with the following parameters: RMTKMAX=10, RMT=2.5 a.u., 72 k-points 3 in irreducible Brillouin zone and PBE-GGA.7,32) It is clearly seen that the spin-orbital interaction (SOI) should be included in the calculation to produce the exact number of carrier pockets, the position of the band edges, and the band gap. Without the SOI, the conduction band edge and the valence band edge are located at Γ, whereas they move to the symmetry lines with the SOI. The discrepancy between reports may originate from the differences in the lattice parameter used in the calculations and in the simulation codes and parameters such as the exchange-correlation potential and the use of van der Waals interactions.

The Characteristic Features of the Electronic Structure of Bi2Te3, Sb2Te3 and Bi2Te3

Fig. 2

Band structure of Bi2Te3 without (red) and with (blue) spin-orbit interaction. The k-path for the band structure is also shown below.

The electrical transport properties of Bi2Te3 have also been obtained with the Boltzmann transport equations shown in Fig. 3. The BoltzTraP code was used to evaluate the transport properties. In evaluating the Seebeck coefficient, the relaxation time is assumed to be constant (10−14 s), and the rigid band assumption is used to calculate the properties with respect to the Fermi energy. The anisotropic property is well seen in the power factor (S2σ), which is due to large anisotropic electrical conductivity, whereas a relatively less anisotropic Seebeck coefficient was evaluated. Experimentally, (SxxxSzz)/(Sxxx) < 0.1 and σxxσzz=4~7 were observed.35,36) The Seebeck coefficient is dependent on the density of state effective mass, md=Nv2/3(m1m2m3)1/3, whereas the electrical conductivity is inversely proportional to the effective mass for the conductivity, mcond = 3(1/m1+1/m2+1/m3)−1. Thus, strong anisotropy in the electrical conductivity is partially due to the ellipsoidal energy surfaces of Bi2Te3. Another reason for the anisotropy is the relaxation time difference. The layered-like crystal structure causes an anisotropic scattering rate, where the relaxation time along the trigonal axis is smaller than that along the basal plane.

Fig. 3

The Seebeck coefficient and power factor of Bi2Te3 obtained from the band structure. The relaxation time=10−14 s is used to evaluate the electrical conductivity.

2.2 Beyond the Constant Relaxation Time Approximation

The constant relaxation time approximation can be used for the Seebeck coefficient, but evaluation of the electrical conductivity should be done with some value of the relaxation time, for which we can fit the obtained electrical conductivity with the experiment and then acquire the value of the relaxation time. However, we need to overcome this obstacle in materials design or materials discovery. As shown in Fig. 4, there are various scattering mechanisms including deformation potential acoustic phonon scattering, ionized impurity scattering and so on. Among these scatterings, polar optic phonon scattering and ionized impurity scattering can be partly evaluated with the obtained parameters from the first-principles calculations such as the band gap, effective mass, high frequency and static dielectric constants. 37) The strength of the electron-phonon scattering can be also estimated in part from the first-principles calculations. For example, metals such as Cu and Ag have electron-phonon scattering as the dominant scattering mechanism of the electrical conduction at room temperature and above. In this case, the electrical conductivity can be also described as

Fig. 4

Various scattering mechanisms for electrons.

(13) 1σ=αsl-ph(TθD)50θDTx5(ex-1)(1-e-x)dx.

Here, θD is the Debye temperature; the constant αelph is proportional to λtrωD/ωp2, where λtr is the electron-phonon coupling constant, and ωp and ωD are the plasma and Debye frequencies, respectively.3840) λtr may be considered as a parameter to figure out the strength of the electron-phonon scattering, especially deformation potential acoustic phonon scattering. The value of λtr can be reduced to

(14) τ=2πλtrkBTFth

where Fth(1-0.0038θD2τ2)-1 when the temperature is greater than 0.7θD.4042) With the calculations of the phonon density of states with the first-principles calculations to obtain θD, the value of λtr can be estimated by comparing the results of the Boltzmann transport equation and the experiment.40) The estimated values of λtr for Cu and Ag are shown in Table 2.

Lattice Parameters Used in the Calculations, Estimated Bulk Modulus (B), the Debye temperature (θD), and Electron-Phonon Coupling Constant (λtr) for Ag and Cu40)

3. Thermal Conductivity

In this section, we briefly summarize the method to calculate the phonon thermal conductivity of thermoelectric materials by using the first-principles calculations. The thermal conductivity is composed of the electronic thermal conductivity and phonon thermal conductivity. The electronic thermal conductivity is computed from the electronic transport tensor by equation 12 shown already. Phonon transport properties are calculated based on the Boltzmann Transport Equation (BTE) combined with the Relaxation Time Approximation (RTA).38,46) Phonon frequencies and phonon life times are computed based on the supercell approach calculating the inter-atomic force constant.4753) The inter-atomic force constants(IFC) are obtained from the density functional theory. Then, 2nd and 3rd order phonon Hamiltonian is constructed. Then, the phonon life times and phonon thermal conductivities are calculated using the many body perturbation theory. Recently, the phonon thermal conductivity codes PHONO3PY and ShengBTE have become public to researchers.5458) Herein, the phonon life times and thermal conductivities are computed with the PHONO3PY code.5456)

To calculate the phonon thermal conductivity, we have to calculate the phonon lifetime. When the phonon Hamiltonian is only composed of 2nd order terms, there is no phonon scattering, and the phonon thermal conductivity should be proportional to the number of phonon channels. However, due to the anharmonicity in atomic potential surfaces, the phonon Hamiltonians are expanded with 3rd or higher order terms written as H = E0 + K + H2 + H3 + ···, where EO, K, Hn are the ground state total energy, the atomic kinetic energy, and the n-body crystal potential terms.5456)

Note that HHarmonic = K+H2 is the harmonic Hamiltonian for the phonon and written as a sum of harmonic oscillators, HH=λωλ(12+a^λ+a^λ). And the anharnonic 3rd order Hamiltonian is composed of three phonon processes, H3=λλλ(a^λ+a^-λ+)(a^λ+a^-λ+). The detailed equations can be found in ref. 54. Here, λ represents the phonon mode of the reciprocal vector q and phonon band index n. For 5-atom Bi2Te3 primitive cell, there are N = 5 atoms with 3N-3 = 15–3 = 12 phonon modes. Here, we can obtain the 3rd order phonon Hamiltonians. Once the 3rd order Hamiltonian is set, we can calculate the 3rd order phonon Hamiltonian, and the phonon life time can be calculated in a perturbative way. 5456) Then, finally, we can compute the phonon thermal conductivity. In the following paragraphs, we will introduce the results of the phonon thermal conductivity calculations for Si and Bi2Te3. The results are based on the use of the VASP and PHONO3PY code.

For the Si FCC structure, there are two Si atoms in the primitive cell. However, to compute the 2nd and 3rd order IFCs, we need to generate many supercell configurations containing atomic displacements. For the 2nd order, only single atomic displacements are needed. However, for the 3rd order, two atomic displacements are needed. If we consider the (4 × 4 × 4) supercell of Si FCC, then there are N = 128 atoms. Then, there are 3N displacements for one Si. Because we need two atomic displacements, there are huge numbers of displacement configurations exceeding 106 (since N = 128, 3N × 3N = 9N2 = 147,456). Fortunately, due to the symmetry operation,5457) the number of configurations can be reduced. In the PHONO3PY code, there are only 416 configurations. In the case of the Bi2Te3 hexagonal unit cell with 15 atoms, the number of configurations is very large (5108) due to less symmetry in the lattice of Bi2Te3. For low symmetry materials, there will be too many configurations. Then, it will be a good idea to reduce the number of configurations to cutoff the large distance displacement configurations. For Si, the evaluation of 5th neighborhood or higher neighborhood interactions changes the phonon thermal conductivity only negligibly.

Before generating displacement supercell configurations, one needs to check whether the primitive or unit cell is fully relaxed. Because the IFC analysis needs the atomic forces, the force noise is very harmful in the calculation of the phonon modes and 3rd order calculations. Here, the atomic forces for Si and Bi2Te3 are smaller than 0.00002 eVÅ−1. Also note that the phonon thermal conductivity is very sensitive to the lattice parameters, especially for thermoelectric materials.59) Here, we fix the lattice parameter of Si to be 5.47 Å from the PBE calculations and those of Bi2Te3 to be a = 4.3835 Å and c=30.487 Å from the experimental values. Also note that the choice of exchange correlation energy will alter the result of the phonon thermal conductivity.

We calculate the IFCs using the Si (4 × 4 × 4) FCC supercell and Bi2Te3 (4 × 4 × 1) hexagonal supercell. They contain 128 and 240 atoms. And the corresponding IFC supercell configuration numbers are 416 and 5108. To reduce the computational cost, we use the Γ point calculation for the electron Brillouin zone integration. We choose the PREC = ACCURATE tag for the VASP calculation and self-consistent-field (SCF) loop criterion of EDIFF = 1 × 10−8 eV following the suggestion of Togo.57) As a tip, we suggest storing the electronic wave function for the perfect non-displaced supercell. Using it, the DFT force calculation time is reduced in half. For us, we have 30 nodes, and each node contains 12 cores. For Bi2Te3, because each configuration takes about 400 to 500 seconds, the total calculation time is about 9 days. If we use more k-point, the computational time can take longer than several weeks or a month.

The 3rd order phonon process then can be calculated using the 2nd and 3rd order IFCs. The PHONO3PY code gathers the atomic forces for various configurations. The next time consuming step is to calculate the phonon self-energy. The phonons are scattered by three phonon processes. In this work, we used the phonon reciprocal lattice q mesh of 11 × 11 × 11 for Si and Bi2Te3. Because there are 113 q points, there are 113 × 113 × 113 three phonon processes. Note that it is a very time and memory consuming step if the unit cell is large, and the q mesh is fine. For the Bi2Te3 case, in the irreducible Brillouin zone, there are 146 grid points. Each q grid point phonon self-energy calculation is done in one node. For us, the total time of the calculations was about 4 days. If the unit cell is large, the calculation can crash, or the calculation time will be very long.

Then, the phonon thermal conductivity is obtained by gathering all phonon self-energies. Fig. 5(a) and (b) represent the phonon thermal conductivities of Si and Bi2Te3 with various calculation methods. We consider three different cases: PBE with a lattice parameter a = 5.47 Å, LDA with a = 5.47 Å and LDA with a = 5.43 Å. We find that the room temperature (300K) phonon thermal conductivity is 122 W/m/K for the PBE with a = 5.47 Å. The differences among the phonon thermal conductivities are about 20 – 70% for temperatures ranging from 10 to 100 K. When the temperature is larger than 300 K, the difference becomes less than 10% for PBE and LDA. However, for Bi2Te3, the choice of exchange-correlation energy is very sensitive to the thermal conductivity results. When we compared it to experimental results, the PBE results seems to be more reliable.

Fig. 5

Phonon thermal conductivities of (a) Si and (b) Bi2Te3 are plotted against temperature. Be aware that (a) is a log-log plot, while (b) is not.

Next, we examine the boundary effects on the phonon thermal conductivities, especially for PBE calculations shown in Fig. 6. The relaxation time of the phonon is modified following Matthiessen’s rule written as 1τ=1τPhonon+1τBoundary, where 1τPhonon and 1τBoundary are the phonon transition rate by phonon-phonon scattering and phonon-boundary scattering. The phonon boundary scattering time is assumed to be the boundary size (L) over the phonon group velocity (vg) written as 1τBoundary=vgL. For Si, the effect of boundary scattering is very large. For 100 nm boundary scattering, the phonon thermal conductivity becomes about half, meaning that the nanostructured bulk process will be a very effective way to reduce the phonon thermal conductivity, as previously reported in ref. 6062.

Fig. 6

Phonon thermal conductivities with various boundary mean free paths (Ls) are plotted for (a) Si and (b) Bi2Te3. With increasing L, the thermal conductivities are reduced. For Si, the effective size for the phonon thermal conductivity reduction is less than 1μm. For Bi2Te3, it is relatively small (less than 100 nm). Be aware that (a) and (b) are log-log plots.

For Bi2Te3, the phonon thermal conductivity with the boundary mean free path is examined. The phonon thermal conductivity curves are fitted using the following equations κphonon = κ0Tα, and the prefactor κ0 and exponent α are summarized in Table 3. The phonon thermal conductivity of Bi2Te3 at 300 K is very small (1.4 W/m/K) and when there is no boundary scattering, the phonon thermal conductivity is proportional to the inverse temperature (T−1). When the boundary scattering length is shortened, the phonon thermal conductivity is decreased, and the temperature dependency is weakened. The most effective boundary mean free path size is smaller than 100 nm. Note that when L is between 100 nm and 1 μm, the reduction of thermal conductivity is only 10 %. However, when 10 nm boundary scattering is introduced, the phonon thermal conductivity is severely reduced to 50%. It is thought that, in Bi2Te3, the smaller sized nanostructured bulk is needed for a classical thermoelectric material such as Bi2Te3. When we considered the average electron mean free path of Bi2Te3 (known to be about 10 – 20nm), the optimal nanostructuring size is about 10 – 20 nm to minimize the reduction of electrical conductivity. A similar result is also found for PbTe.63,64)

Phonon Thermal Conductivity Tensors of Bi2Te3 with the Boundary Mean Free Path (L) are Fitted to κphonon = κ0Tα

4. Utilization of Computational Simulation in Discovering Novel Thermoelectric Materials

Based on the first-principles calculations, all thermoelectric properties can be evaluated with the assumption that the relaxation time is constant so that we can screen candidates to discover novel thermoelectric materials. Madsen who developed the BoltzTraP code reported the screening procedures for thermoelectric materials containing the Sb element.65) Gorai et al., reported on the computational exploration of a binary compound of the A1B1 chemical space for thermoelectric performance.66)

5. Conclusions

We have shown that thermoelectric properties such as the electrical conductivity, the Seebeck coefficient, and the thermal conductivity could be evaluated from the electronic band structures and the inter-atomic force constants both obtained from the first-principles calculations within the density functional theory. It was shown that a combinational use of the simulator of electronic structures with the electrical Boltzmann transport equation solver and phonon calculation packages can give quantitative and comparative results. However, there are also obstacles to be overcome. The effect of scattering from various sources in existing real materials could not be involved in most calculations so that a few assumptions are used such as the relaxation time is constant. Despite this rigorous assumption, the computational simulations can be used to screen materials to discover novel thermoelectric materials. With the improvement of computational power, a new method for treating scattering is being developed so that a time will come in the near future when we can design thermoelectric materials with less approximation.


This work was supported by the National Research Foundation of Korea (NRF) Grant funded by the Korean Government (MSIP) (NRF-2015R1A5A1036133). B. Ryu was supported by the Korea Electrotechnology Research Institute (KERI) Primary research program through the National Research Council of Science and Technology (NST) funded by the Ministry of Science, ICT and Future Planning (MSIP) [No. 16-12-N0101-21, Development of design tools of thermoelectric and energy materials].


1. Hsu KF, Loo S, Guo F, Chen W, Dyck JS, Uher C, Hogan T, Polychroniadis EK, Kanatzidis MG. Cubic AgPbmSbTe2+m: Bulk Thermoelectric Materials with High Figure of Merit. Science 303(5659):818–21. 2004;
2. Snyder GJ, Toberer ES. Complex Thermoelectric Materials. Nature Mater 7(2):105–14. 2008;
3. Heremans JP, Jovovic V, Toberer ES, Saramat A, Kurosaki K, Charoenphakdee A, Yamanaka S, Snyder G. Enhancement of Thermoelectric Efficiency in PbTe by Distortion of the Electronic Density of States. Science 321(5888):554–57. 2008;
4. Pei Y, Shi X, LaLonde A, Wang H, Chen L, Snyder GJ. Convergence of Electronic Bands for High Performance Bulk Thermoelectrics. Nature 473(7345):66–9. 2011;
5. Poudeu PFP, Angelo JD, Downey AD, Short JL, Hogan TP, Kanatzidis MG. High Thermoelectric Figure of Merit and Nanostructuring in Bulk p-type Na1− xPbmSby− Tem+2. Angew Chem Int Ed 45(23):3835–39. 2006;
6. Madsen GKH, Singh DJ. BoltzTraP. A Code for Calculating Band-Structure Dependent Quantities. Comput Phys Commun 175(1):67–71. 2006;
7. Blaha P, Schwarz K, Madsen GKH, Kvasnicka D, Luitz J. WIEN2k, An Augmented Plane Wave Plus Local Orbitals Program for Calculating Crystal Properties Vienna University of Technology. Austria: 2001.
8. Kresse G, Furthmüller J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys Rev B 54(16):11169–86. 1996;
9. Kresse G, Furthmüller J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput Mater Sci 6(1):15–50. 1996;
10. Giannozzi P, Baroni S, Bonini N, Calandra M, Car R, Cavazzoni C, Ceresoli D, Chiarotti GL, Cococcioni M, Dabo I, Corso AD, Gironcoli S, Fabris S, Fratesi G, Gebauer R, Gerstmann U, Gougoussis C, Kokalj A, Lazzeri M, Martin-Samos1 L, Marzari N, Mauri F, Mazzarello R, Paolini S, Pasquarello A, Paulatto L, Sbraccia C, Scandolo S, Sclauzero G, Seitsonen AP, Smogunov A, Umari1 P, Wentzcovitch RM. QUANTUM ESPRESSO: a Modular and Open-Source Software Project for Quantum Simulations of Materials. J Phys: Condens Matter 21(39):395502. 2009;
11. Pizzi G, Volja D, Kozinsky B, Fornari M, Marzari N. BoltzWann: A code for the Evaluation of Thermoelectric and Electronic Transport Properties with a Maximally-Localized Wannier Functions Basis. Comp Phys Comm 185(1):422–29. 2014;
12. Nolas GS, Sharp J, Goldsmid HJ. Thermoelectrics: Basic Principles and New Materials Developments Springer-Verlag. Heidelberg: 2001.
13. Oh MW, Wee DM, Park SD, Kim BS, Lee HW. Electronic Structure and Thermoelectric Transport Properties of AgTlTe: First-Principles Calculations. Phys Rev B 77(16):165119. 2008;
14. Jeffrey G, Ursell TS. Thermoelectric Efficiency and Compatibility. Phys Rev Lett 91(14):148301. 2003;
15. Youn SJ, Freeman AJ. First-Principles Electronic Structure and its Relation to Thermoelectric Properties of Bi2Te3. Phys Rev B 63(8):851121. 2001;
16. Larson P, Mahanti SD, Kanatzdis MG. Electronic Structure and Transport of Bi2Te3 and BaBiTe3. Phys Rev B 61(12):8162. 2000;
17. Scheidemantel TJ, Ambrosch-Draxl C, Thonhauser T, Badding JV, Sofo JO. Transport Coefficients from First-Principles Calculations. Phys Rev B 68(12):125210. 2003;
18. Huang B-L, Kaviany M. Ab Initio and Molecular Dynamics Predictions for Electron and Phonon Transport in Bismuth Telluride. Phys Rev B 77(12):125209. 2008;
19. Lee S, von Allmen P. Tight-Binding Modeling of Thermoelectric Properties of Bismuth Telluride. Appl Phys Lett 88(2):022107. 2006;
20. Mishra SK, Satpathy S, Jepsen O. Electronic Structure and Thermoelectric Properties of Bismuth Telluride and Bismuth Selenide. J Phys Cond Matter 9(2):461. 1997;
21. Larson P, Lambrecht WRL. Electronic Structure and Magnetism in Bi2Te3, Bi2Se3, and Sb2Te3 Doped with Transition Metals (Ti-Zn). Phys Rev B 78(19):195207. 2008;
22. Larson P. Effect of p1/2 Corrections in the Electronic Structure of Bi2Te3 Compounds. Phys Rev B 68(15):1551211. 2003;
23. Yu Yavorsky B, Hinsche NF, Mertig I, Zahn P. Electronic Structure and Transport Anisotropy of Bi2Te3 and Sb2Te3. Phys Rev B 84(16):165208. 2011;
24. Hinsche NF, Yu Yavorsky B, Mertig I, Zahn P. Influence of Strain on Anisotropic Thermoelectric Transport in Bi2Te3 and Sb2Te3. Phys Rev B 84(16):165214. 2011;
25. Kim M, Freeman AJ, Geller CB. Screened Exchange LDA Determination of the Ground and Excited State Properties of Thermoelectrics: Bi2Te3. Phys Rev B 72(3):035205. 2005;
26. Pecheur P, Toussaint G. Electronic Structure and Bonding in Bismuth Telluride. Phys Lett A 135(3):223–26. 1989;
27. Pecheur P, Toussaint G. Tight-binding Studies of Crystal Stability and Defects in Bi2Te3. J Phys Chem Solids 55(4):327–38. 1994;
28. Ryu B, Kim BS, Lee JE, Joo SJ, Min BK, Lee HW, Park SD, Oh MW. Prediction of the Band Structures of Bi2Te3-Related Binary and Sb/Se-Doped Ternary Thermoelectric Materials. J Kor Phys Soc 68(1):115–20. 2016;
29. Oh MW, Ryu B, Lee JE, Joo SJ, Kim BS, Park SD, Min BK, Lee HW. Electronic Structure and Seebeck Coefficients of Bi2Te3, Sb2Te3, and (Bi0.25Te0.75)2Te3: A First-Principles Calculation Study. J Nanoelec Optoelec 10(3):391–96. 2015;
30. Nakajima S. The Crystal Structure of Bi2Te3−xSex. J Phys Chem Solids 24(3):479. 1963;
31. Oh MW, Son JH, Kim BS, Park SD, Min BK, Lee HW. Antisite Defects in n-type Bi2(Te, Se)3: Experimental and Theoretical Studies. J Appl Phys 115(13):133706. 2014;
32. Perdew JP, Burke K, Ernzerhof M. Generalized Gradient Approximation Made Simple. Phys Rev Lett 77(18):3865. 1996;
33. Thomas GA, Rapkine DH, Van Dover RB, Mattheiss LF, Sunder WA, Schneemeyer LF, Waszczak JV. Large Electronic-Density Increase on Cooling a Layered Metal: Doped Bi2Te3. Phys Rev B 46(3):1553. 1992;
34. Thonhauser T, Scheidemantel TJ, Sofo JO, Badding JV, Mahan GD. Thermoelectric Properties of Sb2Te3 Under Pressure and Uniaxial Stress. Phys Rev B 68(8):085201. 2003;
35. Scherrer H, Scherrer S. Bismuth Telluride, Antimony Telluride, and Their Solid Solutions. p. 211–238. CRC Handbook of Thermoelectrics In : Rowe DM, ed. CRC Press. Boca Raton: 1995.
36. Goldsmid HJ. The Electrical Conductivity and Thermoelectric Power of Bismuth Telluride. Proc Phys Soc 71(4):633. 1958;
37. Kang Y, Jeon SH, Son YW, Lee YS, Ryu M, Lee S, Han S. Microscopic Origin of Universal Quasilinear Band Structures of Transparent Conducting Oxides. Phys Rev Lett 108:196404. 2012;
38. Ziman JM. Electrons and Phonons p. 288–333. Oxford University Press. Oxford: 1979.
39. Bid A, Bora A, Raychaudhuri AK. Temperature Dependence of the Resistance of Metallic Nanowires of Diameter ≥ 15 nm: Applicability of Bloch-Grüneisen Theorem. Phys Rev B 74(3):035426. 2006;
40. Kim JY, Oh MW, Lee S, Cho YC, Yoon JH, Lee GW, Cho CR, Park CH, Jeong SY. Abnormal Drop in Electrical Resistivity with Impurity Doping of Single-Crystal Ag. Sci Rep 4:5450. 2014;
41. Allen PB, Pickett WE, Krakauer H. Band-Theory Analysis of Anisotropic Transport in La2CuO4-Based Superconductors. Phys Rev B 36(7):3926–29. 1987;
42. Mehta RJ, Zhang Y, Zhu H, Parker DS, Belley M, Singh DJ, Ramprasad R, Borca-Tasciuc T, Ramanath G. Seebeck and Figure of Merit Enhancement in Nanostructured Antimony Telluride by Antisite Defect Suppression through Sulfur Doping. Nano Lett 12(9):4523–29. 2012;
43. Allen PB. Empirical Electron-Phonon λ Values from Resistivity of Cubic Metallic Elements. Phys Rev B 36(5):2920–23. 1987;
44. Allen PB, Beaulac TP, Khan FS, Butler WH, Pinski FJ, Swihart JC. DC Transport in Metals. Phys Rev B 34(6):4331–33. 1986;
45. Savrasov SY, Savrasov DY. Electron-Phonon Interactions and Related Physical Properties of Metals from Linear-Response Theory. Phys Rev B 54(23):16487–501. 1996;
46. Srivastava GP. Physics of Phonons p. 122–174. CRC Press. Boca Raton: 1990.
47. Giannozzi P, De Gironcoli S, Pavone P, Baroni S. Ab initio Calculation of Phonon Dispersions in Semiconductors. Phys Rev B 43(9):7231. 1991;
48. Deinzer G, Birner G, Strauch D. Ab initio Calculation of the Linewidth of Various Phonon Modes in Germanium and Silicon. Phys Rev B 67(14):144304. 2003;
49. Broido DA, Malorny M, Birner G, Mingo N, Stewart DA. Intrinsic Lattice Thermal Conductivity of Semiconductors from First Principles. Appl Phys Lett 91(23):231922. 2007;
50. Esfarjani K, Stokes HT. Method to Extract Anharmonic Force Constants from First Principles Calculations. Phys Rev B 77(14):144112. 2008;
51. Tang X, Dong J. Pressure Dependence of Harmonic and Anharmonic Lattice Dynamics in MgO: a First-Principles Calculation and Implications for Lattice Thermal Conductivity. Phys Earth Planet Inter 174(1):33. 2009;
52. Tang X, Dong J. Lattice Thermal Conductivity of MgO at Conditions of Earth’s Interior. Proc Natl Acad Sci USA 107(10):4539–43. 2010;
53. Chaput L, Togo A, Tanaka I, Hug G. Phonon-Phonon Interactions in Transition Metals. Phys Rev B 84(9):094302. 2011;
54. Togo A, Chaput L, Tanaka I. Distributions of Phonon Lifetimes in Brillouin Zones. Phys Rev B 91(9):094306. 2015;
55. Katre A, Togo A, Tanaka I, Madsen GKH. First-Principles Study of Thermal Conductivity Cross-over in Nanostructured Zinc-Chalcogenides. J Appl Phys 117(4):045102. 2015;
56. Togo A, Tanaka I. First Principles Phonon Calculations in Materials Science. Scr Mater 108:1–5. 2015;
57. Webpage of phono3py Accessed on 11/04/2016.
58. Li W, Carrete J, Katcho NA, Mingo N. ShengBTE: A Solver for the Boltzmann Transport Equation for Phonons. Comput Phys Commun 185(6):1747. 2014;
59. Hellman O, Borido DA. Phonon Thermal Transport in Bi2Te3 from First Principles. Phys Rev B 90(13):134309. 2014;
60. Esfarjani K, Chen G, Stokes HT. Heat Transport in Silicon from First-Principles Calculations. Phys Rev B 84(8):085204. 2011;
61. Tian Z, Esfarjani K, Shiomi J, Henry AS, Chen G. On the Importance of Optical Phonons to Thermal Conductivity in Nanostructures. Appl Phys Lett 99(5):053122. 2011;
62. Qiu B, Tian Z, Vallabhaneni A, Liao B, Mendoza JM, Restrepo OD, Ruan X, Chen G. First-Principles Simulation of Electron Mean-Free-Path Spectra and Thermoelectric Properties in Silicon. EPL(Europhysics Letters) 109(5):57006. 2015;
63. Tian Z, Garg J, Esfarjani K, Shiga T, Shiomi J, Chen G. Phonon Conduction in PbSe, PbTe, PbSe1−xTex from First-Principles Calculations. Phys Rev B 85(18):184303. 2012;
64. Seklton JM, Parker SC, Togo A, Tanaka I, Walsh A. Thermal Physics of the Lead Chalcogenides PbS, PbSe, and PbTe from First Principles. Phys Rev B 89(20):205203. 2014;
65. Madsen GKH. Automated Search for New Thermoelectric Materials: The case of LiZnSb. J Amer Chem Soc 128(37):12140–46. 2006;
66. Gorai P, Parilla P, Toberer ES, Stevanović Vladan. Computational Exploration of the Binary A1B1 Chemical Space for Thermoelectric Performance. Chem Mater 27(18):6213–21. 2015;

Article information Continued

Fig. 1

Crystal structure of Bi2Te3 in hexagonal (left) and rhombohedral (right) cell.

Fig. 2

Band structure of Bi2Te3 without (red) and with (blue) spin-orbit interaction. The k-path for the band structure is also shown below.

Fig. 3

The Seebeck coefficient and power factor of Bi2Te3 obtained from the band structure. The relaxation time=10−14 s is used to evaluate the electrical conductivity.

Fig. 4

Various scattering mechanisms for electrons.

Fig. 5

Phonon thermal conductivities of (a) Si and (b) Bi2Te3 are plotted against temperature. Be aware that (a) is a log-log plot, while (b) is not.

Fig. 6

Phonon thermal conductivities with various boundary mean free paths (Ls) are plotted for (a) Si and (b) Bi2Te3. With increasing L, the thermal conductivities are reduced. For Si, the effective size for the phonon thermal conductivity reduction is less than 1μm. For Bi2Te3, it is relatively small (less than 100 nm). Be aware that (a) and (b) are log-log plots.

Table 1

The Characteristic Features of the Electronic Structure of Bi2Te3, Sb2Te3 and Bi2Te3

Eg (eV) Position # of carrier pockets References

Bi2Te3 0.11 Z–F Z–F Indirect gap 6 6 15
0.37 0.13 Γ–Z a–Γ Indirect gap 2 6 16
0.29 0.11 a–Γ a–Γ Direct gap 6 6 17
0.39 0.09 Z–F a–Γ Indirect gap 6 6 29
0.27 0.11 Z-F Z-F Direct gap 6 6 28
0.15 (experimental) Indirect gap 6 6 33

Sb2Te3 0.28 Γ–Z a–Γ Indirect gap 2 6 34
0.28 0.11 Γ–Z Z–F Indirect gap 2 6 29
0.027 0.13 Γ–Z Z-F Indirect gap 2 6 28

Bi2Se3 0.185 0.260 Γ Z-F Indirect gap 1 6 28

Eg: Band gap; SOI: spin-orbit interaction; CBM: Conduction band minimum; VBM: Valence band maximum.

Table 2

Lattice Parameters Used in the Calculations, Estimated Bulk Modulus (B), the Debye temperature (θD), and Electron-Phonon Coupling Constant (λtr) for Ag and Cu40)

a0 [Å] B [MPa] ΦD [K] λtr
Ag 4.165 88.6 201.8 0.12140)

Cu 3.637 136.2 327.1 0.13540)

Table 3

Phonon Thermal Conductivity Tensors of Bi2Te3 with the Boundary Mean Free Path (L) are Fitted to κphonon = κ0Tα

Component Fitting parameter Boundary mean free path (L)

0.1 nm 1 m 10 100 m 1 m 10 μm 100 m 1 m
κxx(L) κ0 0.025 0.711 24.640 235.250 405.870 432.610 435.580 436.040
α 0.000 −0.229 −0.609 −0.908 −0.987 −0.996 −0.997 −0.997
R2 0.9867 0.9977 0.9998 1.0000 1.0000 1.0000 1.0000

κzz(L) κ0 0.014 0.438 15.642 132.230 215.420 227.920 229.630 229.470
κ0 0.000 −0.247 −0.638 −0.923 −0.993 −1.002 −1.003 −1.003
α 0.9884 0.9982 0.9990 1.0000 1.0000 1.0000 1.0000

κxx(L)/κzz(L) @ 300K 1.786 1.796 1.862 1.940 1.961 1.964 1.963 1.964

κxx(L)/ κxx(L=1 m) @ 300K 0.017 0.131 0.521 0.897 0.988 0.999 0.999 1.000