# Fileset

[59_E-M2018813.pdf](https://mdr.nims.go.jp/filesets/1a67ad7c-58f8-406e-b5a9-e266f381fd5c/download)

## Creator

[Yukari Katsura](https://orcid.org/0000-0002-8905-2995), Hidenori Takagi, Kaoru Kimura

## Rights

[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[Roles of Carrier Doping, Band Gap, and Electron Relaxation Time in the Boltzmann Transport Calculations of a Semiconductor’s Thermoelectric Properties](https://mdr.nims.go.jp/datasets/c4c015ad-8671-4290-a96a-e940933edf0b)

## Fulltext

Roles of Carrier Doping, Band Gap, and Electron Relaxation Time in the Boltzmann Transport Calculations of a Semiconductor’s Thermoelectric PropertiesRoles of Carrier Doping, Band Gap, and Electron Relaxation Timein the Boltzmann Transport Calculations of a Semiconductor’sThermoelectric PropertiesYukari Katsura1,2, Hidenori Takagi3,4 and Kaoru Kimura11Department of Advanced Materials Science, Graduate School of Frontier Sciences, The University of Tokyo, Kashiwa 277-8561, Japan2Center for Materials Research by Information Integration (CMI2), Research and Services Division of Materials Dataand Integrated System (MaDIS), National Institute for Materials Science (NIMS), Tsukuba 305-0047, Japan3Department of Physics, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan4Max Planck Institute for Solid State Research, Heisenberg Str. 1, 70569, Stuttgart, GermanyAlthough there is a growing demand for first-principles predictions of the thermoelectric properties of materials, the contribution of variouserrors in Boltzmann transport calculations is not negligible. We conducted a typical first-principles calculation and a Boltzmann transportanalysis on a typical semiconductor (Si) at various temperatures Twhile varying the band gap ¾g, electron relaxation time ¸el, and phonon thermalconductivity ¬ph to demonstrate how the calculated thermoelectric properties, which are functions of the carrier doping level, are affected bythese parameters. Bipolar conduction drastically decreased zT via a degradation of the Seebeck coefficient S and an increase in the effectiveLorenz factor Leff, indicating the importance of a wide enough ¾g (several multiples of kBT or higher) for high zT. Thus, the underestimation of ¾g,which frequently happens in first-principles calculations, could induce large errors in calculations for narrow-gap semiconductors. Thecalculation of the electron thermal conductivity without Peltier thermal conductivity was found to limit the zT of typical semiconductors to below1. A small value of ¬ph/¸el, where ¬ph/¸el is the degree to which a material is a phonon-glass electron-crystal, was necessary to achieve a high zT.Fitting the calculations with experimental thermoelectric properties showed that ¸el can vary by an order of magnitude from 10¹15 to 10¹14 s,depending on both T and the samples. This indicates that the use of a fixed relaxation time is inappropriate for thermoelectric materials.[doi:10.2320/matertrans.E-M2018813](Received November 22, 2017; Accepted January 15, 2018; Published April 27, 2018)Keywords: thermoelectric properties, Boltzmann transport calculation, first-principles calculation1. IntroductionDiscoveries of new environmentally-friendly thermo-electric materials for solid state power generation and Peltiercooling devices are highly anticipated.1) Although thermo-electric materials have been intensively sought after for overhalf a century, many potentially thermoelectric compoundsremain unstudied or understudied. Experiments tend tounderestimate the potential of compounds for use asthermoelectric materials, since the optimum doping leveland the optimum method for introducing phonon scatteringare rarely achieved during the first set of the experiments.The best-known and the most popular benchmark for theefficiency of a thermoelectric material is the dimensionlessthermoelectric figure of meritzT ¼ S2·T¬¼ S2·T¬el þ ¬ph; ð1Þwhere S is the Seebeck coefficient, · is the electricalconductivity, ¬el is the electron thermal conductivity, ¬ph isthe phonon thermal conductivity, and T is the temperature.Thermoelectric conversion efficiencies for power generationand Peltier cooling increase monotonically with zT and reachthe Carnot efficiency at zT ¼¨. While the criterion forapplications has long been zT > 1, a number of bulk materialswith zT > 2 have recently been reported. The numerator ofeq. (1) is referred to as the power factor P Ô S2·, whichshould be high, especially for Peltier cooling applications.Equation (1) indicates that a high S, high ·, and low ¬ areneeded to achieve a high zT. Of these parameters, S, ·, and¬el cannot be controlled independently because of their strongdependence on the carrier doping level n Ô ne ¹ nh, where neand nh are the electron and hole concentrations. In dopedsemiconductors, a decrease in the carrier density increasesS and decreases ¬el but also decreases ·. As a result, zT islow both in low-n and high-n limits, and there is an“optimum n” (usually between 1019 and 1021 cm¹3) where zTis maximized.1) Comparisons between different materialsusually rely on the maximum reported zT values. However,the maximum zT of a material can only be achieved afterdeveloping effective methods to (1) dope carriers up to theoptimum n and (2) to decrease ¬ph while avoiding a decreasein the electron relaxation time ¸el. Therefore, researchershave to find the correct impurity elements or good defectsto control n and scatter phonons without also scatteringelectrons. Some materials already have intrinsic phononscattering mechanisms in their pure state and thus exhibitfairly high zT values without additional processes beingrequired to scatter phonons. The search for such intrinsicallylow-¬ph materials is currently being conducted with the helpof sophisticated phonon calculations and machine learning.2)Data-driven searches for intrinsically high-zT materialswithin vast lists of materials are rapidly progressing.Recently, high-throughput first-principles calculations3­7)have been combined with Boltzmann transport calculationcodes such as those in the Boltzmann Transport Properties(BoltzTraP8)) program package. In the TE Design Labdatabase4,5) the electronic properties that are related to thethermoelectric properties of over 2,300 compounds havebeen calculated theoretically. The database was carefullydesigned to prevent misleading people into believing that asingle representative value for zT can be predicted for eachMaterials Transactions, Vol. 59, No. 7 (2018) pp. 1013 to 1021Special Issue on Thermoelectric Conversion Materials X©2018 The Thermoelectrics Society of Japanhttps://doi.org/10.2320/matertrans.E-M2018813compound. Instead of directly predicting the maximum zT,the database provides the values of the ¢-factor,5) which isclosely related to zT via zT = v¢/(u¢ + 1). Properties thatare intrinsic to the electronic structure are included in ¢,and the extrinsic variables related to chemical potential andscattering mechanisms are excluded (they are represented byv and u, respectively). Compounds with a larger ¢ value havehigher chances of exhibiting a high zT after optimization ofn and suppression of the electron scattering.In such calculations of thermoelectric materials it isrelatively easy to reproduce the experimental zT values,because we can set arbitrary values for n, ¸el, and ¬ph ifnecessary. Observing the trends is also possible, because wecan fix these unknown variables as constants, such as thecommon approximation ¸el = 10¹14 s. However, predictingzT is not easy, because all these unknown variables need to bedetermined accurately.Since S, ·, and ¬el are strongly dependent on the carrierconcentration, we need to determine the chemical potential ®that represents the value of n in the sample. We also need tomake sure that the calculated value for ¾g is correct, becausedetermining ¾g from first-principles calculations often yieldsincorrect values, even though many thermoelectric materialsare narrow-gap semiconductors.The solutions of the Boltzmann transport equations8,9)yield ·, S, and ¬ as a function of both ® and T by integratingover energy ¾ such that·ð¾; T Þ ¼ e23Dð¾Þv2ð¾Þ¸elð¾; T Þ; ð2Þ·ð®; T Þ ¼Z·ð¾; T Þ � @fð¾; T Þ@¾� �d¾; ð3ÞSð®; T Þ ¼ � 1eT·ð®; T ÞZð¾� ®Þ·ð¾; T Þ � @fð¾; T Þ@¾� �d¾;ð4Þ¬elð®; T Þ ¼1e2TZð¾� ®Þ2·ð¾; T Þ � @fð¾; T Þ@¾� �d¾� Sð®; T Þ2·ð®; T ÞT : ð5ÞHere, e is the charge of an electron, ·(¾,T ) is the transportdistribution function (spectral conductivity), D(¾) is theelectron density of states, v(¾) is the electron group velocity,and f (¾,T ) is the Fermi-Dirac distribution. The actualcalculation is conducted in the form of tensors, and thevalues of these tensors are obtained by taking the average ofthe diagonal components of the tensors.8) Although ¸el(¾,T )is included in ·(®,T ) and ·(¾,T ), if ¸el(¾,T ) does not dependon ¾, then ¸el(¾,T ) cancels out during the calculation thusyielding exact values for the S(®,T ) tensor. This cancellationof ¸el(¾,T ) does not occur when calculating the ·(®,T ) and¬el(®,T ) tensors. As a result, we can only obtain therepresentative values for · and ¬el, as ·/¸el and ¬el/¸el asfunctions of ® and T.Note that there are two terms for ¬el in eq. (5). The firstterm ¬0 represents the normal diffusion of heat. The secondterm corresponds to the self-Peltier effect: the electric currentinduced by the thermoelectric voltage carries the heatbackwards from the low-T to the high-T end. In this paper,we refer to this second term as the Peltier thermalconductivity ¬P Ô ¹S2·T. When calculating ¬el, ¬P is oftenneglected because it is negligibly small in ordinary metals.9)However, this term may not be negligible in thermoelectricmaterials, because they are designed to maximize P.By using these calculable parameters ·/¸el, S, and ¬el/¸el,the full expression for zT becomeszT ¼ S2·T¬el þ ¬ph¼ S2ð·=¸elÞTð¬el=¸elÞ þ ð¬ph=¸elÞ: ð6ÞIn this expression, we can see that an unknown term (¬ph/¸el)still remains. While there have been intensive studies on ¬ph,few studies have ever been conducted on ¸el. Sophisticatedphonon calculations can predict the ¬ph values of purecompounds;10) however, the ¬ph values in high-zT samplesare significantly lower than those of pure materials owingto the extrinsic phonon scattering centers introduced inten-tionally or unintentionally during the fabrication process.Experimental ¬ph values estimated using the Wiedemann-Franz law ¬el = L·T are also unreliable, because the lawis not applicable to semiconductors, whose transportdistribution function is not linear as a function of energynear ®.9) It is also known that L in semiconductors tends todeviate from the theoretical L by ³20%1) or more.11)Here, the electron scattering time ¸el is the interval betweenelectron scattering events. Because of a lack of knowledgeregarding the real value of ¸el, a constant relaxation timeapproximation with a fixed ¸el is used for all T to roughlyestimate · and ¬el. A widely used approximation is¸el = 1 © 10¹14 s.3­8) Compared with their respective purecompounds, high-zT thermoelectric materials are expected tohave a higher concentration of electron scattering centers asthey originate from the materials that were introduced tocontrol n and scatter phonons.Scattering theories predict that ¸el has a T-dependence thatdepends on the dominant scattering mechanisms. Whenacoustic phonon scattering is dominant, ¸el is proportionalto T¹1/2 in covalent compounds, T¹1 in metals, and T¹3/2in ionic compounds. When ionic impurity scattering isdominant, ¸el is proportional to T3/2. There are also otherscattering centers such as grain boundaries and crystaldefects. The overall ¸el can be calculated from the inversesummation of ¸el of every scattering mechanism (­¸el¹1)¹1.There are a few theoretical approaches to calculate the¸el values of Si12,13) and PbTe14) from first-principles. Theirresults suggest that ¸el lies between 10¹15 and 10¹13 s.Additionally, there have been attempts to estimate the ¸elof Mg2Si1¹xSnx,15,16) skutterudites,17) SnSe,18) Ga2Ru,19) andAgGaTe220) by fitting the Boltzmann transport equationswith experimental · and carrier concentration values. Further,the electron mobility RH· = e¸el/m*, where RH is theexperimental Hall coefficient and m* is the calculatedeffective mass, can be used to estimate ¸el in experimentalsamples, if the correct value for m* is used. Chen et al. reporta list of the ratio of both the experimental and calculatedmobilities of 31 compounds21) under the assumption that¸el = 1 © 10¹14 s; the ratio ranged from 0.0005 to 10,corresponding to ¸el values between 7 © 10¹18 s and 9 ©10¹14 s. Although this mobility calculation has other possiblesources of errors such as effective mass inaccuracy andanisotropy, this implies that ¸el can vary by orders ofmagnitude between different samples.Y. Katsura, H. Takagi and K. Kimura1014When we ignore this unknown term ¬ph/¸el, the expressionbecomes zelT: the theoretical upper limit for zT in the limit of¬ph ¼ 0 is given byzelT ð®; T Þ �S2·¬el¼ S2ð·=¸elÞT¬el=¸el¼ Sð®; T Þ2Leffð®; T Þ: ð7ÞThe effective Lorenz factorLeffð®; T Þ �¬elð®; T Þ·ð®; T ÞT ¼ ¬elð®; T Þ=¸elð·ð®; T Þ=¸elÞTð8Þis another parameter that can be calculated without knowing¸el.In this study, we attempted to demonstrate how theseunknown parameters and calculation errors affect theBoltzmann transport calculation of a material’s thermoelectricproperties as this calculation is widely used in high-throughputfirst-principles calculations. We selected Si (Si­Ge alloy) asthe test material for the demonstration. It should be notedhere that more precise analysis models exist for the theoreticaland experimental transport properties of Si­Ge alloys.11,22­24)2. Calculation MethodThe electronic structure of Si was calculated usingWIEN2k code,25) which employs the full-potential linearizedaugmented plane wave (FLAPW) method. The initial latticeparameters and atomic coordinates were obtained at roomtemperature. For the exchange-correlation potential, thegeneralized gradient approximation (GGA) by Perdew,Burke, and Ernzerhof 26) was applied. In comparison, theTran-Blaha modified Becke-Johnson (TB-mBJ)27) potentialwas applied to the calculation via the GGA to achieve moreprecise values for the band gaps.28) Spin-polarization andspin-orbit interactions were ignored. The number of k-pointswas set to 8000 in the first Brillouin zone in the initialcalculation, and the mesh was later interpolated by a 5©denser mesh in reciprocal space using code in the BoltzTraP8)program package. The muffin-tin radius RMT was set to2.10 a.u., while the maximum k was set such that RMTkmax =7. The core-valence cutoff energy was set to Ecut = ¹6.0Ry.The convergence criteria for energy and charge were0.0001Ry and 0.0001 e, respectively. The Boltzmann trans-port analysis was conducted using BoltzTraP8) version 1.5.3;the code was slightly modified29) to include ¬P during thecalculation of the transport tensors.3. Results and Discussion3.1 Electronic structuresFigure 1 shows the density of states DOS and energydispersion curves for Si calculated using the GGA26) andTB-mBJ27) potentials. When GGA was used the ¾g wasunderestimated (around 0.571 eV), which is nearly half theexperimental ¾g for Si (1.1701 eV).27) This underestimationis a well-known, common problem of first-principlescalculations based on density functional theory (DFT).28)The use of TB-mBJ resulted in a better value of ¾g of³1.156 eV. Although the raw calculation result provides theFermi energy ¾F at the top of the valence band, ¾F at n = 0(undoped) is shifted to the middle of the band gap duringtransport calculations by BoltzTraP.29) In reality, ¾F moves upwith electron doping and moves down with hole doping.Figure 1(b) shows the band structure of Si alongrepresentative reciprocal lattice vectors. The valence bandtop is composed of one heavy band and one light band, whichare degenerate at ¥ point, so that the number of hole pocketsis two in p-type Si. On the other hand, the conduction bandbottom is located along the ¥-X line, so there are six electronpockets in n-type Si.It can be seen that the band structures obtained usingTB-mBJ and GGA are essentially similar, except for thevalue of ¾g. The concavity of the bands or the effectivemasses m* were slightly greater with TB-mBJ than withGGA. We used the calculation results based on theexperimental lattice parameter30) a = 5.4310¡ at 300K sincethere were negligible differences in the calculated value of S.Changes in the lattice volume from ¹5% to 10% resulted indifferences in S within «0.5 µV/K at n µ 1020 cm¹3 in theGGA calculations. An optimization of the lattice parametersfrom GGA calculations resulted in an increased latticeparameter a = 5.4797¡ (+2.7% in volume); using this valueof a, ¾g increased up to 1.212 eV in the TB-mBJ calculation.Since the agreement with the experimental ¾g was poorerwith the optimized a, we employed the results with theexperimental a in the following calculations.3.2 Carrier doping levelIn this study, we plotted the thermoelectric propertiesagainst the carrier doping level n in units of inverse cubiccentimeters by dividing the amount of additional charge bythe unit cell volume.32) This corresponds to the carrier dopinglevel or net carrier concentration, i.e., the difference betweenthe electron and hole concentrations. This n is expected to beconstant within a sample, because n is related to the initialdopant concentration in the sample and is not affected by thethermal excitation of electron­hole pairs across the band gap.We plotted the materials’ thermoelectric properties againstn on a logarithmic scale. This representation is moreinformative than using a linear scale, which shrinks theinformation around the optimum n into a “spiky” regionaround n = 0.3.3 Seebeck coefficientFigure 2(a) shows the values of S of Si at different Tplotted against log n. Both p-type and n-type Si showed very0 0.51 1.52 2.533.5-5 0 5D(ε)/ states eV-1ε−εF / eV(a) TB-mBJGGA43210-1-2-3-4-5εFε−ε / eVF(b)Si5ΓW L X W KFig. 1 (a) Density of states DOS and (b) band structure of Si calculatedusing the FLAPW method using GGA (black) and TB-mBJ (red)exchange-correlation potentials.Roles of Carrier Doping, Band Gap, and Electron Relaxation Time in the Boltzmann Transport Calculations ... 1015similar curves, even though the band structures of the valenceband and the conduction band are different. At equal n,the n-type Si showed slightly greater «S« than p-type Si,especially in bipolar states. This is possibly due to thedifference in number of carrier pockets: more carriers existnear the band edge in n-type Si, so that «S« of electronsbecomes greater than «S« of holes.At 300K, the absolute value of «S« decreased monotoni-cally with increasing n from 1016 to 1022 cm¹3. Extremelyhigh values of «S« were predicted at low n. When we compare«S« at fixed n, the «S«­n curves shifted upwards with increasingT. Above 700K, peaks started to appear in the «S«­n curves.The peak position moved to higher n with increasing T.These behaviors in S are in agreement with bipolarconduction occurring; this conduction is caused by thethermal excitation of electron­hole pairs across the band gap.The drastic decrease in «S« is induced by the increase in thetotal carrier concentration and by the cancellation of S bycarriers of the opposite sign. This effect is significant atlow n, because thermal excitation adds equal concentrationsof electrons and holes, which dominate over the carrierconcentration induced by dopants.The degradation of «S« by bipolar conduction was moresignificant when ¾g was smaller. Figure 2(b) shows thecalculated S at 900K for artificially set ¾g values (from theexperimental value 1.17 eV down to 0 eV) using the‘scissors-cut’ operation implemented in BoltzTraP. Once ¾gwas set to 1.17 eV in GGA, it roughly matched the value ofS calculated using TB-mBJ. However, when ¾g was set tosmaller values, «S« was reduced on the left side of the plot,indicating that the decrease in «S« at low n is influenced bybipolar conduction. When ¾g was decreased down to 0.75 eV(10 kBT ), it started to influence S at n µ 1019 cm¹3, which isthe optimum range of n for thermoelectric materials. When ¾gwas decreased down to 0.25 eV (3.3 kBT ), the influence ofbipolar conduction on S was as high as n µ 1020 cm¹3. When¾g was set to 0, the peak of «S« was around 1021 cm¹3;however, the peak value of «S« was very small.These results tell us two things about ¾g. First, we shouldselect a parent compound for thermoelectric materials thathas a wide enough ¾g. This concept is known as the 10 kBT-rule.33) At least for Si, bipolar conduction lowers the «S« ofthe optimum range of n when the ¾g is less than 10 kBT.Second, when the calculated value of ¾g is small, thepredictions of the thermoelectric properties will containa large uncertainty. This highlights the difficulty of usingfirst-principles calculations to determine a material’s thermo-electric properties. Although many thermoelectric materialsare narrow-gap semiconductors, their band gaps are oftenunderestimated in DFT calculations.Figure 2(c) shows the agreement of the experimental S andcalculated S(n) at 300K. The figure shows the calculated Sfor T = 300K in comparison with the experimental S ofp-type and n-type Si0.8Ge0.2 alloys.22) Although there wasa high agreement between calculated and experimental «S«,there were slight differences in «S«, especially for n-typeSi0.8Ge0.2. One of the possible origins for this is thatelectronic structure of Si0.8Ge0.2 is different from the onefor pure Si, especially in the conduction band.3.4 ·/¸el, ¬el/¸el, and the effective Lorenz factorFigure 3(a) and (b) show the n-dependence of thecalculated value of ·/¸el of Si. Almost no difference in ·/¸el was observed between p-type and n-type Si. At 300­600K, we can see that ·/¸el is almost proportional to n.The quantity ·/¸el was independent of T, as expected fromthe Drude model ·/¸el = ntotale2/m*, where ntotal = ne +nh ’ n under monopolar conduction. Above 700K, the·/¸el­n curves showed a bend at low n. This is because ofthe increase in ntotal caused by bipolar conduction. It shouldbe noted that the thermal excitation of carriers does notchange n.Figure 3(c) and (d) show the n-dependence of ¬el/¸elcalculated using two different expressions for ¬el. At thetemperature where the bipolar conduction is negligible (300­600K), ¬el/¸el increased monotonically with n, just as ·/¸eldoes (see Fig. 3(a) and (b)). This similarity can be explainedvia the Wiedemann-Franz law (¬el = L0·T ).Large discrepancies in the calculated results were observedwhen we excluded or included the Peltier term ¬P = ¹S2·T.When ¬el/¸el was calculated according to ¬el/¸el = ¬0 + ¬P,the value of ¬el/¸el was lower than when ¬el/¸el was calculated0 200 400 600 800 10001016 1017 1018 1019 1020 1021 1022|S|/μVK-1n / cm-3(a) Si      TB-mBJp-typen-type1300120011001000900800300 K0 200 400 600 800 10001016 1017 1018 1019 1020 1021 1022|S|/μVK-1n / cm-3(b)SiGGA withshifted εg900 Kp-typen-type(TB-mBJ)1.17eV1.000.750.500.250.000 100 200 300 4000 1 2 3 4|S| / μVK-1n / 1020 cm-3Si300 K(c) TB-mBJ(p)GGA(p)TB-mBJ(n)GGA(n)experiment (p)(n)Fig. 2 The dependence of the absolute values of calculated Seebeck coefficient «S« on the carrier doping level n in p-type and n-type Si (a)at various T between 300K and 1300K with intervals of 100K, (b) at 900K using scissors-cut band energies with ¾g values between 0and 1.17 eV in GGA calculation, with TB-mBJ calculation (¾g ³ 1.16 eV) as a reference. (c) The comparison between calculated andexperimental values of «S« for Si0.8Ge0.2 alloys (reported by Vining et al.22)) plotted against Hall carrier densities for p-type and n-type Siat 300K. The S­n curves of Si were calculated by using TB-mBJ, and GGA with unshifted ¾g µ 0.571 eV. The dark solid circlescorrespond to the samples, which were used in the calculation of ¸el in the following sections.Y. Katsura, H. Takagi and K. Kimura1016according to ¬el/¸el = ¬0. This difference was the largestaround n µ 1020 cm¹3, where P = ¹¬P/T tended to reach itsmaximum. This shows that the approximation ¬el = ¬0 isunsuitable, especially for thermoelectric materials that aredesigned to achieve a maximum P. In the metallic limit(n > 1021 cm¹3) and in the bipolar limit (n < 1018 cm¹3 forT = 1300K), the value of ¬el/¸el calculated with ¬P wasalmost equal that calculated without ¬P. The bipolar effectincreases ¬el not only by increasing ntotal, but also byincreasing the energy per carrier particle by ¾g. Thermallyexcited carriers possess additional energy equal to ¾g fromthe high-T end to the low-T end of the sample. Therefore,zT is suppressed under bipolar conduction both because ofthe cancellation of S and because of the increase in ¬el.By using both the calculated ¬el/¸el and ·/¸el values,we calculated the effective Lorenz factor Leff = ¬el/(·T ) =(¬el/¸el)(·/¸el)¹1T¹1 and compared the results to those of thetheoretical Lorenz factor L0 in the free-electron modelL0 ¼³2kB23e2� 2:45� 10�8 V2 K�2: ð9ÞFigure 3(e) and (f ) show the values of Leff calculated usingtwo different expressions for ¬el. When we calculated ¬elusing ¬0, Leff varied greatly and differed significantly from L0,except for in the metallic region where n > 1021 cm¹3.However, when we included ¬P, Leff was almost constant overa much wider range of n, except in the region influenced bybipolar conduction. When bipolar conduction was takingplace, the increase in ¬el/¸el resulted in an exponentialincrease in Leff. The Leff in the semiconducting region wascalculated to be around 1.67 © 10¹8V2K¹2 for p-type Siand 1.86 © 10¹8 V2K¹2 for n-type Si, which were 32% and24% lower than L0. Such decreases in Leff have beentheoretically predicted for Si by Bera et al.11) They reportedthat the Leff of Si is dependent on calculation methods, andit may be further decreased by using exact solutions for theBoltzmann transport equations beyond the relaxation timeapproximation.11)3.5 zelT: the electronic upper limit for zTThe n-dependence of zelT is thought to be material-dependent and a function of both n and T; it is believed thatthis dependence can be determined purely from its electronicstructure. However, as presented in Fig. 4(a) and (b), then-dependence of zelT calculated using different expressionsfor ¬el exhibited very different behaviors. When we used¬el = ¬0 + ¬P, the zelT of Si could significantly exceed 1 atlow n. The simple expression zelT = S2/L0 also yieldedsimilar zelT curves to the ones obtained with ¬P, althoughsome discrepancies were observed in the low-n region wherebipolar conduction was dominant.However, when we used the approximation ¬el = ¬0, astrange upper limit of zelT appeared at zelT = 1, indicating thatthere was no chance of achieving zT > 1. In our calculationsthis limit was observed for almost any semiconductors withsteep band edges. This is consistent with past theoreticalcalculations that included ¬P9,34­37) or used the Wiedemann-Franz law18) to successfully achieve zelT > 1, while studiesthat only used ¬07,38­41) only found values of zelT below 1.The origin of this false “wall of zelT = 1” can be explainedas follows. In the low-n limit of non-degenerate semi-conductors, ® (the center of the Fermi-Dirac distribution) islocated nearly at the middle of the band gap, where no hole/electron states exist. In normal semiconductors the DOS atthe band edges are very steep. As a result, the hole/electronstates are concentrated near the band edge energy ¾0, which iseither the top of the valence band (¾V) for holes or the bottomof the conduction band (¾C) for electrons. At low T, only afew carriers can be thermally excited up to these states, andthe difference in ¾ of each carrier from ¾0 is negligibly small,compared to the large ¾0-®. Then we can treat ¾-® as aconstant (¾0-®) for all carriers, which can be taken out of theintegrals of eqs. (3)­(5).1013101410151016101710181019102010211016 1017 1018 1019 1020 1021 1022στel-1  / Sm-1s-1n / cm-3(a) p-type Si TB-mBJ1300 K1200 K1100 K1000 K900 K1013101410151016101710181019102010211016 1017 1018 1019 1020 1021 1022στel-1  / Sm-1s-1n / cm-3(b) n-type Si TB-mBJ800 K700 K600 K500 K400 K300 K109101010111012101310141015101610171016 1017 1018 1019 1020 1021 1022κ elτ el-1  / Wm-1K-1s-1n / cm-3(c) p-type Si TB-mBJκ0κ0+κP109101010111012101310141015101610171016 1017 1018 1019 1020 1021 1022κ elτ el-1  / Wm-1K-1s-1n / cm-3(d) n-type Si TB-mBJκ0κ0+κP10-810-710-61016 1017 1018 1019 1020 1021 1022Leff / V2 K-2n / cm-3(e)p-type SiTB-mBJκ0κ0+κPL0σT10-810-710-61016 1017 1018 1019 1020 1021 1022Leff / V2 K-2n / cm-3(f)n-type SiTB-mBJκ0κ0+κPL0σTFig. 3 Carrier density dependences of calculated values of (a) ·/¸el in p-type Si, (b) ·/¸el in n-type Si, (c) ¬el/¸el in p-type Si, (d) ¬el/¸el in n-typeSi, (e) Leff in p-type Si, and (f ) Leff in n-type Si between 300 and 1300Kin intervals of 100K. 0.11 101001016101710181019102010211022z elTn / cm-3(a)p-type SiTB-mBJκ0+κPL0σTκ0 0.11 101001016101710181019102010211022z elTn / cm-3(b)n-type SiTB-mBJκ0+κPL0σTκ01300600700800900100011001200300 K1300600700800900100011001200300 KFig. 4 Dependence of zelT on the carrier doping level of (a) p-type and(b) n-type Si calculated using three different expressions for ¬el betweentemperatures of 300 and 1300K in intervals of 100K.Roles of Carrier Doping, Band Gap, and Electron Relaxation Time in the Boltzmann Transport Calculations ... 1017Sð®; T Þ �������!j¾0�®j��¾ � 1eT·ð®; T Þ ð¾0 � ®Þ·ð¾0; T Þ! � 1eTð¾0 � ®Þ ð10ÞSð®; T Þ2·ð®; T ÞT �������!j¾0�®j��¾ 1e2Tð¾0 � ®Þ2·ð¾0; T Þ ð11Þ¬0ð®; T Þ �������!j¾0�®j��¾ 1e2Tð¾0 � ®Þ2·ð¾0; T Þ ð12Þ¬Pð®; T Þ �������!j¾0�®j��¾ � 1e2Tð¾0 � ®Þ2·ð¾0; T Þ ð13Þ) zeT ¼ Sð®; T Þ2·ð®; T ÞT¬0ð®; T Þ þ ¬Pð®; T Þ�������!j¾0�®j��¾ 1 ð¬0 onlyÞ1 ð¬0 þ ¬PÞ�ð14ÞWhen we include ¬P, the denominator ¬0 + ¬P becomes 0 andzelT diverges to infinity in the low-n, low-T limit. However,when we neglect ¬P, zelT does not exceed 1. This implies thatzT > 1 can only be achieved when a large power factor (S2·)cancels out a part of ¬0 via the self-Peltier effect. Mahan andSofo have proposed that when the transport distributionfunction ·ð¾Þð� @fð¾;T Þ@¾ Þ is a Dirac delta function and ® islocated away from ¾, then ¬P completely cancels out ¬0 toyield an extremely high value of zT.34) We consider thatordinary semiconductors in the low-n limit may havetransport distribution functions that are like delta-functionsbecause of their steep band edges, thus fulfilling thiscondition for high zT.3.6 Effects of ¬ph and ¸el on zTAlthough knowing zelT is useful for predicting the upperlimit for zT, we need to estimate the actual values of zT byincluding the contribution of the ignored term ¬ph/¸el. Atheoretical estimation of this parameter is difficult, since both¬ph and ¸el are expected to strongly depend on the materialand its microstructure since phonons and electrons scatter atimpurities, grain boundaries, and crystal defects.As shown in Fig. 5, the introduction of ¬ph/¸el results in alarge reduction of zT at low n. The strangely high zelT at lown only occurs when neglecting ¬ph, as it is the dominant termof ¬ in the low-n region. Therefore, the value of zelT shouldnot be used to estimate zT in the low-n region. However, athigh n the value of zelT is more related to zT, because then zTis less influenced by introducing ¬ph/¸el as ¬el is dominant.This unknown term ¬ph/¸el can be understood to represent thematerial’s degree of phonon-glass electron-crystal, a conceptthat was first proposed by Slack.42)3.7 Estimation of ¬ph and ¸el from the experimentalvaluesIt is difficult to calculate ¬ph and ¸el from first-principlescalculations. Even if the values of ¬ph and ¸el of pure crystalswere calculated, these values will be decreased in the actualsamples owing to extrinsic electron and phonon scatteringmechanisms. In this study, we therefore analyzed whether¬ph and ¸el should instead be obtained by comparing thecalculation results with the experimental values.By combining the calculation results with the reportedexperimental values for Sexp(T ), ·exp(T ), and ¬exp(T ) we firsttried to find the dataset of values of ¾F that best agrees withthe experimental data. This method is useful when the carrierdensity is not available, although it requires that thecalculation reproduces the actual electronic structure. TheBoltzmann transport calculation can provide thousands ofdatasets for different values of (¾F,T ). The parameter ncalc,which is the doped carrier concentration, is determined by ¾F.The other parameters Scalc, (·/¸el)calc, and (¬el/¸el)calc are allobtained as functions of (¾F,T ). We searched for the ¾F thatbest reproduced Sexp and nexp. Instead of using the value ofnexp from Hall measurements, the corresponding dataset forT was found by searching for a dataset with Scalc = Sexp,because in general nexp at higher T was not reported. Twoor more datasets were found to have Scalc = Sexp when theS­n curve peaked; we selected the dataset with the mostreasonable value for n.Figure 6(a) shows the T-dependence of ¾V-® for p-typeand ¾C-® for n-type semiconductors that was estimated bycomparing the value of S calculated using TB-mBJ to thevalues of Sexp of three samples of Si0.8Ge0.2 alloys that werereported by Vining et al.22) and Dismukes et al.41) Here, theDOS of Si0.8Ge0.2 and Si was assumed to be the same.Vining’s p-type (No. 75) and n-type (No. 93) samples wereobtained via ball-milling of the alloys followed by hot-pressing. A Dismukes’ n-type sample (No. 1834) wasobtained using the zone-levelling crystal growth techniqueso that the sample would have a higher crystallinity andlarger grains compared with the Vining’s samples. Since ¾V-®for the p-type semiconductor yielded almost 0meV and ¾C-®for the n-type yielded about 6meV, we can see that ® islocated almost at the band edge, as expected for semi-conductors near the degenerate regions. The n estimatedfrom S in Fig. 6(b) was on the order of 1020 cm¹3, which isconsistent with the experimental values. The cause of theslight increase with T thus appears to be either calculationerrors or the thermal excitation of localized carriers. 0.01 0.11 101001016 1017 1018 1019 1020 1021 1022zTn / cm-3(a)Sip-type300 KzelT10111012101310141015 0.01 0.11 101001016 1017 1018 1019 1020 1021 1022zTn / cm-3(b)Sin-type300 KzelT10111012101310141015 0.01 0.11 10 1001016 1017 1018 1019 1020 1021 1022zTn / cm-3(c)Sip-type900 KzelT101310141015 0.01 0.11 10 1001016 1017 1018 1019 1020 1021 1022zTn / cm-3(d)Sin-type900 KzelT101310141015Fig. 5 Dependence of zT on the carrier doping level for different valuesof ¬ph/¸el in [Wm¹1 K¹1 s¹1] for (a) p-type Si at 300K, (b) n-type Si at300K, (c) p-type Si at 900K, and (d) n-type Si at 900K. The data for¬ph/¸el = 0 are presented as zelT in Fig. 4.Y. Katsura, H. Takagi and K. Kimura1018The effective relaxation time ¸effel was obtained as¸effel ¼ ·expð·=¸elÞcalcð15Þby using experimental values of the electrical conductivity·exp. Although ¸effel differs from the value of ¸el of the idealcrystals,12,13) ¸effel and ¸idealel should satisfy Matthiesen’s law inthe Boltzmann transport equations,1¸effel¼ 1¸idealelþ 1¸extel; ð16Þwhere ¸extel is the extrinsic relaxation time, which isdetermined by impurity doping and the microstructure.Figure 6(c) shows the T-dependence of ¸effel estimated forthe three samples. All samples exhibited a value of ¸effel on theorder of 10¹15 s, while a monotonic decrease of ¸effel wasobserved with increasing T. We also estimated the value of ¸elof the other thermoelectric materials by using the values ofSexp, ·exp, and ¬exp at their best values of zT. Although thecalculations were not reliable enough, the values of ¸el wereestimated to be between 10¹15 and 10¹14 s; single crystals43)tended to exhibit longer ¸el values than polycrystallinesamples.44)This analysis of ¸effel enabled us to determine the dominantscattering mechanism in the samples. Figure 6(d) shows thatin the two n-doped samples ð¸effel Þ�1 and T had a linearrelationship to each other. Therefore, the dominant scatteringmechanism at high T in these samples was expected to bescattering caused by an acoustic distortion potential (inmetals), where ¸el is proportional to T¹1. The scatteringmechanism in the p-type sample appears to be moreindependent of T.In Fig. 6(e), the effective phonon thermal conductivity wascalculated according to¬effph ¼ ¬exp � ¸effel ð¬el=¸elÞcalc; ð17Þusing the experimental thermal conductivity ¬exp. Since thevalue of Leff of the compound was smaller than of L0, thevalue of ¬effph calculated by this method was higher than thatusing the free-electron model ¬FEMph ¼ ¬exp � L0·T . TheDismukes’ sample with its large grains exhibited a higher¬effph than the Vining’s sample with its small grains. Theseresults imply that ¬ph is strongly dependent on the density ofthe grain boundaries, which act as phonon scattering centers.4. Remaining Sources of ErrorsThe sources of errors in first-principles calculation ofthe thermoelectric properties of materials can be classifiedinto three categories: over-simplification of the crystalstructure, inaccuracy in the electronic structure calculation,and inaccuracies in the transport calculation.An over-simplification of the crystal structure can occur incalculations of almost any thermoelectric material since thehighest zT values are usually achieved in extremely complexstructures.1) For carrier concentration control and phononscattering, the crystals are solid solutions of various impurityelements that are randomly distributed throughout the crystal.Defects and disorders can act in a similar way, making itimpossible to define unit cells. A number of thermoelectricmaterials have incommensurate or misfit-layer structures,which are difficult to define by periodic unit cells. Thepresence of such structures would modify the electronicstructure near the band gap so that it is different from that ofthe parent compound.The second type of error is related to the electronic structurecalculations. Underestimation (and sometimes, overestima-tion) of ¾g is an unavoidable error in DFT. This error is criticalfor narrow-gap semiconductors. However, if the bandgap wasat least open, the scissors-cut operation of the band energiesin BoltzTraP8) could be used to improve the accuracy.21) Theexperimental values of the band gaps can be obtained fromdatabases such as citrination.45) The use of other high-costexchange-correlation potentials that reproduce ¾g may also behelpful, although the errors in m* need to be small enough.The introduction of various interactions, such as spin-polarizations for magnetic elements and spin-orbit interactionsfor heavy elements, is also necessary to reduce errors.Magnetic insulators exhibit a finite ¾g only in spin-polarizedcalculations, and Mott insulators may require an electronlocalization potential U to open the band gaps. Spin-orbitinteractions should be dominant in many heavy-elementthermoelectric materials; these split the bands into two similarbands separated by several hundred millielectronvolts. At thecrossing-points of those bands, band reconstruction occurs togenerate bands with smallm* that are separated by a small gap.In some compounds, the electronic structures are sensitive tothe lattice volume and to atomic positions. In such cases acareful selection of the calculation conditions is necessary.-202468 400  600  800 1000 1200μ-ε VB orε CB- μ / meVT / K(a)n (BM)n (ZL)p (BM)101910201021 400  600  800  1000 1200n /  cm-3T / K(b)02468 10 400  600  800 1000 1200τ el / 10-15sT / K(c)024680  400  800  1200τ el-1 / 1014s-1T / K(d)01234 400  600  800 1000 1200κ ph / Wm-1K-1T / K(e)Fig. 6 Temperature dependence of (a) the chemical potential, where ¾V is the top of the valence band and ¾C is the bottom of theconduction band; (b) the carrier doping level; (c) the effective relaxation time; (d) the inverse of the effective relaxation time, with linearfits of the data points; and (e) the phonon thermal conductivity in comparison with ¬total­L0·T (displayed as pale lines and symbols) thatwas estimated using experimental values of S, ·, and ¬ of 0.08% boron-doped Si0.8Ge0.2 alloys reported in the literature. The ball-milled,pressure-sintered samples prepared by Vining et al.22) is shown as BM, and the zone-levelled, high-crystallinity sample by Dismukes isshown as ZL.Roles of Carrier Doping, Band Gap, and Electron Relaxation Time in the Boltzmann Transport Calculations ... 1019The third type of error is related to the Boltzmann transportcalculations. This type of error can be usually be mitigated byincreasing the density of the k-point mesh. In BoltzTraP, theuse of a coarse k-point mesh results in oscillations and noisein the thermoelectric properties at high n. The introductionof anisotropy would improve the calculations of anisotropiccompounds. A more precise treatment of ¸el, which should bedependent on energy and bands, may also improve accuracyof the Boltzmann transport calculations.5. ConclusionWe investigated how the Boltzmann transport calculationof the thermoelectric properties of a typical semiconductor(Si) from a first-principles calculation is affected by unknownvariables such as n, ¾g, and ¬ph/¸el. We determined that toachieve a high zT the parent compound should have a wide¾g, low ¬ph, and long ¸el. We demonstrated that calculatedthermoelectric properties such as S, ·/¸el, ¬el/¸el, zelT, andzT can vary by orders of magnitude when n is varied. Theprediction of an extremely high value for S was not foundto be special, as it was achieved in a normal semiconductorwith a wide enough ¾g in the low-n limit. These propertieswere strongly affected by ¾g, even though the estimation of¾g is difficult in DFT. This error is critical for thermoelectricmaterials, because narrow-gap semiconductors are favoredfor ease of carrier doping. We also found that ¬P is requiredfor calculations of ¬el when we need to show that zT canexceed 1. The calculated Leff was nearly equal to L0 inthe metallic limit of n > 1021 cm¹3; however, it was ³20%lower than L0 when n < 1020 cm¹3 and increased by ordersof magnitude when bipolar conduction became dominant.The values of zT were largely dependent on ¬ph and ¸el,although both parameters are strongly influenced by typicalengineering processes such as impurity doping and micro-structural control.We attempted to estimate the T-dependence of ®, n, ¸el, and¬ph by fitting the calculation results with the T-dependencesof Sexp, ·exp, and ¬exp provided in experimental papers onthe thermoelectric properties of certain materials. We showedthat ¸el was on the order of several 10¹15 s and varied with Tas expected for the dominant electron scattering mechanism.It was also shown that the values of ¸el differed betweensamples fabricated in different ways.To predict semiconductor’s thermoelectric propertiescorrectly, it is necessary to reduce many errors in the valuesof parameters related to the crystal structures, electronicstructure calculations, and Boltzmann transport calculations.We also found that it is necessary to careful check that thecalculated results are in agreement with the experimentaldata for any calculations of the thermoelectric properties ofsemiconductors from first-principles.AcknowledgmentsThis work was supported by KAKENHI grants (grantnumbers 23760647 and 16K14379) from the Ministry ofEducation, Culture, Sports, Science and Technology and bythe “Materials research by Information Integration” Initiative(MI2I) project of the Support Program for Starting UpInnovation Hub from Japan Science and Technology Agency(JST). We thank Jim Bailey, PhD, from Edanz Group(www.edanzediting.com/ac) for editing a draft of thismanuscript.REFERENCES1) G.J. Snyder and E.S. Toberer: Nat. Mater. 7 (2008) 105.2) A. Seko, A. Togo, H. Hayashi, K. Tsuda, L. Chaput and I. Tanaka:Phys. Rev. Lett. 115 (2015) 205901.3) G.K.H. Madsen: J. Am. Chem. Soc. 128 (2006) 12140.4) P. Gorai, D. Gao, B. Ortiz, S. Miller, S.A. Barnett, T. Mason, Q. Lv, V.Stevanović and E.S. Toberer: Comput. Mater. Sci. 112 (2016) 368.5) J. Yan, P. Gorai, B. Ortiz, S. Miller, S.A. Barnett, T. Mason, V.Stevanović and E.S. Toberer: Energy Environ. Sci. 8 (2015) 983.6) Q. Hao, D. Xu, N. Lu and H. Zhao: Phys. Rev. B 93 (2016) 205206.7) M. Núñez-Valdez, Z. Allahyari and A.R. Oganov: arXiv:1610.07824v1.8) G.K.H. Madsen and D.J. Singh: Comput. Phys. Commun. 175 (2006)67.9) T. Takeuchi: Mater. Trans. 50 (2009) 2359.10) J.E. Turney, E.S. Landry, A.J.H. McGaughey and C.H. Amon: Phys.Rev. B 79 (2009) 064301.11) C. Bera, M. Soulier, C. Navone, G. Roux, J. Simon, S. Volz and N.Mingo: J. Appl. Phys. 108 (2010) 124306.12) B. Qiu, Z. Tian, A. Vallabhaneni, B. Liao, J.M. Mendoza, O.D.Restrepo, X. Ruan and G. Chen: Europhys. Lett. 109 (2015) 57006.13) M. Fiorentini and N. Bonini: Phys. Rev. B 94 (2016) 085204.14) S. Ahmad and S.D. Mahanti: Phys. Rev. B 81 (2010) 165203.15) X.J. Tan, W. Liu, H.J. Liu, J. Shi, X.F. Tang and C. Uher: Phys. Rev. B85 (2012) 205212.16) H. Wang, W. Chu and H. Jin: Comput. Mater. Sci. 60 (2012) 224.17) J. Yang, P. Qiu, R. Liu, L. Xi, S. Zheng, W. Zhang, L. Chen, D.J. Singhand J. Yang: Phys. Rev. B 84 (2011) 235205.18) S. Hao, F. Shi, V.P. Dravid, M.G. Kanatzidis and C. Wolverton: Chem.Mater. 28 (2016) 3218.19) Y. Takagiwa, K. Kitahara and K. Kimura: J. Appl. Phys. 113 (2013)023713.20) W. Wu, K. Wu, Z. Ma and R. Sa: Chem. Phys. Lett. 537 (2012) 62.21) W. Chen, J.-H. Pöhls, G. Hautier, D. Broberg, S. Bajaj, U. Aydemir,Z.M. Gibbs, H. Zhu, M. Asta, G.J. Snyder, B. Meredig, M.A. White, K.Persson and A. Jain: J. Mater. Chem. C 4 (2016) 4414, SupplementaryMaterial.22) C.B. Vining, W. Laskow, J.O. Hanson, R. Van der Beck and P.D.Gorsuch: J. Appl. Phys. 69 (1991) 4333­4340.23) A.J. Minnich, H. Lee, X.W. Wang, G. Joshi, M.S. Dresselhaus, Z.F.Ren, G. Chen and D. Vashaee: Phys. Rev. B 80 (2009) 155327.24) G.A. Slack and M.A. Hussain: J. Appl. Phys. 70 (1991) 2694.25) P. Blaha, K. Schwarz, G.K.H. Madsen, D. Kvasnicka and J. Luitz:WIEN2k. An augmented plane wave plus local orbitals program forcalculating crystal properties, Vienna University of Technology, Austria(2001).26) J.P. Perdew, K. Burke and M. Ernzerhof: Phys. Rev. Lett. 77 (1996)3865.27) F. Tran and P. Blaha: Phys. Rev. Lett. 102 (2009) 226401.28) D.J. Singh: Phys. Rev. B 82 (2010) 205102.29) After the line “thermal(1:3,1:3)=kappa(1:3,1:3)” in fermiintegrals.f inBoltzTraP,8) we added the following code: “DO i=1,3; DO j=1,3; DOialp=1,3; thermal(i,j)=thermal(i,j)+temp*seebeck(ialp,i)*nu(ialp,j);ENDDO; ENDDO; ENDDO”.30) G. Celotti, D. Nobili and P. Ostoja: J. Mater. Sci. 9 (1974) 821.31) W. Bludau, A. Onton and W. Heinke: J. Appl. Phys. 45 (1974) 1846.32) In BoltzTraP,7) the number of additional charge per unit cell is given as“N” in case.trace. The unit cell volume, which needs to be the reducedcell volume, is written in case.outputtrans.33) G.D. Mahan: Solid State Phys. 51 (1998) 81.34) T. Fang, S. Zheng, T. Zhou, L. Yan and P. Zhang: Phys. Chem. Chem.Phys. 19 (2017) 4411.35) Y. Takagiwa, K. Kitahara, Y. Matsubayashi and K. Kimura: J. Appl.Phys. 111 (2012) 123707.36) G.D. Mahan and J.O. Sofo: Proc. Natl. Acad. Sci. USA 93 (1996)Y. Katsura, H. Takagi and K. Kimura1020http://www.edanzediting.com/achttps://doi.org/10.1038/nmat2090https://doi.org/10.1103/PhysRevLett.115.205901https://doi.org/10.1021/ja062526ahttps://doi.org/10.1016/j.commatsci.2015.11.006https://doi.org/10.1039/C4EE03157Ahttps://doi.org/10.1103/PhysRevB.93.205206https://doi.org/10.1016/j.cpc.2006.03.007https://doi.org/10.1016/j.cpc.2006.03.007https://doi.org/10.2320/matertrans.M2009143https://doi.org/10.1103/PhysRevB.79.064301https://doi.org/10.1103/PhysRevB.79.064301https://doi.org/10.1063/1.3518579https://doi.org/10.1209/0295-5075/109/57006https://doi.org/10.1103/PhysRevB.94.085204https://doi.org/10.1103/PhysRevB.81.165203https://doi.org/10.1103/PhysRevB.85.205212https://doi.org/10.1103/PhysRevB.85.205212https://doi.org/10.1016/j.commatsci.2012.03.039https://doi.org/10.1103/PhysRevB.84.235205https://doi.org/10.1021/acs.chemmater.6b01164https://doi.org/10.1021/acs.chemmater.6b01164https://doi.org/10.1063/1.4775602https://doi.org/10.1063/1.4775602https://doi.org/10.1016/j.cplett.2012.04.025https://doi.org/10.1039/C5TC04339Ehttps://doi.org/10.1039/C5TC04339Ehttps://doi.org/10.1063/1.348408https://doi.org/10.1103/PhysRevB.80.155327https://doi.org/10.1063/1.349385https://doi.org/10.1103/PhysRevLett.77.3865https://doi.org/10.1103/PhysRevLett.77.3865https://doi.org/10.1103/PhysRevLett.102.226401https://doi.org/10.1103/PhysRevB.82.205102https://doi.org/10.1007/BF00761802https://doi.org/10.1063/1.1663501https://doi.org/10.1016/S0081-1947(08)60190-3https://doi.org/10.1039/C6CP07897Dhttps://doi.org/10.1039/C6CP07897Dhttps://doi.org/10.1063/1.4729772https://doi.org/10.1063/1.4729772https://doi.org/10.1073/pnas.93.15.74367436­7439.37) Z. Fan, H.Q. Wang and J.C. Zheng: J. Appl. Phys. 109 (2011) 073713.38) Y. Katsura and H. Takagi: J. Electron. Mater. 42 (2012) 1365.39) L.F. Wan and S.P. Beckman: Phys. Chem. Chem. Phys. 16 (2014)25337.40) B.-Z. Sun, Z. Ma, C. He and K. Wu: Phys. Chem. Chem. Phys. 17(2015) 29844.41) S. Kessair, O. Arbouche, K. Amara, Y. Benallou, Y. Azzaz, M.Zemouli, M. Bekki, M. Ameri and B.S. Bouazza: Indian J. Phys. 90(2016) 1403.42) G.A. Slack: CRC Handbook of Thermoelectrics 407 (1995).43) J.P. Dismukes, L. Ekstrom, E.F. Steigmeier, I. Kudman and D.S. Beers:J. Appl. Phys. 10 (1976) 2899.44) Y. Katsura: Mater. Trans. 57 (2016) 1035.45) J. O’Mara, B. Meredig and K. Michel: JOM 68 (2016) 2031. URL:https://citrination.com/.Roles of Carrier Doping, Band Gap, and Electron Relaxation Time in the Boltzmann Transport Calculations ... 1021https://doi.org/10.1073/pnas.93.15.7436https://doi.org/10.1063/1.3563097https://doi.org/10.1007/s11664-012-2226-zhttps://doi.org/10.1039/C4CP03328Khttps://doi.org/10.1039/C4CP03328Khttps://doi.org/10.1039/C5CP03700Jhttps://doi.org/10.1039/C5CP03700Jhttps://doi.org/10.1007/s12648-016-0876-zhttps://doi.org/10.1007/s12648-016-0876-zhttps://doi.org/10.2320/matertrans.MF201612https://doi.org/10.1007/s11837-016-1984-0https://citrination.com/