Optimal Band Structure for Thermoelectrics with Realistic Scattering and Bands

Understanding how to optimize electronic band structures for thermoelectrics is a topic of long-standing interest in the community. Prior models have been limited to simplified bands and/or scattering models. In this study, we apply more rigorous scattering treatments to more realistic model band structures - upward-parabolic bands that inflect to an inverted parabolic behavior - including cases of multiple bands. In contrast to common descriptors (e.g., quality factor and complexity factor), the degree to which multiple pockets improve thermoelectric performance is bounded by interband scattering and the relative shapes of the bands. We establish that extremely anisotropic"flat-and-dispersive"bands, although best-performing in theory, may not represent a promising design strategy in practice. Critically, we determine optimum bandwidth, dependent on temperature and lattice thermal conductivity, from perfect transport cutoffs that can in theory significantly boost optimized $zT$ to values exceeding 10. Our analysis should be widely useful as the thermoelectric research community eyes $zT>3$.


Introduction
Thermoelectricity enables clean electricity generation and fluid-free cooling. The ultimate goal of basic thermoelectric materials research is to design or discover materials with high figure of merit zT , commonly expressed as Here, the thermoelectric power factor (PF) is the product of Ohmic charge conductivity (σ) and the Seebeck coefficient (α) squared. The total thermal conductivity κ is the sum of electronic thermal conductivity (κ e ) and lattice thermal conductivity (κ lat ). A major challenge in achieving high zT and PF is that the electronic transport quantities are linked by a set of anti-complementary correlations: [1][2][3][4][5][6][7][8][9][10][11] σ and κ e are positively correlated whereas σ and α are negatively correlated. Only κ lat , a lattice property, is relatively independent, though it too exhibits some positive correlation with σ through structural symmetry. These interrelations make it difficult to determine the effect of various design strategies to optimize zT . Equations based on the single parabolic band (SPB) model often underpin intuition about thermoelectric behavior. However, they tacitly assume that there is always enough (infinite) dispersion in all directions to cover the entire energy range relevant to thermoelectric phenomena. Instead, a true solid-state band is bounded by the Brillouin zone (BZ) and has a finite dispersion. Also, in most cases of practical interest, a band's dispersion must change in curvature from positive to negative and cross the boundary orthogonally. In addition to band shape considerations, thermoelectric properties can widely vary depending on what is assumed of the scattering behavior. Typical models and descriptors assume a behavior that is dominated by intraband/intravalley, elastic acoustic phonon scattering, and can be derailed when other scattering mechanisms and interband/intervalley transitions have large effects [12][13][14][15][16]. Several studies have analytically investigated thermoelectricity using model band structures and scattering [15][16][17][18], but they had one or both of the following limitations: 1) the bands were purely parabolic or parabolic-like with infinite dispersion; 2) models for scattering and/or transport were based on constant lifetimes, constant mean free paths, or at best scattering proportional to the density of states (DOS).
In this study, we create more realistic solid-state model band structures and more faithfully model carrier scattering due to multiple sources. It thereby not only generally addresses the question of optimal band structure but also fine-tunes some conclusions drawn from the simpler models. We generate model band structures properly confined to a finite BZ by smoothly inverting upward parabolic bands on the way to the BZ boundary -a key for retaining both generality and approximate compatibility with established scattering formalism. The main objective is to monitor how thermoelectric properties of one or more bands respond to variation of directional effective masses of a band (see Fig. 1). We fix the band mass of the principal band in the transport direction (m x ) -except when the principal band itself evolves isotropically -to compare the performance of various band structures on an equal footing. We apply established formulas for various scattering mechanisms, modified for our band structures, to compute carrier lifetimes under deformation-potential scattering (DPS), polar-optical scattering (POS), and ionized-impurity scattering (IIS). We examine thermoelectric behavior of a variety of band structures and design  Note that two-dimensional (2D) band structures are shown for graphical purposes. In the study where three-dimensional (3D) bands are used, two types of anisotropic evolution are considered: one where a band grows heavy in one direction and another where the band grows heavy in two directions. Each band is an upward paraboloid smoothly inflecting to an inverted paraboloid halfway to the BZ boundary.
strategies, including anisotropic bands and resonance levels. Finally, we determine the optimum bandwidths as a function of temperature and κ lat , which can critically improve zT beyond what is typically accessible.

Results
To set the stage for the discussions to follow, we rewrite Eq. 1 as to better reflect fundamental physics: In this formulation, the key role is played by ζ, a quantity for which there appears to be no conventional name. We refer to it as the "thermoelectric conductivity"; in the Onsager-Callen formulation of coupled charge-and-heat conduction [19][20][21], ζ is the quantity responsible for the thermal-gradient-to-charge-current conversion. That is, ζ represents the charge conductivity due to thermal driving force (J c = σE − ζ∇T ), which is the essence of thermoelectricity. An immediate benefit of Eq. 2 is that it removes the hidden coupling between α and σ (α = ζ/σ). This correctly identifies ζ as the quantity that must be high and σ as a quantity that must be low, counter to conventional thermoelectric interpretation. That is, we desire high thermoelectric conductivity, not Ohmic conductivity -a correction to the routine but ambiguous thermoelectric adage that "electrical conductivity" must be high.
To identify the factors for optimal zT , Eq. 2 must be examined through the lens of fundamental transport equations. Under Boltzmann transport formalism [22][23][24][25], where V is the cell volume, E F is the Fermi level, f (E) is the Fermi-Dirac distribution, and Σ(E) = v 2 (E)τ (E)D(E) is the spectral conductivity, composed of group velocity (v), lifetime (τ ), and DOS (D). The three integrands share in common the term Σ(E) − ∂f ∂E , the source of the positive correlations between σ, ζ, and κ e . The integrands differ only in the power relation (E F − E) p as p = 0, 1, 2. This juxtaposition states that, in relative terms, low energy carriers contribute most to σ, highenergy carriers contribute most to κ e , while it is the medium-energy carriers that are most responsible for ζ. That is to say, if one wished to increase ζ relative to σ and κ e , then Σ(E) should be high in some medium-energy range and low elsewhere. The results that follow are interpreted with this picture in mind.
In the following, we examine the thermoelectric behaviors of various band shapes as directional effective masses a b  Fig. 1a. b) Two bands exist where the second band evolves as depicted in Fig. 1b. Each zone indicates certain characteristic evolution: isotropic increase in m from 0.05 to 500 in Zone 1, anisotropic increase in my from 0.05 to 500 in Zone 2, anisotropic increase in both my and mz in Zone 3 from 0.05 to 500. Four different scattering regimes are considered: the POS limit (blue), the IIS limit (green), the DPS limit (red), and the overall effect (black). The black dashed horizontal line marks the isotropic value.
are varied as depicted in Fig. 1. For zT , we use low κ lat of 0.5 W m −1 K −1 that is representative of highperformance thermoelectrics. Details of the scattering models used to calculate the thermoelectric properties for these band shapes are presented in the Methods section.
The Seebeck Coefficient with a Fixed Fermi Level. To begin, we focus on the behavior of the Seebeck coefficient α with a fixed E F at the band minimum. In later sections, we examine the thermoelectric behavior with E F optimized for maximum zT . However, we begin here as it helps clarify the difference in behavior of our model versus previous models, and also clearly illustrates the concept of improving ζ relative to σ via band structure alone. The results are plotted in Fig. 2, where three "zones" are displayed: 1) isotropic increase in m, 2) anisotropic increase in m in one direction ("unidirectional anisotropy"), and 3) anisotropic increase in m in two directions ("bidirectional anisotropy").
The results for single band evolution as depicted in Fig.  1a are given in Fig. 2a. The main takeaway from Zone 1 is the problem of insufficient dispersion for critically heavy bands. We observe that α starts on a plateau that corresponds to the classic SPB model behavior. When the band becomes critically heavy, however, α starts to decrease below the SPB value. The origin of this deviation is that the band becomes heavy enough to terminate at the BZ boundary before gaining enough energy to fully trigger ζ, and misses out on some high-energy states that would otherwise contribute relatively more to ζ than σ. Typical SPB models overlook this problem of insufficient dispersion owing to the finite-sized BZ, breaking down for heavy enough dispersions. More discussion of SPB models is provided in in supplemental section IIA [26]. Of note, α goes to 0 as the band completely flattens out, which is explained in supplemental section IIC.
The main lesson from Zones 2 and 3 in Fig. 2a, where the band evolves anisotropically, is the role of group velocity. Under DPS and to a lesser extent under POS, we observe that moderate anisotropy gives the highest α, represented by the α peak in the middle of the two zones (m y , m z = 5). Because Σ(E) ∝ v 2 (E) under DPS, thermoelectric behaviors are determined entirely by the average group velocities, v 2 x (E) . For an isotropic parabolic band, v 2 (E) = 2mE 3 ; thus, the contribution to Σ(E) scales linearly with E. For a moderately anisotropic band, v 2 x (E) develops a kink. That is, it abruptly steepens in slope (see supplementary Figs. S6a-b [26]). Specifically, for unidirectional anisotropy (heavy only in z), v 2 x (E) ∝ mE post-kink which is the 2D parabolic velocity scaling. For bidirectional anisotropy (heavy in y and z), v 2 x (E) ∝ 2mE post-kink which is the even steeper 1D parabolic velocity scaling. Relative to the scaling of an isotropic band, the kinked v 2 x (E) profiles weight ζ more than σ because the velocities increase more steeply at higher energies that at lower energies. This allows α to peak at some moderate anisotropy, and the peak is higher for bidirectional anisotropy. For extremely anisotropic bands mimicking low-dimensional bands, which are also popularly known as "flat-and-dispersive" bands [27,28], v 2 x (E) reverts to linear scaling but with the steeper, post-kink slopes.
Another subtlety regarding extremely anisotropic bands is that they exhaust "low-energy voids". That is to say, wherever carriers line up along the heavy direction(s), their dispersion in the light direction starts from essentially the band minimum energy. This poses a stark contrast to a less anisotropic band or an isotropic band, for which some dispersion towards the light direction may start from higher energies leaving behind a void of states at lowest energies (see supplemental Fig. S6c). Because low-energy states contribute relatively more to σ than to ζ, their absence is a clear benefit to thermoelectrics. Ex-tremely anisotropic bands exhaust these voids. Therefore, the overall α for extremely anisotropic bands is somewhat lower than the isotropic case. POS retains some of the same α signatures of DPS, while under IIS α only decreases with anisotropy.
Next, we examine the case of two bands whose results are in Fig. 2b. Here, we fix the first band in shape and evolve the second band mass according to Fig. 1b. The results are largely similar to the single-band results but for a prominent peak in α in the middle of Zone 1 (under DPS and POS) where the second band flattens out isotropically. This peak represents the second band acting as a resonance level [29] that performs energy-filtering due to interband scattering. Although the isotropically heavy band has negligible direct contribution to transport, it can act as a localized scattering partner for the dispersive principal band where their energies overlap, or "resonate," thereby preferentially scattering low-energy carriers. This increases ζ relative to σ because the low-energy states that had previously contributed much more to σ than ζ are now selectively scattered by the narrowed second band. As a result, α is able to well exceed its singledispersive-band value. As the second band completely flattens out, however, its width becomes too narrow to filter enough states and hence α is reduced again. We observe that IIS is not a good agent of energy-filtering.
In summary, whereas α is constant for any band under the SPB model as long as E F is fixed, our revised model correctly reflects its fluctuating response to changes in a band structure especially as it approaches extreme shapes.
Optimal Performance -Single Band. We now attempt to identify the form of a single band that is ultimately capable of yielding the highest zT with E F optimized for each case to extract their optimal zT . The optimized E F is plotted in Fig. 3a. E F optimizes entirely below the band minimum, i.e., below zero. This can be understood by noting that generally κ e >> κ lat in our model, which in turn leads to zT ≈ α 2 L where L is the Lorenz number. Because α is higher at lower E F and L is relatively constant with respect to E F , zT generally favors low E F and non-degenerate doping. However, because κ e is lower for heavier bands, higher PF is required for optimizing zT . Because the PF is maximized when E F is near the band minimum, E F increases to access it as a band turns heavier. However, as the band becomes critically heavy and narrow, E F must again be placed at a distance from the band minimum in order to generate finite α for reasons described in supplemental section IIC.
We first consider the case where the band mass varies isotropically (Zone 1 in Figs. 3b-d). As expected, a light band is preferred: the PF and zT decrease with increase in m, as numerous studies agree upon [30][31][32]. The key is that a lighter band has higher mobility (µ) and thus is less needy of carrier concentration (n) in providing a given value of σ (σ = nµ). Because adjusting E F to raise n reduces α, a heavier band is at a disadvantage. For a parabolic band, µ ∝ m − 5 2 ∼− 3 2 (the exponent depends on scattering) whereas n ∝ m 3 2 , meaning σ ∝ m −1∼0 . As a result, performance degrades with increasing m.
Anisotropic bands are deemed beneficial primarily for their capability to accommodate, along the heavy direction(s), a very large number of carriers with high velocity in the light direction, for given E F . However, increase in n via anisotropy comes at a price of enlarged phase space for scattering. In fact, if scattering is exactly proportional to DOS as in DPS, then the gain in n ∝ D(E) is exactly canceled out by decrease in τ ∝ D −1 (E), forcing the benefits of anisotropy to draw purely from the behavior of the average group velocity in the transport direction, v 2 x (E) . Our results for anisotropic band evolution are plotted in Zones 2 and 3 in Fig. 3.
We observe that, even under DPS, anisotropy is immensely beneficial upon optimizing E F . This is rooted in the aforementioned steepening slope of the v 2 (E) profile of anisotropic bands that favors ζ over σ, which simultaneously lowers optimal E F . We make three major observations. First, in terms of zT , bidirectional anisotropy (one light, two heavy directions) outperforms unidirectional anisotropy (two light, one heavy direction) due to steeper v 2 (E) . Second, toward the extreme limit, both types of anisotropy plateau in performance. This occurs for the two reasons described in the previous section: v 2 (E) under extreme anisotropy converges to the low-dimensional linear limits, and extreme anisotropy exhausts low-energy voids. Third, because IIS and POS are less dependent on D(E) than DPS, anisotropy is even more beneficial under them. Eventually though, owing to the closest relation to D(E) of DPS, it inevitably overtakes POS and IIS as anisotropy grows larger; the overall effect thus most closely resembles the DPS limit.
Though anisotropic bands can theoretically deliver dramatic performance, with zT nearing 10 being modeled for bidirectional anisotropy, we later discuss some practical issues with anisotropy in the Discussion section.
Optimal Performance -Multiple Bands. Realistic band structures often feature multiple bands near the Fermi level. Conventionally, one of the best designs for increasing σ without paying a penalty on α is the presence of multiple band pockets aligned at or near E F [33][34][35][36][37]. Band multiplicity comes in various forms, however; we therefore examine the effects of (i) multiplicity of identical bands, (ii) coexistence of inequivalent bands (with varying the second band shape), and (iii) bipolar transport in the presence of valence and conduction bands (with varying valence band shapes). These band structures are illustrated in Fig. 1b-c. We remind the reader (see Methods) that our modeling of interband/intervalley scattering (henceforth inter-scattering) is rather phenomenological, defined only through the enlarged phase space and the factor s int = 0.5 making it half as strong as intraband scattering (henceforth intra-scattering). Nevertheless, it is sufficient for broadly describing the effects of the above cases. For comparison purposes, we also provide results obtained without inter-scattering, s int = 0, in supplemen- tal Fig. S8. In practice, the magnitude of s int will depend on multiple factors, such as wavefunction symmetries and the relative locations of band pockets in the BZ.
We generate two bands with aligned band minima, the first of which is isotropic and fixed while the second band starts identically but then evolves, according to Fig. 1b. E F is again optimized for zT . The results are plotted in Fig. 4. The left edge of Zone 1 for each plot, where the two bands are identical, represents symmetry-degenerate band pockets. This offers higher zT and PF as compared to the case of a single band (Fig. 3) though less than by twofold. Two identical bands result in essentially identical α as σ draws benefits from doubled n somewhat negated by inter-scattering. One question of interest is the effect of increasing the number of identical carrier pockets. It is generally known that the more pockets the better, though it is straightforward even from our simplified analysis that doubling their number does not double the PF or zT due to inter-scattering.
, which as N v grows saturates to s −1 int . For example, with s int = 0.5 that we assume, the maximum PF gain even with an infinite number of identical pockets is a factor of 2. In fact, if inter-scattering is stronger than intra-scattering, s int > 1, then N v is detrimental. As such, the benefit of increasing pockets is bounded by the degree of inter-scattering, and minimizing the degree of inter-scattering should be a priority of multipocket strategies. Furthermore, in the κ e >> κ lat regime, the relative benefit of additional pockets for zT is even more limited because more pockets would provide a higher κ e as well as a higher PF (through σ). If κ lat = 0, hypothetically, then N v would have no effect on zT at all because it will cancel exactly for the PF and κ e .
Next, keeping the principal first band fixed in shape and maintaining s int = 0.5, we make the second band heavier. As it turns heavier isotropically (Zone 1 in Fig.  4), zT and the PF increasingly suffer until they sink well below even the values that the lighter principal band alone generates (Fig. 3). This means that non-symmetryrelated, accidentally degenerate pockets harms zT if their band masses in the transport direction (m x ) are sufficiently different. Two main reasons account for this. As the second band grows heavier isotropically, its direct contribution to transport diminishes, and it indirectly sabotages the lighter principal band by triggering heavier inter-scattering within the energy range relevant to thermoelectric transport. That is, until the second band's width becomes narrow enough as to function as a reso- nance level, whereby zT and the PF rebound. They do not fully recover the values generated by the original twin degenerate bands unless DPS or POS dominates. The presence of strong IIS, due to the high impurity concentration required for doping a very heavy band, could limit the resonance level effect. If the second band evolves anisotropically in the y and/or the z directions, the thermoelectric response is largely similar to what is seen for a single band turning anisotropic. Anisotropy increases α and the PF as well as zT until they plateau. Also, zT is not noticeably higher here than in the case of single anisotropic band because the anisotropic band dominates transport and κ e > κ lat . This again is a nod to the increasing unimportance of the presence of additional bands as κ e > κ lat .
Another two-band situation we model is a semimetallic one in which there exists a "conduction band" and a "valence band" with no gap in between, triggering bipolar transport. Bipolar effect is a significant suppressor of the Seebeck coefficients of metals and small-gap semiconductors. Extrapolating the lessons from above, it is rather straightforward that for ζ to be large in magnitude (whether positive or negative), Σ(E) must be highly asym-metric about the Fermi level: either juxtaposing mobile and anisotropic "conduction" bands against isotropically heavy "valence" bands or vice-versa. We confirm this by fixing the conduction band and evolving the valence band as described in Fig. 1c. Detailed results and discussions are provide in supplemental section IIB. It is therefore no surprise that high-performing semimetals and narrow-gap semiconductors feature quite drastic band asymmetries about the Fermi level [38][39][40][41][42][43].
Optimum Bandwidth from Perfect Transport Cutoff. The relative energy ranges from which σ, ζ, and κ e draw contributions imply that the best performance would be obtained by suppressing both low-energy contributions (to suppress Ohmic current and improve α) and veryhigh-energy contributions (to suppress thermal current and lower L). Accordingly, we investigate the scenario in which the contribution to transport abruptly vanishes at some optimum energy, which we define as the optimum bandwidth (see the inset in Fig. 5a). It would arise for a band that is abruptly crossed by numerous, flat, energyfiltering states acting as perfect resonance levels, or a band that sharply and discontinuously flattens out at some energy. Admittedly neither is achievable to perfection in real life, but that theoretical limit is of our interest.
We stress that our definition of bandwidth differs from the "full-width" definition adopted in a previous study [18] whereby bandwidth was bounded by effective mass and the BZ boundary (increasing width with decreasing m). In that study of tight-binding dispersions, it was determined that there exists no optimum full-width under the τ (E) ∝ D −1 (E) model as zT continues to increase with larger full-width (smaller m) [18]. Our bandwidth, in contrast, essentially represents the finite transport distribution width and is not inherently bound to band mass, although the optimum bandwidth (W opt ) may still depend on it. Mathematically, we solve for the following as to maximize zT : We consider an isotropic parabolic band under DPS only and optimize E F . Within this set-up, we find that there exist temperature-and-κ lat -dependent W opt delineated by Fig. 5. Achieving a perfect energy filter at W opt would be a tremendous boost for zT , capable of elevating it well beyond 10 -higher than any value attainable through plain band structures of previous sections. For given κ lat , Fig.  5a shows that W opt generally increases with temperature, as expected as the larger range of carrier excitation at higher temperatures. This implies that achieving W opt is particularly consequential for low temperatures (T ≤ 300 K) where a difference of 0.1 ∼ 0.2 eV can force a shift in zT by an order of magnitude beyond 5.
As κ lat decreases to 0, W opt also shrinks to 0, and zT diverges. This would be the Mahan-Sofo limit, so named after the seminal paper in which they deduced that a fully localized, widthless band would deliver maximal thermoelectric performance [44], if only κ lat = 0. Our recovery of this limit is also evidenced by supplementary Fig. S4f, in which the electronic-part zT diverges to infinity as the band completely flattens out. Because κ lat > 0 in real materials, a widthless band would in practice yield zT = 0 as v(E) and the PF vanish alike. At the other end of the spectrum, as κ lat increases to very large values, W opt also diverges for all temperatures such that the bandwidth becomes virtually irrelevant for zT .

Discussion
Theoretical Interpretations. Commonly used descriptors for thermoelectric efficiency include the quality factor, which under DPS is [45] and the Fermi surface complexity factor [32] Both metrics promote small effective mass (m or m ) and high band multiplicity (N v ); the latter further promotes band anisotropy m ⊥ m . This study serves as a general assessment of these well-known blueprints in thermoelectrics, confirming some while offering fresh perspectives and more complete physical pictures to others.
1. Small m in the transport direction is always beneficial. Band anisotropy is also very beneficial, but the extent depends on its type. The advantage of anisotropy draws largely from the fact that v 2 x (E) rises to a steeper, low-dimensional slopes. Bidirectional anisotropy mimicking 1D band structure is particularly beneficial, capable of increasing maximum zT by more than twofold for a given band mass in the light direction.
2. Although not captured by Eq. 7 or Eq. 8, the gains from pocket multiplicity and band convergence depend on the relative shapes of the bands and what is assumed of interband scattering. A heavier pocket in the presence of a lighter pocket can be detrimental. Metrics such as QF or C always predict better performance in the presence of more bands because they do not possess any component that accounts for inter-scattering or differential intrinsic transport of each bands. The metrics ought take these effects into account by bounding the gain from N v , e.g., using a term such as N v (1 + s int (N v − 1)) −1 as was previously described.
3. Within the limits of our investigation, the type of scattering mechanism does not play a pivotal role in determining what band structure is optimal for zT , except in the context of resonance levels. In other words, the best-performing band structure is for the most part the same under the DPS, POS, or IIS. The type of scattering decides how much zT increases as a band transforms, but no transformation is decisively beneficial under one scattering regime but decisively detrimental under another. Exceptionally, resonance levels can be very beneficial if the dominant scattering mechanism is efficient at energyfiltering -DPS or POS in our model. If an ineffective filtering mechanism, such as IIS, activates comparably or dominates, then resonance levels lose merit.
4. There exist optimum bandwidths to a plain parabolic band if the transport contribution can be, albeit hypothetically, abruptly curtailed at some energy. Optimum bandwidth arises because low-energy states are undesired owing to their large contribution to σ and very high-energy states are undesired owing to their large contribution to κ e . W opt for a given dispersion depends on temperature and κ lat ; it is typically small (< 0.3 eV) so long as κ lat < 1 W m −1 K −1 , and can push zT beyond 5 at low temperatures and beyond 10 at high temperatures.
5. Though more of a philosophical point, we propose that analysis of thermoelectrics be more frequently framed in terms of the "thermoelectric conductivity," ζ, which can offer more straightforward insights than that those framed in terms of the Seebeck coefficient and Ohmic conductivity. By juxtaposing ζ against σ and κ e , it becomes clear that a band must develop high Σ(E) in the mid-energy region to be optimal for thermoelectric application.

Flat-and-Dispersive Band Structures in Real Materials.
In spite of the theoretically remarkable performance of extremely anisotropic, flat-and-dispersive band structures, they in practice would be subject to a disadvantage due to symmetry considerations and polycrystallinity of real commercial-scale materials. Indeed, candidate materials have thus far not achieved the high zT modeled here, and we project why in light of our modeling. We distinguish bands in a cubic cell from those in a non-cubic cell.
1. In a non-cubic cell, a flat-and-dispersive band limits light transport to only certain direction(s). Assuming polycrystallinity, conductivity through a series of differently oriented grains is best described by the lower Wiener bound for composite media, i.e., the harmonic average of directional conductivities [46]. Due to poor conductivities along the heavy branch(es), the harmonic average seriously hampers the overall performance under our model, as described in supplementary Fig. S5. An anisotropic band is then never as good as its isotropic counterpart whose polycrystalline-averaged conductivities are identical to those in any principal direction.
Bi 2 PdO 4 [47] and BaPdS 2 [48] are good examples as neither compound is cubic but exhibit bidirectionalanisotropic flat-and-dispersive valence bands. The DOS profiles are characterized by peak-like protrusions near the band edges followed by decays, confirming the 1D-like band structure. Polycrystalline Bi 2 PdO 4 has been experimentally synthesized and investigated, but recorded rather disappointing p-type PF (1 mW m −1 K −2 ) and zT (0.06) [49]. Because electronic transport is mobile only in one direction and inhibited in the two heavy directions by design, it is unlikely that their presumably high thermoelectric potential in the light direction would shine through unless the sample is a single crystal. BaPdS 2 has not yet been tested, but it is reasonable to hypothesize that it may exhibit similar behavior.
2. Under cubic symmetry, all three principal directions are guaranteed the same number of light and heavy branches. Polycrystallinity may be irrelevant here, but the coexistence of light and heavy branches in the direction of transport (as opposed to light only in transport direction, heavy only in other directions) with interscattering between them would be inherently detrimental (recall the multiband discussion and Zone 1 in Fig. 4). Relevant cases are Fe-based full-Heuslers [50] and perovskite SrTiO 3 [51], which exhibit unidirectional flat-anddispersive conduction bands, confirmed by the 2D-like, precipitous DOS at the band edge.
Fe 2 TiSi, a member of the former family, is particularly intriguing because its flat-and-dispersive conduction bands are opposed by triply-degenerate, isotropically dispersive valence bands, offering a direct comparison of the respective n-type and p-type performances of this band structure. According to our DFT [52] band structure calculation using the PBE functional [53], the lowermost conduction band is very flat along Γ − X whose energy width is 0.05 eV (m ≈ 41) and dispersive in other directions (m ⊥ ≈ 0.7). A second isotropic conduction band (m ≈ 0.4) is degenerate at Γ. Opposing them are three isotropic valence bands with comparable 0.4 ≤ m ≤ 0.75. Theoretical thermoelectric properties were studied with rigorous first-principles treatment of electron-phonon scattering [54], but the n-type PF (5 mW m −1 K −2 at 300 K) was predicted to be only barely higher than the p-type PF (4 mW m −1 K −2 at 300 K), with no sign of the dramatic performance promised by the flat-and-dispersiveness.
As for SrTiO 3 , according to our DFT band structure calculation, the width of the heavy branch of the lowermost conduction band along Γ − X is approximately 0.1 eV (m ≈ 7) while in the dispersive directions m ⊥ ≈ 0.8. Two additional relatively isotropic conduction bands (0.4 ≤ m ≤ 0.7) disperse from Γ at the CBM. Multiple experimental reports exist for SrTiO 3 on single crystals, which should be the best-performing and the most comparable to theory results. Although respectable n-type PFs of 3.6 mW m −1 K −2 [55] and 2.3 mW m −1 K −2 [56] have been recorded at room temperature, neither are these values anywhere near what Fig. 3 promises.
These observations collectively suggest that cubic symmetry indeed limits the full potential of the flat-anddispersive bands in real materials. As a separate point, it would help if the band masses in the light direction of the both compounds were much smaller.
The Optimal Band Structure for Thermoelectrics. As a final deliberation, we address the question, what then is the optimal band structure, all things considered? The literature has convincing cases for both an extremely anisotropic, flat-and-dispersive band and a band with multiple dispersive pockets at off-symmetry points. Reflecting on our modeling, we conclude the following. For a single band, if a bidirectional flat-and-dispersive band attaining m ⊥ m ≈ 1000 but with m < 0.1 can be realized in a single-crystal, it would constitute the optimal single band structure. Otherwise, a band with multitude of dispersive pockets at off-symmetry points with weak inter-scattering would be the best targets, as they provide moderate anisotropy and would be more immune to polycrystallinity. For a given single band, presence of additional bands with equally light mass in the transport direction would increases performance, though this benefit is limited in the κ e < κ lat regime and/or if s int ≈ 1. For a given overall intrinsic band structure, engineering resonance levels and optimum bandwidth would further improve performance, the latter being capable of boosting zT to unprecedentedly high values and particularly consequential for low-temperature operation.

Conclusion
As efforts to discover and design thermoelectrics with zT > 3 continue, the blueprints for high performance grow increasingly influential. Common rules regarding beneficial band structures for bulk thermoelectrics are largely drawn from simple band models without realistic scattering. Using a straightforward but improved approach, we herein fine-tune those blueprints while proposing optimal band structures and design principles along the way. Our generalized findings from this modeling study are mutually supportive of and consistent with the findings from recent targeted studies of high-performing materials with high-fidelity first-principles computations [38,40,54,57,58]. We hope the theoretical investigations of the present study help the community navigate rationally towards next-generation thermoelectrics.

Methods
The Hartree atomic units ( = m e = a o = q = 4π 0 = 1) are used throughout the methods section. All calculations are performed on a set of in-house Mathematica codes.
Band Structure. To generate a realistic solid-state band structure, a band is created by smoothly connecting an upward paraboloid to an inverted paraboloid at a selected inflection point. The advantages of such a band structure include: 1) it remains faithful to solid-state band theory that requires a band to cross the zone-boundary orthogonally, save for when crystal and orbital symmetries allow band crossing or degeneracy at the zone boundary (e.g. in graphene); 2) relatively simple analytic models of scattering can be directly applied to the upward-parabolic portion, and can be applied with modification for the inverted-parabolic portion; 3) it formally maintains the validity of effective-mass-based description of the band throughout. The equation for this band structure is Here,k x denotes the inflection point in the x-direction, and G x is the reciprocal lattice vector in the x-direction, i.e., the BZ boundary in the x-direction. If inflection occurs at halfway to the BZ boundary, thenk x = G x /2.
The last terms are subtracted (−) ifk x ≥ G x /2, and added (+) ifk x < G x /2. Effective masses of the inverted, downward parabolic portion (m down ) are obtained by enforcing derivative continuity at the inflection point in every direction. Although a broad range of band shapes could be explored by changing the inflection pointk in any three Cartesian directions, in this work we limit ourselves to bands that inflect halfway to the zone boundary in all three directions. Under this assumption, the effective mass (inverse curvature) profiles of the upward paraboloid portion and the inverted paraboloid portion are identical (m up = m down ), rendering the entire band structure describable with one common set of directional effective masses. We note that this band could also serve as a firstorder approximation of the tight-binding cosine band. We create these model bands centered at Γ in a simple cubic Brillouin zone corresponding to an arbitrary lattice parameter of 15 a o (Bohr radius) ∼ 7.9Å, which is a reasonable lattice parameter for a real thermoelectric. The primary parameters to be tuned in this study are the directional effective masses, or dispersions, of the band(s), which result in different types of electronic structure. The band(s) may evolve isotropically or anisotropically and may be offset in energy. The types of band structures investigated in this work are illustrated in Fig. 1. The density of states (DOS) is calculated using the tetrahedron method [59]. Details on the band structure and DOS are given in the Supplementary Information [26].
Carrier Scattering and Transport. Thermoelectric properties are computed by numerically integrating Eqs. 3-5. The BZ is sampled to convergence with a k-point mesh of 40 × 40 × 40. To calculate τ , various scattering mechanisms, namely deformation-potential scattering (DPS) by (acoustic) phonons, polar-optical scattering (POS), and ionized-impurity scattering (IIS), are treated according to well-established formalisms [23,[60][61][62][63][64] with appropriate adjustments to account for non-parabolicity and the BZ-bounded nature of our band structures, with details to follow in the next subsections. Once the lifetimes under the three mechanisms are calculated, the overall τ is estimated by Matthiesen's rule [65].
Deformation-Potential Scattering. Within the theory of deformation potential [66] and the assumption of elasticity for long-wavelength acoustic phonon scattering, the carrier lifetimes can be cast in the k-dependent form as follows: where ∆ is the deformation potential. Here, the energydependence is encoded purely in the DOS term. This describes intraband/intravalley scattering. When more than one bands/pockets are present, interband/intervalley scattering must be accounted for. A general expression for inter-scattering matrix elements is not obtainable. However, it is reasonable to assume that inter-scattering is generally weaker than intra-scattering because it relies on zone-boundary phonons (which are less populated than zone-center phonons) and because the wavefunction overlap between distinct bands or pockets is generally weaker than that within a band. To account for the strength of inter-scattering separably, we introduce a parameter that modulates its strength relative to intraband scattering, s int , and set s int = 0.5 to reflect the usually weaker inter-scattering. Though admittedly a phenomenological treatment, it does capture the essence of interband DPS. This allows for trivial extension of Eq. 10 by accounting for the added phase space due to the second band proportional to its DOS, and the overall DPS lifetime, for band 1, now becomes where the subscripts 1 and 2 indicate each of two bands. The addition of inter-scattering phase space is determined by the presence of second band states at given given E 1 . That is, if D 2 (E 1k ) = 0, then there is no inter-scattering. We note that because our DPS lacks the inelastic component, heavy bands or branches in this model have diminished effect on scattering of the rest of the band structure.
Polar-Optical Scattering. For a parabolic band, there exists an analytic formula for energy-dependent lifetime due to inelastic POS [60,61], The first (second) term in the square brackets represents emission (absorption). However, derivation of Eq. 12 uses an energy-momentum phase-space integral simplified with a parabolic dispersion relation, rendering its direct application to non-parabolic bands unjustified. The hyperbolic arcsin terms account for the availability of DOS (∼ √ E for a parabolic band) for carriers to be scattered into (final states). Important corrections to Eq. 12 can be made by using its k-dependent prior form and exchanging √ E with our tetrahedron-integrated D(E): where use of the group velocity norm (|v k |) obviates the dependence on the effective mass m and energy E through the relation |v k | = 2E m . The optical phonon frequency is represented by ω o , and b(ω o ) is the Bose-Einstein population. The left term in the square brackets accounts for phonon emission whereas the right term accounts for phonon absorption. For our band structures, Eq. 13 is exact for the upward-parabolic portions up to the inflection point, and past it, closely approximates the true lifetimes. When more than one band pocket exists at one k-point, nothing prohibits interband POS from occurring. To account for this, we take the same approach as we do with DPS and use D(E) = D 1 (E) + s int D 2 (E) to enlarge the phase space.
Ionized-Impurity Scattering.The established energydependent formalism for IIS of a parabolic band is the Brooks-Herring formula [63,64], where γ k is a screening term defined as, and F z (E F ) is the Fermi-Dirac integral Its derivation involves k-space integration over spherical isoenergy surfaces of a parabolic band, rendering it also non-trivial to extend to non-parabolic bands. However, much of the necessary corrections can again be made by using a k-dependent prior form: where N i is the concentration of ionized impurities, and Equations 17-18 are then evaluated using our tetrahedronintegrated DOS. For our band structures, this corrected formula is again exact for the upward-parabolic portions  from Γ to the inflection point, and past it, closely approximates the true lifetimes. When more than one band pocket exists at the same k-point, interband IIS can take place. To account for this, we again use to modify the phase space. Of note, we assume that one impurity donates or accepts one charge, meaning the effective impurity charge of Z = 1. This choice forces that the carrier concentration (n) effectively equals the impurity concentration (n = N i ) at appreciable doping levels. Fig. S3 in Supplementary Information plots the scattering and lifetime profiles of the three processes based on Eqs. 10, 13, and 17. For the isotropic cases, they exhibit excellent agreement with the corresponding parabolic band formulations. DPS exhibits the expected τ −1 ∝ D(E) relation. Polar-optical scattering is characterized by the emission onset near the band edge followed by a gradual decay. Ionized-impurity scattering faithfully reproduces the Brooks-Herring behavior.

Materials parameters.
To calculate specific values of scattering, material-dependent quantities that control the relative strength of the three scattering mechanisms must be chosen. These include the deformation potential, the dielectric constants, and the optical phonon frequency among others. We select plausible values for these quantities that occur prevalently in real materials, as listed in Table I. In emphasis, the choice of these values renders the relative strength of each scattering channel arbitrary. What is not arbitrary is the characteristic thermoelectric behavior of bands under a given scattering regime.

Data Availability
The Mathematica notebooks used for computation for this study will be uploaded online with the manuscript.

A. On the Single-Parabolic Band Model
A review of the single parabolic band model and the Pisarenko formulas for the Seebeck coefficient is in order owing to their generality, wide reference, and limitations our study improves upon. The Seebeck coefficient is, in the degenerate or metallic limit, and in a non-degenerate case, where n is carrier concentration, m is the band effective mass, k B is the Boltzmann constant, and s is the power of energy to which carrier lifetimes are proportional (τ ∝ E s ). Contrary to a widely held notion, a heavier band (high m) does not by default generate higher Seebeck coefficient under this model since α can be written as a function of either E F or m-and-n, where the trivial exchange of variables is allowed by their interrelations: in the degenerate case, and in the non-degenerate case. Equations S1-S2 stipulate that, if E F is fixed, then α is constant whatever the m value because n would change accordingly. A light band with a small m would produce the same Seebeck coefficient as a heavier band with a larger m provided that E F is kept fixed.
The lead-up to Eqs. S1-S2 bears the assumption that there is always enough (infinite) dispersion in all directions to cover the entire energy range relevant to thermoelectric phenomena. However, it breaks down when a band becomes critically heavy that it reaches the BZ boundary before gaining enough energy to cover the entire relevant energy range for transport. Also, because a true solid-state band must at some point change in curvature from positive to negative and cross the BZ boundary orthogonally, there comes a point where assuming constant positive curvature introduces additional unrealistic effects. These effects lead to some conclusions that deviate from what would otherwise be drawn from Eqs. S1-S2. Other assumptions in typical parabolic band models include absence of opposing bands (bipolar effect) and a monotonically behaved, or at least slow-varying, Σ(E) such that the Sommerfeld expansion is valid, a requirement for arriving at Eqs. S1-S2.

B. On Metals & Bipolar Effect
For a simple study of bipolar effect, we fix the band gap to 0 such that the two bands are if there is an asymmetry in Σ(E) about the Fermi level. The greater the contrast, the better, though the benefit effectively saturates past a point. From the n-type point of view, a completely flat valence band would not carry any hole current but only function as a potential resonance level for low-energy electrons (would require inelastic processes). The desired effect in this picture is essentially a hybrid of the light-band-over-heavy-band rule and energy-filtering: keep holes heavy and filter them further out with as much resonance scattering as possible, while keeping electrons mobile and scattering-free except at very low energies. It is worth recognizing that, in the limit of completely flat valence band, one essentially has a degenerate-semiconducting state with a resonance level and a "gap" below.
This indicates that, within the conventional setting herein assumed, the ideal limit for a metallic thermoelectric is precisely the semiconducting limit with only DPS being present. and S9a proves the point.
Lastly for metals and tiny-gap semiconductors, κ e is frequently far higher than κ lat , where small Lorenz number (L = κe σT ) becomes critical. Even for typical semiconductors, once κ lat is reduced and the PF is improved, realizing small L would be the final piece of the puzzle.
Mirroring the way in which bipolar transport is fought, it is theoretically rather clear what must be done to achieve small L: filter out very high-energy states because they contribute relatively more to the thermal current than the thermoelectric or Ohmic currents. They could, in theory, be either 1) filtered with additional states that locally accommodate heavy scattering at high energies, or 2) better yet by simple absence of high-energy states.

C. More on Optimal Bandwidth
Investigation of optimal electronic structures for thermoelectrics can be traced back to the seminal work by Mahan and Sofo 1 . They took a purely mathematical approach to formulate zT in terms of the energy integrals presented in the main text, and derived the ideal spectral conductivity for maximization of zT . They determined that a Dirac delta function near the Fermi level, say at E † ≈ E F , is the ideal functional form: where Σ 0 is some pre-factor. D(E) and v 2 (E) are directly determined by the electronic structure, while τ (E) is only indirectly related to it and heavily depends on electron scat-tering mechanisms. However, their approach was purely mathematical in nature with an implicit assumption that κ lat = 0. Eq. S5 must be interpreted with caution when applied to reality.
For Eq. S5 to be satisfied, mathematically, at least one of D(E), v 2 (E), or τ (E) must be δ E,E † . Physically, however, the terms cannot be independently reduced to a delta function.
Firstly, v 2 (E) cannot be δ E,E † while D(E) is not, because v 2 (E) is categorically zero without some band dispersion around E F , while no band dispersion can arise at all if D(E) has no width. Secondly, provided there is some dispersion, τ (E) cannot be a delta function unless electrons are perfectly scattered everywhere but at single E F , which is next to impossible.
Then the only plausible way in which τ (E) can be a delta function is if D(E) is also. These considerations indicate that the only way for Eq. S5 to hold is for D(E) = N v δ E,E † , reflecting one or more (N v ) perfectly localized states, or perfectly flat bands. The factor N v arises from the fact that eDOS of each band must integrate to 1 (or 2 if spin-degenerate) to conserve the number of electrons: This in turn forces τ (E) = τ 0 δ E,E † for some finite τ 0 limited by elastic scattering, but more importantly it forces v 0 → 0 and therefore by default Σ(E) → 0.
The implications of Σ(E) → 0 are as follows. First, the conductivity is immediately 0, as has also been pointed out by a previous study 2 , Second, the Seebeck coefficient is expressible as the following limit as Σ(E) → 0 pointby-point: Generally speaking, this is a non-trivial limit to evaluate because Σ(E) is a function, not a scalar. However, with the knowledge that Σ(E) is widthless (due to D(E)), and so it would approach zero at a single point, we can reformulate the limit as lim Σ(E)→0 where n ≥ 1 is an integer, and evaluate Eq. S8 as and this behavior is graphically verified in Fig. 2a and Fig. 3c in the main text. In Fig. 2a, where E † tends to E F at the band minimum as the band narrows, α tends to 0. In Fig. 3c, where E F is away from the band minimum that E † tends to, α tends to some finite value corresponding to (E F − E † ).
Third, by the same token as above, the Lorenz number can be shown to tend to 0 for a widthless band: and this behavior is graphically suggested in Fig. S4e. By Eqs. S10 and S11, given some E † , a perfectly localized band of widthless D(E) would lead to divergence in "electronic-part zT ," or zT without κ lat : This behavior, consistent with the conclusions of the Mahan-Sofo theory, is graphically suggested in Fig. S4f. However, because of finite Seebeck and vanishing conductivity, the PF would vanish, and compounded by κ lat > 0 in real materials, zT would vanish to 0 alike.
Ergo, even if a set of perfectly localized states could exist in real materials, it is not the physically ideal structure for zT or the PF in real materials. The fundamental barrier is, again, that the components of Σ(E) cannot be independently widthless. The value of z e T as a metric for thermoelectric performance improves only as v(E) becomes finite and large and as κ lat is kept minimal.
In supplement to the optimum bandwidth of a parabolic band, we also consider an isotropic quartic band, E = c(k 2 x + k 2 y + k 2 z ) 2 . The quartic dispersion coefficient c is selected such that it gives the quartic band the same energy as the parabolic band at the inflection point. Fig. S10 shows the results. Between parabolic and quartic dispersions, we see that the latter performs better by about 20%. Since τ ∝ D −1 (E) under DPS, it is again the average group velocities distribution that is responsible for the better quartic performance in the transport direction: v 2 x (E) = 16 √ c 3 E 1.5 . This obviously grows faster with E than that for the parabolic v 2 x (E) = 2 3m E, thereby weighting ζ relatively more than σ.
We also perform the optimization on an isotropic quartic band of comparable dispersion.
We find that the optimium quartic zT exceeds the optimum parabolic zT by around 20%.
Quartic dispersion would perform better under other scattering mechanisms than DPS as well, since the said benefit from v 2 x (E) is still just as effective, and the D(E) ∝ E − 1 4 profile would generally supply higher lifetimes to higher energies than to lower energies.