Calculation of thermoelectric properties of Bi2Te3

This section describes the calculation method to obtain the Bi₂Te₃ thermoelectric properties used in STAR Methods (design of high efficiency graded legs using Bi2Te3).

The thermoelectric properties are calculated using the density functional theory (DFT) (Hohenberg and Kohn (1964); Kohn and Sham (1965)) combined with the Boltzmann transport equation. For the DFT calculations, we use the generalized gradient approximation (GGA) parameterized by PBE (Perdew, Burke, and Ernzerhof) (Perdew et al. (1996)), and the projector augmented wave (PAW) pseudopotential (Blöchl (1994)); both of them are implemented in the VASP code (Kresse and Furthmüller (1996); Kresse and Joubert (1999)). The experimental lattice parameters for Bi₂Te₃ are used, while the internal coordinates are fully relaxed. The electronic band structure is calculated using the spin-orbit interaction. The k-point mesh of 36×36×36 is used. The electronic transport properties are predicted using the DFT band structure coupled with the Boltzmann transport equation within a rigid band approximation and the constant relaxation time approximation; they are implemented in BoltzTraP code (Madsen and Singh (2006); Ryu et al. (2017)). Note that we use the experimental band gap of 0.18 eV. The phonon thermal conductivity is calculated using phono3py code (Togo et al. (2015); Ryu and Oh (2016)). The force constants are obtained from the 240-atom supercell with the two-atom displacements using VASP code with the single k-point Γ and then the third-order phonon Hamiltonian is constructured. The three phonon scattering rates are calculated using the Fermi's golden rule. We also include the effective boundary scattering of 10 nm in addition to the three-phonon scattering. Then the thermal conductivity is calculated by integrating the conductivity on the phonon q-point mesh of 11×11×11.

Here we derive the integral Equation 7. For simplicity, we denote the term with Thomson heat and Joule heat in Equation 5 by fT(x):

Then the Equation 5 is ddx(κdTdx)+fT=0. If the solution Tsol of Equations 5, 6, and 26 is known, we may put κ(x):=κ(Tsol(x)) and f(x):=fTsol(x) to find a linear differential equation

Since this equation is linear, we can find the solution T by decomposing it into a homogeneous solution T1 and a particular solution T2 where T=T1+T2; see Figure S1. The T1 and T2 are solutions of

To solve the Equation 24, we integrate it over x to yield κ(x)dT1dx(x)=C for some constant C. Dividing both sides by κ and integrating from 0 to x, we have T1(x)−T1(0)=C∫0x1κ(x)dx. Imposing the boundary conditions yields C=−KTh−TcA and

To solve the Equation 25, we integrate it from 0 to x to yield κ(x)dT2dx(x)−C=−∫0xf(s)ds=:−F(x) for some constant C. Dividing both sides by κ and integrating from 0 to x, we have T2x−T20=−∫0xFsκsds+C∫0x1κsds. Imposing the zero boundary conditions yields

where δT:=∫0LF(x)κ(x)dx is a scalar quantity.

Summing up, we can check that the solution T=T1+T2 of Eqautions 23 and 6, and its gradient can be written as Equation 7.