First principles search for n-type oxide, nitride, and sulfide thermoelectrics Kevin F. Garrity∗ Material Measurement Laboratory, National Institute of Standards and Technology, Gaithersburg MD, 20899 (Dated: June 24, 2016) Oxides have many potentially desirable characteristics for thermoelectric applications, including low cost and stability at high temperatures, but thus far there are few known high zT n-type oxide thermoelectrics. In this work, we use high-throughput first principles calculations to screen transitionmetaloxides,nitrides,andsulfidesforcandidatematerialswithhighpowerfactorsandlow thermal conductivity. We find a variety of promising materials, and we investigate these materials indetailinordertounderstandthemechanismsthatcausethemtohavehighpowerfactors. These 6 materialsallcombineahighdensityofstatesneartheFermilevelwithdispersivebands,reducingthe 1 trade-off between the Seebeck coefficient and the electrical conductivity, but they do so for several 0 different reasons. In addition, our calculations indicate that many of our candidate materials have 2 low thermal conductivity. n u I. INTRODUCTION oxide thermoelectrics has focused on a relatively small J number of oxides, mostly binaries and perovskites, leav- 2 ingopenthepossibilitythatbetteroxidethermoelectrics 2 The need for clean efficient power generation has led to a renewed interest in thermoelectric materials, which exist. ] In this work, we use high-throughput density func- can directly convert a temperature gradient into electri- i c cal power. Thermoelectrics can take advantage of a va- tional theory (DFT) calculations[12–23] to identify s promising n-type thermoelectric oxides and related ma- - riety of heat sources, including solar or waste heat, to l cleanly generate electricity[1–5]. Conversely, they could terials from the Inorganic Crystal Structure Database r t be used in cooling applications via the Peltier effect. (ICSD)[24]. The large amount of work required to syn- m thesize, optimize, and measure thermoelectrics experi- There has been an extensive effort over recent years to . mentallymakethistypeoftheoreticalscreeningofcandi- t discover and optimize materials with high zT, a dimen- a date materials particularly desirable. Similar techniques sionless thermoelectric figure of merit. While there has m have been used successfully to study the thermoelectric been significant progress in this area, existing materials - have not yet managed to provide a combination of high behavior of a variety of materials, including oxides[25– d 36]. While a fully first principles theoretical calculation n zT, low materials cost, and high durability that would of zT remains challenging, especially for oxides, which o result in widespread adoption. Much of the research on often have partially localized carriers, we can neverthe- c thermoelectricshasfocusedonhighmobilitysemiconduc- [ tors with small band gaps. Unfortunately, many of the less screen materials for both electronic and vibrational 2 most promising candidate materials have practical con- propertiesthatarenecessaryforgoodthermoelectricper- formance. Inthiswork,weperformsuchascreeningpro- v cerns(cost,toxicity,stability)whichhavethusfarlimited cedure,identifyingmanycandidatematerialswithcalcu- 2 their use in applications[6]. 2 latedthermoelectricpropertiesthataresimilartoorsur- In this work, we focus on the less explored group of 6 passexperimentallystudiedn-typeoxides. Furthermore, wide band gap transition metal oxides, as well as re- 1 we analyze the mechanisms behind the high thermoelec- lated nitrides and sulfides. While oxides are not usually 0 tricperformanceofthesematerials,findingthattheyfall . thought of as promising for thermoelectric applications, 1 into a small number of groups with similar properties. due to their typically low mobilities, the discovery of 0 goodthermoelectricperformanceinp-typeNa CoO and 6 x 2 1 otherlayeredCo-basedmaterialsresultedinanincreased II. METHODS : interest in this class of materials[1–3, 7]. n-type mate- v rials such as ZnO and SrTiO have also displayed high i 3 X powerfactors,buttheirzT valueshavethusfarbeenonly A. Calculating Thermoelectric Performance r moderate, due to high thermal conductivities[8–10]. De- a spitelimitedsuccessthusfar,oxidesprovidemanypoten- The dimensionless figure of merit for thermoelectrics tial advantages as thermoelectrics: 1) high thermal and can be written as chemical stability in air, 2) chemical versatility, allowing for extensive substitutions and doping, 3) low thermal zT =σS2T/κ, (1) conductivity,and4)lowcostmaterialsandprocessing[1– 3, 6, 11]. Thus far, much of the experimental work on where σ is the electrical conductivity, S is the Seebeck coefficient, κ is the total thermal conductivity (electrical pluslattice),andT isthetemperature. Thepowerfactor, which determines the electrical response of a material to ∗Electronicaddress: [email protected] a temperature gradient, is S2σ. 2 Unfortunately, the components of zT are not all easy calculating band structure related quantities. We find to calculate using first principles techniques. Within the that for most materials this correction has a relatively constant relaxation time approximation, which is used minor effect beyond increasing the band gap, and larger in this work, S can be calculated from a band struc- gapshavenodirecteffecton thermoelectric performance ture calculation without any adjustable parameters[37]. as long as the gap is already large enough to avoid Within the same approximation, it is possible to cal- significant thermal carrier excitation. We perform culate σ/τ , where τ is the electronic relaxation time. phonon calculations using DFT perturbation theory[51] e e Unfortunately, calculating τ from first principles re- without the +U correction. e mains challenging[38, 39]. This problem is especially Our main results are done on fully relaxed struc- severe for oxides, which often display complicated con- tures with initial coordinates from the ICSD. We use duction mechanisms and polaronic effects at low dop- PYMATGEN[52] to manipulate files from the ICSD to ing and low temperatures. In this work, we are con- setuptheinitialstructuresforrelaxation. Becausecalcu- cerned primarily with the opposite regime of high tem- lationsarerunatafixednumberofplanewaves,changes peratures and high doping, where the carrier mobilities in the unit cell during relaxation can effectively modify of oxides are typically larger[1–3]. Because we are com- the basis set. To ensure consistency between the ba- paring materials which are chemically similar, we expect sis set and the final structure, we run each relaxation them to have broadly similar electron scattering mech- three times, with a force convergence tolerance of 0.001 anisms. Therefore, we will use the quantity S2σ/τe to Ry/Bohr, an energy tolerance of 1 × 10−4 Ry, and a rank our candidate materials for suitability as thermo- stress tolerance of 0.5 Kbar. For phonon calculations, electrics. This estimate, which has been used in many we decrease the force tolerance to 5 × 10−5 Ry/Bohr. previousworks[26,34,40],shouldbesufficienttoatleast The BFGS algorithm as implemented in QUANTUM screen materials for those with band structures that are ESPRESSO was used for relaxations. promising for thermoelectric applications, even if deter- We use maximally localized Wannier functions as im- mining the final ranking of materials will require experi- plemented in WANNIER90[53–55] to interpolate band mental input. structures and BOLTZWANN, the WANNIER90 trans- Forreference,firstprinciplestechniquescanreproduce portmodule,tocalculatetheSeebeckcoefficientandcon- the thermoelectric properties of SrTiO with τ ≈ 4 fs 3 e ductivity under the relaxation time approximation[56, at room temperature[33], a typical value for oxides, but 57]. The use of Wannier interpolation allows us to per- some high mobility oxides like ZnO have much longer form accurate calculations of thermoelectric quantities scattering times[41]. All wide band materials have to be starting from relatively sparse first principles k-point doped in order to be used as thermoelectrics. In this grids, which we then interpolate to a k-point spacing of work, we use the rigid band filling model to estimate the 0.02 ˚A−1. This density is about ten times as dense as effects of doping, and we rank materials by S2σ/τ at e the first principles calculation along each direction in k- their optimum doping. space. The use of Wannier functions also allows us to Afteridentifyingmaterialswithpromisingbandstruc- calculate band structure derivatives analytically, which tures, we perform more computationally expensive accuratelytreatsdegeneratepointsintheBrillouinzone. phonon calculations for a limited number of candidate InordertouseWannierinterpolationforthiswork,we materials to estimate the lattice thermal conductiv- had to develop a procedure for automating the construc- ity, which is the dominant contribution to the thermal tion of localized Wannier functions. Because we are in- conductivity for most thermoelectrics, as described in terestedinthepropertiesofbothvalenceandconduction Sec. IIC. states, we normally include all possible orbitals which could contribute to states near the Fermi level (see sup- plemental materials for list). This is in contrast to many B. Band Structure Calculations applications of Wannier functions, which are concerned with either only the occupied bands or only a localized All of our calculations are based on DFT subspace of bands (e.g. d-orbitals). In these cases, the calculations[42, 43], as implemented in QUANTUM Wannier functions extend over several atoms and may ESPRESSO[44] and using the GBRV high-throughput be sensitive to the details of the localization procedure. ultrasoft pseudopotential library[45]. We use a plane For our application, we include all the relevant orbitals, wave cutoff of 40 Ryd for band structure calculations which results in Wannier functions that are atomic-like and 45-50 Ryd for phonon calculations. For Brillouin and strongly localized, even before the iterative localiza- zoneintegration, weuseaΓ-centeredgridwithadensity tion procedure, making the final result more robust. To of 1500 k-points per atom. calculate the Wannier functions, we use an inner ’frozen’ We use the PBEsol exchange-correlation window default of 4.5 eV above and below the conduc- functional[46], which provides more accurate lattice tion/valence band edges in order to ensure an accurate constants and phonon frequencies than other GGA interpolation of the band structure. In testing, our cal- functionals. We use the DFT+U technique[47–49], with culatedthermoelectricpropertiesareinsensitivetominor a U value of 3 eV for transition metal d-states[50], when variations in this window. 3 WhiletheWannierizationprocedureweoutlinedabove tothemodelinBjerget. al.,wheretheGruneisenparam- is relatively robust, there are a few situations that can eter (γ) and the Debye temperature (Θ )are calculated D resultinfailuresintheWannierization, whichareidenti- from the first principles phonon dispersion fied by monitoring the spread of the Wannier functions. First,iftherearesemicorestatesthatwereexcludedfrom (cid:115)5h¯2 (cid:82)∞ω2g(ω)dω Θ = n−1/3 0 (2) the Wannierization that overlap in energy with the va- D 3k2 (cid:82)∞g(ω)dω lence states, it will be necessary to include those states B 0 in the valence. Second, sometimes there are problems γ2 = (cid:80)i(cid:82) 8dπq3γi2qCiq (3) including orbitals with high energy (e.g. Sr d-states in (cid:80) (cid:82) dq C i 8π3 iq SrTiO ),asthesestatescanbecomedifficulttodisentan- 3 V ∂ω gle from the free electron-like bands when their energy γ = − iq (4) iq ω ∂V becomes too high. In both of these cases, we simply ad- iq just the orbitals that we include in the Wannierization where n is the number of atoms per unit cell, ω is the iq procedure by hand to fix these problems. Another is- angularfrequencyofphononmodeiatq-pointq,g(ω)is sue can arise if the ’frozen’ window overlaps with free the phonon density of states, γ is the mode Gruneisen iq electron-like bands. This can be fixed by adjusting this parameter, C is the mode specific heat, and V is the iq window downward to avoid overlap. We encountered all volume. The sum for the Gruneisen parameter is only of these problematic cases only rarely, and we adjust for performed over modes with ¯hω < k Θ . As per the iq B D them when necessary. discussion in Refs. 63, 69, we square γ to avoid cancel- iq One potential drawback of the Wannierization ap- lationbetweenpositiveandnegativeanharmonicitywhen proach is the necessity of including a large number of calculating γ. emptybandsinanon-self-consistentDFTcalculation,in Using the Debye temperature and Gruneisen parame- order to construct well-localized conduction band Wan- ter calculated above, we then insert them into the Slack nier functions. However, these extra bands are only re- model [63, 67, 68, 70, 71]: quiredonthethesparsek-pointgrid,andinpracticethe √ computationalcostofthisstepissmallerthantheinitial 0.849×334 κ (T)= structural relaxation. l 20π3(1−0.514γ−1+0.228γ−2) ×(cid:16)kBΘD(cid:17)2 kBMV 31 ΘD (5) ¯h ¯hγ2 T C. Thermal conductivity whereM istheaverageatomicmass. Wefindthatthe For typical thermoelectrics, the thermal conductivity DebyetemperatureandGruneisenparameterusedinthis is dominated by the lattice thermal conductivity (κl). waycontainalmostalloftheinformationofthefullBjerg First principles calculations of the thermal conductiv- model, as shown in table II, which presents correlations ity have been shown to be accurate for a wide variety of various models and quantities with our reference set ofmaterials[58–62]. Unfortunately,thesecalculationsre- of thermal conductivities. In fact, as shown in table II, quiretheanharmonicforceconstants,whicharetoocom- the quantity Θ /γ also has a high rank correlation with D putationally expensive to use as an initial screening tool thereferencethermalconductivities,althoughthereisno for high-throughput calculations, especially as many of computational advantage in using Θ /γ instead of Eq. D the materials we consider have large unit cells with rela- 5. tively low symmetry. The reason we do not use the full Bjerg model is that There have been various recent attempts to model the in some cases we found difficulty in fully converging the lattice thermal conductivity without performing a full acoustic modes for large unit cells. These modes depend calculationoftheanharmonicforceconstants[31,62–66]. on careful cancellation between all of the force constants Yan et. al. use a Debye-Callaway model with a con- to produce modes with zero eigenvalues at q =Γ, which stant Gruneisen parameter[31, 66]. Toher et. al.[65] is challenging to achieve numerically. This cancellation demonstrated that a modeled the Debye temperature can be enforced at Γ by using the acoustic sum rule to and the Gruneisen parameter, combined using the Slack modify the force constant matrix in various ways, but model[67, 68] for thermal conductivity, is useful as a we sometimes found results which depend on how the screening method for thermal conductivity. Another rule was enforced. Therefore, we opted to use a more screening method by Bjerg et. al.[63, 69] incorporates computationallyrobustprocedureappropriateforahigh- aspects of the first principles phonon band structure to throughputstudybyusingtheSlackmodelinsteadofthe approximate the lattice thermal conductivity. Bjerg model. In this work, we want a method which is both com- We find empirically in testing that using this combi- putationally feasible to apply to a few dozen compounds nation of the Slack model with the Bjerg definition of to use as a secondary screening procedure, and accurate theDebyetemperatureandGruneisenparameteroveres- enoughtoprovideareasonableorderingofcompoundsto timatesthethermalconductivity,sowereport70%ofthe consider for further study. We employ a method similar modelvalue,whichimprovesthequantitativeaccuracyin 4 principles thermal conductivities for a variety of half- Heusler compounds [62, 65, 72]. In addition, we com- pare with experimental thermal conductivities for a few ) additional oxides[73]. The results are shown graphically K 0 in Fig. 1 and correlations are given in table II, see the 0 3 supplementary materials for details. We find that our at chosen method is sufficient for screening materials for -1K those likely to have low thermal conductivity. As shown 1 in table II, the Spearman rank correlation[65] between -m W thereferenceandmodeledthermalconductivitiesforthe entire test set is 0.91, indicating we are able to identify ( ck promising materials. If we limit the dataset to materi- a sl als with κ < 50 W m−1K−1 at 300 K, a more realistic κ l rangeforcomplexoxides,therankcorrelationforthefull model drops to 0.83, which is still reasonable for select- ing materials to study further. We also note that when κ (W m-1 K-1 at 300 K) expt considering the entire dataset, the Debye temperature (see Eq. 2) alone has a rank correlation of 0.82 with the FIG.1: Comparisonofreferenceexperimentalandfirstprin- thermalconductivity,andwemakeitasaninitialscreen- ciplesthermalconductivities(x-axis)andtheSlackmodel(y- ing tool, as it is less computationally expensive than the axis). full model. However, directly calculating the Gruneisen parameter, rather than estimating it[31, 65, 66], signifi- cantly increases the accuracy of the model. Finally, ex- TABLE I: Correlations of various quantities with the ref- perimental thermal conductivities are sensitive to many erence thermal conductivities at 300 K (see Supplementary factors beyond the scope of this work, including defects Materials for list). κ is the model used in this work (see Slack and grain boundaries, which both makes comparisons Eq. 5) ; Bjerg refers to the full model of Ref. [63]. Pear- with experiments difficult but increases the possibility son and Spearman refer to the standard Pearson correlation of engineering materials to have lower thermal conduc- and the Spearman rank correlation, respectively. The first twocolumnsincludeallmaterials,thenexttwoarelimitedto tivities. materials with κ <50 W m−1K−1 at 300 K. l Quantity Pearson Spearman Pearson Spearman D. Effective Masses Low κ Low κ l l κ 0.83 0.91 0.65 0.83 Slack Bjerg[63] 0.93 0.92 0.69 0.88 In order to understand the conductivity and Seebeck Θ 0.74 0.82 0.43 0.71 coefficient, we consider several definitions of the effec- D 1/γ 0.39 0.75 0.66 0.60 tivemass. UsingourWannierinterpolation,wecalculate ΘD/γ 0.71 0.89 0.60 0.83 derivatives of the band structure analytically[56] at the conductionbandminimatofindtheeffectivemasstensor (m )−1 = 1 d2E . We will sometimes concentrate on ij h¯2dkidkj our testing but makes no difference in a ranking of com- mmin, the smallest eigenvalue of mij, which helps deter- pounds for those with the lowest thermal conductivity. mine the largest value of the conductivity tensor at low For materials with unstable phonon modes, we cannot temperature. We will also consider miso =(m1m2m3)31, calculateaGruneisenparameterinameaningfulwayus- the isotropic effective mass, where m1, m2, and m3 are ing purely harmonic calculations, as the unstable modes the eigenvalues of mij. must be stabilized anharmonically at finite temperature. A related band structure descriptor we calculate is a Therefore, we do not estimate a thermal conductivity versionoftheeffectivemassbasedonthedensityofstates for those compounds. In practice, we expect many ma- (DOS):[31, 74, 75] terials with unstable modes at zero temperature to have m (E) = ¯h2(cid:112)3 π4g(E)g(cid:48)(E) (6) lowthermalconductivityduetoanharmonicinteractions, DOS so the observation of unstable modes is already a useful (cid:82) dEg(E)m (E)(−df ) m (T,n ) = DOS dE . (7) indicator of anharmonicity. We present our calculated DOS d dEg(E)(−df ) thermal conductivities at 300 K, even though we expect dE thesematerialstobeusedathighertemperatures,where In this expression, g(E) is the DOS at energy E, f(E) the thermal conductivity will be lower. is the Fermi function, and T and n are the temper- d We establish the validity of this method for screen- ature and doping. This definition of m matches DOS ing the thermal conductivity by comparing the model m for a single parabolic band, but it is higher for eff with the experimental thermal conductivities of a vari- non-parabolic bands or when multiple bands contribute ety simple binary semiconductors, as well as the first to the conduction. These features allow m to give DOS 5 a good description of the Seebeck coefficient for many S2σ/τ materials[31, 74, 75]. E. Materials selection and screening procedure s) f m We are interested in discovering new n-type thermo- / S electric oxides, nitrides, or sulfides. As there are over 3 80,000 entries with oxygen in the ICSD, there is a need -0 1 to significantly limit our search space before proceeding. ( Previous experimental and theoretical work on thermo- τ / electricshassuggestedthatgoodthermoelectricstendto σ have anisotropic and non-parabolic bands and high den- sities of states, all of which can be created by empty d-orbitals[33, 40, 76–78]. Furthermore, materials with emptydorbitalscanusuallybedopedn-type,withsome |S| (μV/K) carriersoccurringnaturallyduetooxygenvacancies[1–3]. Therefore,inthiswork,wefocusonmaterialscontaining FIG. 2: |S|(µV/K) versus σ/τ (103 S ) for the entire at least one of Y, Sc, Ti, Zr, Hf, Nb, Ta, Mo, or W as e mfs dataset, at T=700K and doping at 1021 cm−3. The color well as at least one of O, N, or S. In order to limit com- scale indicates size of power factor, S2σ/τ . putational time, we restrict our search set to structures e with primitive unit cell volumes of less than 300 ˚A3. Our screening procedure proceeds in several steps, a) b) with each step becoming more computationally expen- sinivgem. Fatiresrti,alwse(7c6a6lcucolamtepothuendbsa)nadngdaepliomfianlaltoeutrhoeumrsettaarlts-. V/K) S / m fs) FSeoerbmecaktecroiaelffiswciietnht,aSg,apan(d59e2lemctartoenriicalcs)o,nwdeucctailvciutyla,tσe/tτhee. S (μ -3 τ (10 For materials with a high power factor (191 materials), σ / we then calculate the Debye temperature using Eq. 2. Ftoirnaalnlyd,/foorr alowsubDseebtyoef 1te9mmpaetreartiuarlse,wwitehchailgchulpaoteweκrSflaacck- c) m D O S d) 1/mmin l withEq. 5. ResultsarepresentedintableIIwithfurther details in the supplementary materials. τ σ/ τ 2S σ/ 2S III. RESULTS AND DISCUSSION We performed our screening procedure starting with 766 compounds from the ICSD as discussed in Sec. IIE, DOS at E (m )2/m F DOS min consistingof661oxides,60nitrides,and71sulfides(some compoundscontainmultipleanions). Ofthatlist,wefind FIG. 3: Thermoelectric properties of the entire dataset, at 592 materials with band gaps according to DFT+U (551 fixedT=700Knd=1021cm−3. a)mDOSversusS,b)1/mmin oxides, 53 nitrides, and 25 sulfides). For these materials, versus σ/τe, where mmin is the smallest component of mij, c) DOS at E versus the power factor, and d) m2 /m wecalculateS andσ/τe foravarietyoftemperaturesand F DOS min versus the power factor. dopings. If we sort this list of candidate materials by esti- mated power factor, S2σ/τ , at 700 K and with opti- e mizeddoping,werediscoverseveralcompoundsthathave perior to existing materials. In table II, we list some of previously been measured to have good n-type thermo- our most promising candidate materials, including those electric properties. For example, doped TiO2, SrTiO3, with high S2σ/τe and a few with moderate S2σ/τe and KTaO3,andTiS2haveallbeenmeasuredtohavepromis- low κl. We remove structures that are minor distortions ing power factors and show up highly in our list[1– of other structures on the list, that have missing atoms 3,9,10,33,35,79–82]. Thisgivesusconfidencethatour in the ICSD, or that are only theoretically proposaled screening procedure is useful. In addition to these pre- structures; full results are presented in the Supplemen- viously measured materials, there are a variety of com- tary Materials. pounds which have not been studied for thermoelectric We begin our analysis by looking for patterns in the applications and which may have properties that are su- entire dataset. First, we note that under the rigid band 6 TABLE II: Thermoelectric properties of the most promising compounds. The first five columns consist of the compound name, its space group number, the DFT+U band gap (eV), the isotropic effective mass. The next several columns are the DOS,m (seeEq. 6-7),S,andσ/τ ,allevaluatedatat700Kandfixed1021cm−3doping. ThenextcolumnisS2σ/τ atthe DOS e e optimizeddopingandat700K,andfollowingcolumnisthatoptimizeddoping. Thelastcolumnisthemodeledlatticethermal conductivity at 300 K. Materials with − as κ have unstable phonon modes or we were unable to converge the Gruneisen l parameter. Material Space Grp. Band Gap m DOS m S σ/τ S2σ/τ Opt. Doping κ iso DOS e e l (eV) (103 eV−1˚A−3) (µV/K) 10−3 S 103 W (1021cm−3) W mfs mK2fs mK CaTaAlO 15 4.0 2.8 4.08 5.0 -189 42.9 1.7 3 − 5 TiO 225 1.7 1.1 1.53 17.0 -305 9.6 1.7 8 − 2 LiNbO 161 3.5 2.1 3.92 7.0 -245 23.4 1.6 2 40 3 TiO 136 2.1 2.2 3.63 7.5 -250 20.7 1.6 3 <1 2 HfS 164 1.6 1.6 5.51 2.7 -165 56.9 1.5 1 3.4 2 NaNbO 63 1.8 3.8 1.89 0.8 -268 75.6 1.4 4 − 3 Ba TaInO 225 4.3 2.3 2.75 10.0 -285 12 1.4 5 − 2 6 YClO 129 5.1 1.1 3.20 6.6 -182 32.2 1.4 3 6.0 LiTaO 161 3.8 2.9 5.50 3.4 -204 36.3 1.4 1 34 3 Li ZrN 164 1.9 0.6 5.10 3.5 -178 42.5 1.4 2 19 2 2 CaTiSiO 15 3.2 4.0 4.28 6.4 -229 25.5 1.3 2 1.1 5 HgWO 15 2.4 1.5 4.22 3.5 -171 46.5 1.3 2 − 4 P WO 12 2.3 11.8 6.01 3.5 -150 56.1 1.3 0.7 6.5 2 8 ZrS 164 1.1 1.6 4.88 3.3 -172 45.1 1.3 2 22 2 TaPO 85 3.7 2.6 4.38 4.1 -160 49.5 1.3 2 − 5 LaTaO 36 3.4 2.9 5.13 3.7 -183 38.1 1.3 2 26 4 NbTl S 217 2.3 0.9 3.59 8.6 -246 19.2 1.3 2 − 3 4 SrTaNO 140 0.9 0.7 3.31 1.8 -109 122.6 1.3 0.3 − 2 TiS 164 0.4 0.4 4.73 3.7 -171 42 1.3 2 5.3 2 PbTiO 99 2.0 0.6 3.68 1.7 -489 35.3 1.2 7 16 3 Sr TaInO 225 4.3 2.1 2.65 10.0 -266 10.4 1.2 6 − 2 6 SrTiO 140 1.2 1.3 3.56 1.6 -165 39.5 1.2 6 − 3 NaNbN 166 1.1 0.8 3.70 3.3 -142 42 1.1 4 − 2 HfTaNO 25 2.0 0.8 4.22 2.4 -122 73 1.1 0.6 − 3 KNbO 99 1.6 3.0 1.82 1.3 -367 75.5 1.1 3 − 3 HfSiO 141 5.8 2.7 4.36 6.3 -276 18.8 1.1 2 190 4 Y O 164 4.2 1.0 3.44 3.7 -122 41.8 1.1 4 5.8 2 3 ZrO 225 3.7 1.1 3.78 4.6 -145 18.7 1.1 10 − 2 CdTiO 62 2.5 1.0 5.2 4.4 -205 20.7 0.9 2 <1 3 CaTiO 62 2.7 1.1 5.23 3.5 -170 26.8 0.8 1 <1 3 Y Ti O 227 3.1 1.1 2.21 10.7 -224 11.4 0.6 2 <1 2 2 7 model used in this work, most materials have optimal the following sections, we explore in more detail how dopings of about 1021 cm−3, which corresponds to dop- some of these individual materials achieve this higher ings on the order to 10%. While this is much higher than expected combination of S and σ/τ . e than typical semiconductor thermoelectrics, it is consis- Thetrade-offbetweenS andσ/τ makesfindingasim- e tent with the behavior of oxides like SrTiO3[1–3, 9, 10]. ple descriptor of the power factor in terms of features of Reaching such high doping values may be difficult in the band structure very difficult, even though we can practice, and will require further experimental and theo- relateσ andS individuallytofeaturesinthebandstruc- retical work (see for instance [83]). In this work we con- ture. IntableIII,wepresentSpearmanrankcorrelations centrate on identifying promising materials for further betweenseveralthermoelectricpropertiesandvariousde- optimization. scriptors ofthe band structure, and several ofthese rela- In Fig. 2, we plot the values of S versus σ/τ , at 700K tionships are plotted in Fig. 3. e and for a fixed doping of 1021 cm−3. The color scale For example, we find that the smallest value of the ef- shows the value of S2σ/τ [86]. There is a clear trade-off fectivemasstensor,m ,ishighlycorrelatedwithσ/τ , e min e between S and σ/τ , which is consistent with the be- and we plot this relationship in Fig. 3b. Unsurprisingly, e havior of simple parabolic bands where S ∝ m and materials with small effective masses usually have high eff σ ∝ 1/m , where m is the effective mass[4]. The conductivities, although this relationship can be compli- eff eff best materials do not maximize either S or σ/τ , but in- cated by anisotropy in the effective mass tensor or by e stead have S and σ/τ values in the center of observed many bands contributing to the conduction. Similarly, e range, but with a larger combination than is typical. In as shown in Fig. 3a, we find that we can model the See- 7 3 TABLE III: Spearman rank correlation matrix of various band structure descriptors and thermoelectric quantities at 2 700Kandatafixeddopingof1021cm−3. S2σ/τ isthepower e factor,|S|istheabsolutevalueoftheSeebeckcoefficient,σ/τ e ) 1 is the electrical conductivity, mmin is the minimum value of eV mDOijS,tihsetheffeedcetinvseitmyaosfsstteantesosra,tmtDheOSFeirsmdiefilenveedl.inCoErqr.el7a,tiaonnds gy ( 0 r e with absolute value above 0.65 are in bold. n E 1 Quantity S2σ/τ |S| σ/τ m m DOS e e min DOS S2σ/τ − 0.44 0.07 -0.03 0.22 0.67 e 2 |S| 0.44 − -0.81 0.66 0.84 0.43 σ/τ 0.07 -0.81 − -0.81 -0.83 -0.21 e 3 m -0.03 0.66 -0.81 − 0.74 0.28 R Γ X M Γ min m 0.22 0.84 -0.83 0.74 − 0.32 DOS FIG. 4: Band structure of cubic SrTiO . Energies are rela- DOS 0.67 0.43 -0.21 0.28 0.32 − 3 tive to the Fermi level. beck coefficient as S ∝m (T,n ) (see Eq. 7). three classes of materials, although some materials fall DOS d Despite these relatively strong relationships for S and into several classes. In general, the mechanisms for high σ in terms of certain definitions of the effective mass, powerfactorsconsistofcombiningalargenumberofflat we find that the combination of m2 /m has only a bands near the conduction band minimum with at least DOS min weakcorrelationwithS2σ/τ ,asplottedinFig. 3d. The somehighlydispersivebands. Thiscombinationallowsa e problem is that as shown in table III, the two definitions large number of carriers, some of which are in dispersive of the effective mass are strongly correlated with each bands, increasing σ, while at the same time keeping the other, and dividing one by another does not produce a Fermi level near the conduction band minimum, where usefuldescriptor. Otherquantitieslikem havesimilar |S| is largest, mitigating the typical trade-off. iso problems. Furthermore,alloftheeffectivemasses,aswell as S and σ/τ , are strongly correlated or anti-correlated e witheachother,butnonearebythemselvesstronglycor- A. Symmetry driven degeneracy related with the power factor. This can be understood in part by looking at Fig. 2, which shows that the best This group of promising thermoelectrics contains ma- materials do not lie at the extreme of either S or σ/τ , terials that have symmetries (or near symmetries) which e buttheyinsteadhaveanatypicalrelationshipbetweenS cause degeneracies in the conduction band minimum. and σ/τ . Finding a simple descriptor of that relation- This degeneracy increases the DOS for any given Fermi e ship is difficult when the stronger trend is the trade-off level, relative to a material without degeneracies. These between S and σ/τ . In addition, in many anisotropic degeneracies can be due to a single degenerate minimum e materials, the Seebeck coefficient and conductivity are inBrillouinzone, ortheycanbedueto abandstructure not maximized in the same direction, which makes find- with a conduction band minimum which is repeated due ing a simple descriptor for the maximum power factor to symmetry. more difficult. Finally, as we will see below, all of our In addition to having degeneracies, the conduction best materials have unusual band structures with some bands of these materials all consist of empty tran- combination of high anisotropy, non-parabolic behavior, sition metal d-orbitals that have highly anisotropic and multiple bands contributing to the conduction, and dispersions[33,35,40,78]. Theseanisotropicbandstruc- identifyingthesequalitiesrequiresgoingbeyondasimple tures allow the material to have both low m and eff effective mass description of parabolic bands. high m bands at the same minimum, combining large eff The best simple descriptor we found for S2σ/τ does Seebeck coefficients with the high conductivity. Simi- e not include any effective mass, but instead is just the lar degeneracies and anisotropic bands are behind the DOS evaluated at the relevant doping and temperature, high power factors of several semiconducting materi- as shown in Fig. 3c. While there is a clear relationship als which rely of empty p-orbitals instead of empty d- between the DOS and the power factor, many high DOS orbitals[27, 31, 40, 78]. materials have low power factors, making a high DOS a All of these features are present in the band struc- useful design criterion but not a sufficient condition for ture of cubic SrTiO , as shown in see Fig. 4, which is 3 good thermoelectric performance. knowntobeagoodn-typethermoelectric. SrTiO hasa 3 In the following sub-sections, we will investigate the single triply degenerate conduction band minimum at Γ band structures of some of the materials in table II, in due to the t states originating from the Ti-d orbitals. 2g ordertoevaluatethemechanismsthatallowthesepartic- These bands have highly anisotropic dispersions, with ular materials to minimize the trade-off between S and one nearly flat band (m = 6.3) and two highly dis- eff σ. We find that these materials separate roughly into persive bands (m = 0.4) going from Γ to X. The eff 8 a) a) b) c) Al b) c) Hf O S Ta Ca FIG. 5: a) Band structure of HfS . b-c) Side and top views FIG. 6: a) Band structure of CaTaAlO . b-c) Top and 2 5 of HfS . Larger gray atoms are Hf, smaller yellow atoms are sideviewsofCaTaAlO . LargeyellowatomsareCa,medium S 5 S. gray atoms are Ta, smaller magenta atoms are Al, smallest red atoms are O. combination of high degeneracy and high effective mass bands with very low effective mass bands, which allow by reducing its dimensionality and therefore increasing for high conductivity, is what allows SrTiO to escape its DOS is well-known, and has been shown in SrTiO 3 3 the normal trade-off between S and σ/τ . superlattices[76, 85]. We note that here we are consider- e Similar features are present in many of the other per- ing thermodynamically stable materials, rather than ar- ovskite variants which we find to be candidate thermo- tificial superlattices, nanowires, or quantum dots, which electrics (SrTiO , PbTiO , NaNbO , LiNbO , KNbO , should reduce manufacturing costs and increase thermo- 3 3 3 3 3 LiTaO ,Ba TaInO ,CaTiO ,Sr TaInO ,SrTaNO ). In dynamic stability. 3 2 6 3 2 6 2 addition, various phases of TiO and ZrO have similar Wepresenttwoexamplesofeffectivelylowdimensional 2 2 featureswhichleadtohighpowerfactors. Manyofthese materials which we predict have high power factors. materials have been studied as thermoelectrics before, First, in Fig. 5, we show the band structure and atomic and the mechanisms leading to their power factors are structure of HfS , which consists of weakly bound two- 2 relatively well-known[27, 31, 33, 35], so we will proceed dimensional hexagonal trilayers. The conduction bands with a discussion of the next two groups. areveryflatfromM totheminimumatL(m =4.5), eff characteristic of two-dimensional materials, but they are much more dispersive in other directions (m =0.3). eff B. Low dimensional conductors Second, in Fig. 6, we show the band structure and atomic structure of CaTaAlO , which consists of TaO 5 6 While all of the structures studied in this work are octahedra arranged into one-dimensional columns that three dimensional, in many cases the atoms which dom- are separated from each other by Ca ions and AlO 4 inate the conduction band minima are arranged in tetrahedra. This arrangement of Ta atoms leads to an two-dimensional layers, one-dimensional lines, or zero- anisotropic band structure with very flat bands from Γ dimensional dots, which leads to effectively low dimen- to Y (m = 2.4) but stronger dispersion from Γ to A eff sional conduction. In some cases, the material itself con- (m = 0.5). There are additional nearly degenerate eff sists of weakly bound layers, while in others there are conduction band minima at Y which also contribute to strong bonds in all three directions, but the transition the DOS. In HfS and CaTaAlO , the high DOS and the 5 metals are arranged in a low-dimensional way. strong anisotropy, which are caused by the low dimen- Reducing the effective dimensionality of a material sionality, create the conditions for a high power factor. results in highly anisotropic conduction bands and an Within our set of candidate thermoelectrics, ZrS , 2 increased DOS at the bottom of bands, which can in- TiS , HfS , YClO, CaTiSiO , WP O , TaPO , NaNbN 2 2 5 2 8 5 2 crease the power factor[77, 78, 84]. The idea of im- have quasi-two-dimensional structures, CaTaAlO , 5 proving the power factor of a candidate thermoelectric HgWO , LaTaO , and HfSiO have quasi-one- 4 4 4 9 3 metald-bands. Inaddition,inbothYClOandY O ,the 2 3 emptysandd-statesoftheYatomsarebothlocatednear 2 the conduction band minimum, which results in similar 1 behavior to the case where the orbitals come from dif- ) V 0 ferent atoms. Depending on the crystal structure and e gy ( 1 tfrhoemanCioun,sZ,nit, mAga,yCbde,pAouss,iborleHtogeonrgpin-beaenrdesmfprotymsI-bna,nSdns, r e n 2 Tl,Pb,orBitobecomedegeneratewithtransitionmetal E bands in this fashion. This type of engineering could al- 3 low one material to take advantage of the high Seebeck 4 coefficientsoftransitionmetaloxideswhileincorporating the higher mobility of semiconductors, which often have 5 L Γ X K Γ empty s or p orbitals from main group elements. The exact alignments of empty states from different atoms FIG. 7: Band structure of Sr TaInO . Bands with greater 2 6 is difficult to predict using DFT+U, so further study of then 35% In content are colored red, others are black. these materials to determine the band alignments more precisely may be necessary. dimensional structures, and NbTl S , Ba TaInO , 3 4 2 6 and Sr TaInO have quasi-zero-dimensional structures, 2 6 as their transition metals are separated from each other. D. Thermal Conductivity There are other possible advantages in using low- dimensional materials as thermoelectrics besides the in- Due to the high computational cost, we were not able creased DOS, including potentially lower thermal con- to calculate the thermal conductivity of our full dataset. ductivity, due to phonon scattering from the atomic For 191 compounds, we calculated the Debye tempera- layers, as well as the ability to physically separate ture (see Eq. 2), which is fairly strongly correlated with dopantsfromconductingchannels,whichcanreduceelec- thermal conductivity (see table II and supplementary tron scattering. One disadvantage is that the thermo- materials). Calculations of the Debye temperature are electric properties of low-dimensional materials will be bothlesssensitivetotheq-pointsamplingofthephonon anisotropic, resulting in reduced efficiency in polycrys- bandstructurethantheGruneisenparameterandrequire talline samples. the phonons at only one volume, making the computa- tions much faster. We find that in our set of transition metal oxides, ni- C. Accidental Degeneracies trides, and sulfides, there is relatively little variation in the Debye temperature (mean of 342 K, standard devia- Onefinalmechanismforincreasingthepowerfactorof tion of 66 K), as compared to our test dataset of simple an n-type oxide is to find or engineer a material with ac- binary and ternary semiconductors (mean 319 K, stan- cidental degeneracies of the conduction band minimum. dard deviation 234 K). This is likely due to the fact that Whilethiscanhappenforphysicallysimilarbandswhich all of our compounds contain ionic bonds between light happen to be degenerate at different points in the Bril- anionsandmediumtoheavytransitionmetals,whilethe louin zone, here we consider cases where the bands come test dataset contains a range of bonds, from covalent to fromdifferentorbitalsandhavedifferenteffectivemasses. ionic, and a range of atom masses. In both datasets, For example, in the double perovskite Sr TaInO , the there is a significant correlation between V−1 and the 2 6 conduction band consists of both Ta d-states and In p- Debye temperature, with a correlation coefficient of 0.68 states,whichhappentobeatsimilarenergies(seeFig.7, in the oxides, and 0.88 in the test set. This suggests whichhighlightstheInstatesinred). TheInstateshave that looking at oxides with larger unit cells could be loweffectivemasses(m =0.2−0.3)whiletheTastates beneficial[31, 66]. eff have much higher effective masses (m =13−62), al- DuetotherelativelyweakvariationintheDebyetem- eff lowing the material to take advantage of both types of peraturethroughoutoursetofoxides, theGruneisenpa- bands, inadditiontotheincreasedDOSprovidedbythe rameter becomes more important to identify the materi- near degeneracy. The large effective masses of the Ta- als likely to have low thermal conductivity. Many of the d bands are caused by small overlap between them and oxides we consider have soft or unstable phonon modes, many of the neighboring In-p orbitals, which results in whichlikelyresultsinstronganharmonicityandlowther- very flat bands, which contribute to a high DOS and malconductivity, butthisisdifficulttoquantifywithout high power factor. more involved calculations. Due to the large computa- Similarmaterialswithtwodifferentatomscontributing tional cost, we are only able to calculate the Gruneisen totheconductionareHgWO ,Ba TaInO andNbTl S , parameter of the materials in table II. We do not have a 4 2 6 3 4 where the Hg(+2), In(+3) and Tl(+1) ions, contribute large enough database of oxide Gruneisen parameters to empty s/p-bands at similar energies to the transition identify any trends which would predict which materials 10 will have soft modes without doing phonon calculations. effective mass, respectively, of the conduction band elec- As shown in table II, many of the materials we have trons. However, to find materials with high power fac- identifiedashavingpromisingpowerfactorsalsohavelow tors, it is necessary to look for materials which are not thermalconductivityaccordingtoourmodel. Mostofthe well described by a single parabolic band, but instead perovskite materials we study have strongly anharmonic have degeneracy, anisotropy, or other features which re- modes, which leads to relatively low thermal conductiv- sult in a high density of states combined with disper- ities, both in our calculations and in experiment[73]. In sive bands at the Fermi level. These materials achieve addition, we find that many of the materials with one or their high power factors due to some combination of twodimensionalbondingalsohavesoftmodes,likelydue symmetry-enforced degeneracies, low dimensionality, or to the fact that many of the atoms are relatively free to accidental degeneracies. In addition, we use phonon cal- vibrate in at least one direction. culations to model the thermal conductivity of our best candidates, and we find many that have low lattice ther- malconductivityorthatrequireanharmonicstabilization IV. CONCLUSIONS of the harmonic modes. We hope further work on these materials, as well as the understanding gained by exam- We have used high throughput first principles calcu- iningthemechanismswhichleadtohighpowerfactorsin lations to search for n-type transition metal oxides, ni- oxides,willleadtoimprovedthermoelectricperformance trides, and sulfides which are promising for thermoelec- in oxides. tricapplications. Wefindmanymaterialswithestimated power factors which are comparable to or surpass previ- ously studied oxide thermoelectrics. Acknowledgments Across the entire sample of compounds, we find the expected correlations between the Seebeck coeffient and We wish to acknowledge discussions with Igor Levin electricalconductivitywiththeeffectivemassandinverse andhelpwiththeICSDfromVickyKarenandXiangLi. [1] S.Walia,S.Balendhran,H.Nili,S.Zhuiykov,G.Rosen- cher,T.Mueller,K.A.Persson,andG.Ceder,Comput. garten, Q. Wang, M. Bhaskaran, S. Sriram, M. Strano, Mater. Sci. 50, 2295 (2011), ISSN 0927-0256. andK.Kalantar-zadeh,Prog.Mat.Sci.58,1443(2013). [13] A. Jain, G. Hautier, C. J. Moore, S. P. Ong, C. C. Fis- [2] M. Backhaus-Ricoult, J. Rustad, L. Moore, C. Smith, cher,T.Mueller,K.A.Persson,andG.Ceder,Comput. and J. Brown, Appl. Phys. A 116, 443 (2014). Mater. Sci. 50, 2295 (2011). [3] J.He,Y.Liu,andR.Funahashi,J.Mater.Res.26,1762 [14] G. Ceder, Y. Chiang, D. Sadoway, M. Aydinol, Y. Jang, (2011). and B. Huang, Nature 392, 694 (1998). [4] M.Zebarjadi,K.Esfarjani,M.S.Dresselhaus,Z.F.Ren, [15] A. R. Akbarzadeh, V. Ozolins, and C. Wolverton, Adv. and G. Chen, Energy Environ. Sci. 5, 5147 (2012), URL Mater. 19, 3233 (2007), ISSN 1521-4095. http://dx.doi.org/10.1039/C1EE02497C. [16] D. Morgan, G. Ceder, and S. Curtarolo, Meas. Sci. & [5] S. Riffat and X. Ma, Applied Thermal Engi- Tech. 16, 296 (2005), ISSN 0957-0233. neering 23, 913 (2003), ISSN 1359-4311, URL [17] O.Levy,G.L.W.Hart,andS.Curtarolo,J.Am.Chem. http://www.sciencedirect.com/science/article/ Soc. 132, 4830 (2010). pii/S1359431103000127. [18] G. H. Johannesson, T. Bligaard, A. V. Ruban, H. L. [6] M. W. Gaultois, T. D. Sparks, C. K. H. Skriver, K. W. Jacobsen, and J. K. Norskov, Phys. Rev. Borg, R. Seshadri, W. D. Bonificio, and D. R. Lett. 88, 255506 (2002). Clarke, Chemistry of Materials 25, 2911 (2013), [19] S.Curtarolo,D.Morgan,andG.Ceder,Calphad29,163 http://dx.doi.org/10.1021/cm400893e, URL (2005), ISSN 0364-5916. http://dx.doi.org/10.1021/cm400893e. [20] W.SetyawanaandS.Curtarolo,Comput.Mater.Sci.49, [7] I.Terasaki,Y.Sasago,andK.Uchinokura,Phys.Rev.B 299 (2010). 56, R12685 (1997). [21] X. Zhang, L. Yu, A. Zakutayev, and A. Zunger, Adv. [8] M.Ohtaki,T.Tsubota,K.Eguchi,andH.Arai,Journal Funct. Mater. 22, 1425 (2012). of Applied Physics 79 (1996). [22] J. W. Bennett and K. M. Rabe, J. Solid State Chem. [9] S. Ohta, T. Nomura, H. Ohta, M. Hirano, H. Hosono, 195, 21 (2012). and K. Koumoto, Applied Physics Letters 87, 092108 [23] S.Curtarolo,G.Hart,M.BuongiornoNardelli,N.Mingo, (2005),URLhttp://scitation.aip.org/content/aip/ S.Sanvito,andO.Levy,NatureMaterials12,191(2013). journal/apl/87/9/10.1063/1.2035889. [24] A.Belsky,M.Hellenbrandt,V.L.Karen,andP.Luksch, [10] T. Okuda, K. Nakanishi, S. Miyasaka, and Y. Tokura, Acta Crystallographica Section B 58, 364 (2002), URL Phys.Rev.B63,113104(2001),URLhttp://link.aps. http://dx.doi.org/10.1107/S0108768102006948. org/doi/10.1103/PhysRevB.63.113104. [25] S.Wang,Z.Wang,W.Setyawan,N.Mingo,andS.Cur- [11] W. Pan, S. Phillpot, C. Wan, A. Chernatynskiy, and tarolo, Phys. Rev. X 1, 021012 (2011). Z. Qu, MRS Bulletin 37, 917 (2012). [26] I. Opahle, A. Parma, E. J. McEniry, R. Drautz, and [12] A. Jain, G. Hautier, C. J. Moore, S. P. Ong, C. C. Fis- G. K. Madsen, New Journal of Physics 15, 105010