We use first-principles calculations to reveal the effects of divalent Pb, Ca, and Sn doping of Bi_{2}Te_{3} on the band structure and transport properties, including the Seebeck coefficient, α, and the reduced power factor, α^{2}σ/τ, where σ is the electrical conductivity and τ is the effective relaxation time. Pb and Ca additions exhibit up to 60%–75% higher peak α^{2}σ/τ than that of intrinsic Bi_{2}Te_{3} with Bi antisite defects. Pb occupancy and Ca occupancy of Bi sites increase σ/τ by activating high-degeneracy low-effective-mass bands near the valence band edge, unlike Bi antisite occupancy of Te sites that eliminates near-edge valence states in intrinsic Bi_{2}Te_{3}. Neither Pb doping nor subatomic-percent Ca doping increases α significantly, due to band averaging. Higher Ca levels increase α and diminish σ, due to the emergence of a corrugated band structure underpinned by high-effective-mass bands, attributable to Ca–Te bond ionicity. Sn doping results in a distortion of the bands with a higher density of states that may be characterized as a resonant state but decreases α^{2}σ by up to 30% due to increases in the charge carrier effective mass and decreases in both spin–orbit coupling and valence band quasidegeneracy. These results, and thermal conductivity calculations for nanostructured Bi_{2}Te_{3}, suggest that Pb or Ca doping can enhance the thermoelectric figure of merit ZT to values up to ZT ∼ 1.7, based on an experimentally determined τ. Our findings suggest that divalent doping can be attractive for realizing large ZT enhancements in pnictogen chalcogenides.

## I. INTRODUCTION

Realizing thermoelectric devices with >20% device energy conversion efficiencies^{1–3} requires materials with simultaneously high Seebeck coefficient α, high electrical conductivity σ, and low thermal conductivity κ. This criterion is met when the material has a thermoelectric figure of merit ZT > 3, where $ZT=\alpha 2\sigma T\kappa $, and T is the absolute temperature. Much work has focused on lowering κ by suppressing phonon transport through crystal design, nanostructuring, and alloying while preserving a high bulklike power factor α^{2}σ.^{1,4,5} However, it is widely recognized that further ZT increases require realizing α^{2}σ values beyond limits imposed by natural bulk properties.

Our recent theoretical work has shown that manipulating the electronic structure is an attractive means to obtain α^{2}σ values multifold higher than those of intrinsic single crystals.^{6} For example, we showed that manipulating lattice strain in Bi_{2–x}Sb_{x}Te_{3–y}Se_{y} alloys can increase α^{2}σ by ∼250% due to favorable distortions of the band structure and resulting density of states (DOS). This is corroborated by experiments showing up to eightfold greater power factors compared to bulk alloys.^{7,8} These exciting results underscore the importance of correlations between the electronic band topology and power factor, motivating us to explore subatomic-percent doping as a low-pressure means to realize such beneficial band topologies.

Here, we reveal the effects of three divalent p-type dopants, namely, Pb, Ca, and Sn, on the electronic structure and power factor of Bi_{2}Te_{3}. We chose these dopants because subatomic-percent additions of Pb or Sn have been shown^{9,10} to result in 10%–50% higher power factor α^{2}σ, but the underlying mechanisms are not yet clear. Ca doping has been shown to allow tuning of the Fermi energy E_{F} to the Dirac point,^{11} but the consequent effects on band structure, α, and σ are not known. Accordingly, for each dopant addition, we computed the Bi_{2}Te_{3} electronic band structures using first-principles density functional theory (DFT) calculations. We extracted power factors using Boltzmann transport equations (BTE). Baseline calculations were carried out for intrinsic Bi_{2}Te_{3} to separate the effects of antisite point defects.

Our results show that subatomic-percent Pb doping leads to a 75% α^{2}σ/τ enhancement because Pb occupancy of Bi sites ($PbBi\u2032$) preserves the high-mobility bands that are otherwise suppressed in intrinsic Bi_{2}Te_{3} due to Bi occupancy of Te sites (i.e., $BiTe6c\u2032$ antisite defects). Ca occupancy of Bi sites ($CaBi\u2032$) creates up to 60% higher α^{2}σ/τ by forming a corrugated band structure from a mix of high and low-effective-mass bands. In contrast, Sn occupancy of Bi sites ($SnBi\u2032$) results in a 30% lower α^{2}σ/τ due to decreased spin–orbit coupling and increased charge carrier effective mass arising from resonant state formation. These insights along with thermal conductivity calculations for nanostructured Bi_{2}Te_{3} suggest that ZT ≥ 1.7 could be attained by Pb or Ca doping of nanostructured Bi_{2}Te_{3}.

## II. FIRST-PRINCIPLES COMPUTATIONAL METHODS

Our DFT calculations used the Vienna *ab initio* simulation package (VASP v5.3.5)^{12,13} with projector-augmented wave (PAW) pseudopotentials and the modified Perdew–Burke–Ernzerhof generalized gradient approximation (PBEsol) exchange-correlation functional.^{14} The plane wave energy cutoff was set to 300 eV. Except ionic relaxation, all our calculations included spin–orbit coupling. Ionic relaxations ceased when interatomic forces fell to < 0.1 eV/nm, and electronic relaxations completed when the difference in ground state energy between successive steps was <10^{−5} eV.

Ionic and electronic relaxations were performed using Γ-centered, 12 × 12 × 12, 4 × 4 × 4, and 2 × 2 × 2 k-meshes for the rhombohedral Bi_{2}Te_{3} unit cell [Fig. 1(a)], 135-atom 3 × 3 × 3 supercell [Fig. 1(b)], and 320-atom 4 × 4 × 4 supercell, respectively. A single defect in the 320-atom supercell allowed the simulation of ∼0.3 at. % doping. Doping levels of ∼0.75 and ∼2.2 at. % were achieved in the 135-atom supercell by introducing 1 and 3 dopant atoms, respectively. For the 2.2 at. % concentration, we varied the separation distances of the 3 defects (see Sec. IV E), and the configuration with the lowest ground state energy was chosen for further analysis.

Any nondefective supercell was otherwise identical to a primitive unit cell (see band structure comparison in the supplemental material). We sampled the electronic structure with 100 k-points between each high-symmetry point along the L-Z-a-Γ-Z-F path [Fig. 1(c)] in the first Brillouin zone (IBZ) to obtain band structure cross sections from ellipsoidal energy surfaces near the valence band edge (VBE). The reciprocal lattice coordinates of these high-symmetry k-points^{6} scale linearly with supercell size. For clarity in comparing the effects of different dopants, we plotted only the band closest to the VBE, since bands within ∼3k_{B}T of E_{F} dominate transport.^{15} We determined the electron localization function (ELF) and density of states (DOS) using the linear tetrahedron method to sample the converged charge density over the irreducible Brillouin zone (IBZ). We used symmetrized 36 × 36 × 36, 8 × 8 × 8, and 6 × 6 × 6 Γ-centered k-meshes to sample the unit cell, 3 × 3 × 3 supercells, and 4 × 4 × 4 supercells, respectively.

We used BoltzTraP,^{16} the version available in 2016, to extract α, σ/τ, and κ_{e}—the electronic contribution to thermal conductivity, at room temperature (300 K) along the high-ZT $[101\xaf0]$ direction, from the band structure data used for DOS calculations, but interpolated onto a k-mesh 5 times as dense as the mesh used in our DFT calculations. We applied a 0.055 eV scissors bandgap shift to nondefective Bi_{2}Te_{3} to match the experimental value E_{G} = 0.15 eV.^{17} Scissors shifts were applied to each of our doped systems to reach the same E_{G}, to enable us to separate the effects of band topology from bandgap. Hole concentrations were obtained by dividing band occupancy at the calculated Fermi level E_{F} by the respective supercell volume.

We extracted α and στ^{−1} from the band structures of Bi_{2}Te_{3} with $BiTe6c\u2032$ defects (see Sec. IV A). We used a charge scattering relaxation time τ = 3.4 × 10^{−14} s to match our calculated α^{2}σ values with intrinsic, single-crystal experimental data.^{9,18–20} This τ value was subsequently used to determine σ values for all systems. For the hole concentration range 10^{19} ≤ h ≤ 10^{20} cm^{−3}, our α and σ values determined from each nondefective supercell were within 2% of those obtained from a unit cell. We report thermoelectric transport properties over a range of carrier concentrations for each dopant. Our aim is to mainly capture the overall trends of how doping and antisite defects separately impact the electronic band structure and hence the thermoelectric properties. Our results could serve as a basis for interpreting and understanding experimental results that can be impacted by a variety of factors.^{9,10,18} We note that the effects of codoping, chemical stability, or formation energy, or doping efficiency are beyond the scope of this work.

## III. BISMUTH TELLURIDE PROPERTIES

### A. Crystal structure

Bi_{2}Te_{3} has a rhombohedral structure [Fig. 1(a)] specified by the $R3\xafm$ space group with Te_{6c}-Bi_{6c}-Te_{3a}-Bi_{6c}-Te_{6c} quintuple layer stacking, wherein Bi_{6c} sites are given by fractional coordinates (±0.400, ±0.400, ±0.400), Te_{3a} sites by (0,0,0), and Te_{6c} sites by (±0.209, ±0.209, ±0.209). Our calculated lattice parameters $aR$ = 1.042 nm and $\theta R$ = 24.322° are within 1% of experimentally determined values.^{21} Since Pb,^{9,22} Ca,^{11,23} and Sn^{10,23} are known to occupy pnictogen sites as acceptors in V_{2}VI_{3} materials, we placed the dopants at Bi sites in the supercell.

### B. Orbital interaction

In Bi_{2}Te_{3}, experimental charge carrier effective mass ranges between 0.06m_{0} ≤ m* ≤ 0.08m_{0} for electrons and 0.08m_{0} ≤ m* ≤ 0.11m_{0} for holes,^{24,25} where m_{0} is the free electron mass and $m\u2217=(mx\u2217my\u2217mz\u2217)1/3$ is the geometric mean of the Cartesian surface curvatures of energy surfaces in the IBZ. Bohr ionization energies $EI=e4m\u22172(4\pi \u03f5\u03f50\u210f)2$ are, therefore, ≤ 0.44 meV for donors and ≤ 0.6 meV for acceptors,^{26} both of which are <0.5% of the bandgap 0.13 ≤ E_{g} ≤ 0.15 eV.^{17} These shallow dopant levels are easily ionizable and strongly interact with the host bands, consistent with experimentally observed low-temperature charge carrier saturation in intrinsic p-Bi_{2}Te_{3}.^{19} The Bohr radii $rB=4\pi \u03f5\u03f50\u210f2m\u2217e2$ ^{26} are >33 nm for donors and >24 nm for acceptors, indicating delocalized effects for defect concentrations ≤ 10^{−4} at. %. Despite a factor-of-two uncertainty,^{26} our estimates of E_{I} and r_{B} indicate that Bi_{2}Te_{3} should have nearly completely ionized dopants at all concentrations investigated here.

## IV. RESULTS AND DISCUSSION

### A. Intrinsic defects

Our DFT calculations reveal that Bi occupancy of Te sites (i.e., $BiTe6c\u2032$ antisite defect) alters the electronic structure and, therefore, the transport properties in Bi_{2}Te_{3}. Our calculated electron localization function (ELF) map (Fig. 2) shows an increased electron localization near the center of Bi–$BiTe6c\u2032$ bonds, indicating increased covalency compared to ionic Bi–Te_{6c} bonds in nondefective Bi_{2}Te_{3},^{27,28} where the ELF is high near the atoms. The increased covalency due to Bi–$BiTe6c\u2032$ bonds should promote high hole mobility $\mu h$ because of increased overlaps in *ppσ* bonds.^{29–31}

The relatively low value of our calculated 200 ≤ $\mu h$ ≤ 300 cm^{2 }V^{−1 }s^{−1} suggests the primacy of a factor other than bonding. Bi antisite occupancy of Te sites decreases the near-edge valence state population built mostly from Te_{6c} *p*-orbitals.^{15} Increasing $[BiTe6c\u2032]$, therefore, tends to retract valence band states farther from the edge in nondefective Bi_{2}Te_{3} (Fig. 3). For example, $[BiTe6c\u2032]$ = 2.2 at. % causes electronic bands along Z-a, a-Γ, Z-F, and $Z\u2212F\u030c$ to retract by 30–70 meV (1.2k_{B}T to 2.7k_{B}T at 300 K), leading to a 50%–75% lower density of states (DOS) near the VBE compared to that in nondefective Bi_{2}Te_{3} (Fig. 4).

Increasing $[BiTe6c\u2032]$ from 0.3 to 2.2 at. % decreases σ by 30% [Fig. 5(a)] due to an equivalent decrease in $\mu h$ [Fig. 5(b)] and decreases α by 20% [Fig. 5(c)]. This result can be understood in terms of the doping-induced movement of the electronic bands near E_{F}. Recalling that individual band degeneracy g represents the number of equivalent energy surfaces in the IBZ near E_{F}^{32} and recognizing multiple inequivalent bands with 2 ≤ g ≤ 6 near the VBE of Bi_{2}Te_{3}, we note that movement of the bands alters the quasidegeneracy, i.e., the total number of bands within ∼3k_{B}T of the VBE. In particular, doping-induced retraction of bands decreases quasidegeneracy, which in turn decreases both α and σ because $\alpha =2kB23e\u210f2mDOS\u2217T\pi 3n23EF$,^{33} where $mDOS\u2217\u221d(g2mx\u2217my\u2217mz\u2217)1/3$,^{6} and $\sigma \u2245\u2211i=1g\sigma i,$^{34} where σ_{i} is the contribution from each band. We obtain the highest α ^{2}σ ∼ 3.9 mW m^{−1 }K^{−2} at h = 3 × 10^{19} cm^{−3} with $[BiTe6c\u2032]$ = 0.3 at. % [Fig. 5(d)]. Increasing [$BiTe6c\u2032]$ to 2.2 at. % reduces α^{2}σ to 1.7 mW m^{−1 }K^{−2}.

### B. Pb doping

The band structure and bonding of Bi_{2}Te_{3} with substitutional $PbBi\u2032$ doping strongly resembles that of nondefective Bi_{2}Te_{3}, and correlates with up to 25% increase in σ over intrinsic Bi_{2}Te_{3} with $BiTe6c\u2032$. In particular, we find that Pb-doped and nondefective Bi_{2}Te_{3} exhibit similar band convergence and quasidegeneracy near the VBE (Fig. 6). For example, for $[PbBi\u2032]$ = 0.75 at. %, the high-degeneracy (g = 6) low-effective-mass Z-F band resides 0.65k_{B}T from the VBE, which is the same as in nondefective Bi_{2}Te_{3}. In contrast, this same low-m* Z-F band is 0.8k_{B}T farther from the VBE for intrinsic Bi_{2}Te_{3} with $[BiTe6c\u2032]$ = 0.75 at. %. This similarity can be understood in terms of the $BiTe6c\u2032$ and $PbBi\u2032$ defect chemistries.

Contrary to the effect of $BiTe6c\u2032$ as discussed above, $PbBi\u2032$ doping does not suppress Te *p*-orbital states. The similar atomic masses of Pb and Bi ensure that $PbBi\u2032$ doping does not significantly alter the spin–orbit coupling effect observed in nondefective Bi_{2}Te_{3}.^{32} However, the higher electronegativity difference between Pb and Te (0.3) compared to Bi and Te (0.2)^{35} suggests that Pb doping increases bond polarity. Indeed, we observe increased electron localization on both Pb and Te atoms (Fig. 7), compared with Bi and Te atoms in nondefective Bi_{2}Te_{3}. This Pb-doping-induced ionicity increase thereby decreases electron wave function overlap, which decreases intraband energy level splitting and thus, $\Delta Ev$. The consequent increase in m* for $[PbBi\u2032]$ = 2.2 at. % manifests as a 0.4k_{B}T retraction of the low-$m\u2217$ Z-F band, and a negligible change in the density of states near E_{F} (Fig. 8).

Preservation of the low-m* Z-F band in $PbBi\u2032$-doped Bi_{2}Te_{3} increases α^{2}σ by enhancing $\mu h$ relative to intrinsic Bi_{2}Te_{3} with equivalent $BiTe6c\u2032$ doping. For example, at $[PbBi\u2032]$ = 0.3^{ }at.^{ }%, this band contributes to a 25% higher σ ∼ 1.9 × 10^{5} Ω^{−1} m^{−1} [Fig. 9(a)] due to an equivalent increase in mobility to $\mu h$ ∼ 400 cm^{2 }V^{−1 }s^{−1} [Fig. 9(b)], at h ∼ 3 × 10^{19} cm^{−3}. Since band contributions are averaged as $\alpha =\u2211i=1g\alpha i\sigma i\u2211i=1g\sigma i$,^{6} α only increases by ∼5% [Fig. 9(c)], suggesting that the contribution of the high-$m\u2217$ bands to α dominates over that of the low-$m\u2217$ Z-F band. Therefore, $PbBi\u2032$ doping of 0.3 at. % provides a ∼30% higher α^{2}σ ∼ 5.1 mW m^{−1 }K^{−2} [Fig. 9(d)]. Increasing Pb doping to $[PbBi\u2032]$ = 0.75 at. % does not further alter α^{2}σ, in contrast to a 25% α^{2}σ decrease seen upon a similar increase in $[BiTe6c\u2032]$ in intrinsic Bi_{2}Te_{3}.

At the highest studied doping level of [$PbBi\u2032$] = 2.2 at. %, the charge carrier mobility decreases by 30% to $\mu h$ ∼ 275 cm^{2 }V^{−1 }s^{−1} [Fig. 9(b)] due to Z-F band retraction from the VBE, while α decreases by <5% [Fig. 9(c)]. Thus, the power factor decreases to α^{2}σ ∼ 3.3 mW m^{−1 }K^{−2} [Fig. 9(d)]. Our results agree well with experimentally measured α^{2}σ ∼ 4.5 mW m^{−1 }K^{−2} in Bi_{2}Te_{3} with $[PbBi\u2032]$ = 0.03^{ }at.^{ }%, which is 50% higher than that measured in Bi_{2}Te_{3} with $BiTe6c\u2032$.^{9} No power factor enhancements over intrinsic Bi_{2}Te_{3} with $BiTe6c\u2032$ defects were experimentally reported for $[PbBi\u2032]$ > 0.2^{ }at.^{ }%.^{9}

### C. Ca doping

We find that substitutional $CaBi\u2032$ doping in Bi_{2}Te_{3} activates new bands, whose topology is sensitive to the Ca doping concentration. Controlling the $CaBi\u2032$ doping level leads to up to 20% higher σ or 30% higher α than in intrinsic Bi_{2}Te_{3} with $BiTe6c\u2032$. Subatomic-percent $CaBi\u2032$ doping shifts multiple secondary low-m* bands toward the VBE. For example, for $[CaBi\u2032]$ = 0.75^{ }at.^{ }%, the low-m* Z-F energy ellipsoid transforms such that an alternate cross section labeled as the $Z\u2212F\u030c$ band (Fig. 10), aligns with the a-Γ band at the VBE at ∼0.8k_{B}T higher than in nondefective Bi_{2}Te_{3}. This ellipsoid transformation also manifests as a movement of the left portion of the Z-F band 1.2k_{B}T closer to the VBE. In addition, two low-m* bands along L-Z and Z-a shift to within 2k_{B}T of the VBE. Collectively, this results in an up to threefold increase in DOS within 3k_{B}T from the VBE (Fig. 11) compared to Bi_{2}Te_{3} with equivalent $[BiTe6c\u2032]$.

Such a convergence of low-m* bands at $[CaBi\u2032]$ = 0.3^{ }at.^{ }% yields a 20% higher σ ∼ 1.8 × 10^{5 }Ω^{−1 }m^{−1} at h ∼ 3 × 10^{19} cm^{−3} [Fig. 12(a)] compared to that with equivalent $[BiTe6c\u2032]$, due to a high $\mu h$ ∼ 375^{ }cm^{2 }V^{−1 }s^{−1} [Fig. 12(b)]. The 10% increase in α [Fig. 12(c)] is twofold greater than that observed for equivalent $PbBi\u2032$ doping. The net effect is a ∼33% increase in the power factor to α^{2}σ ∼ 5.3 mW m^{−1 }K^{−2} [Fig. 12(d)]. At $[CaBi\u2032]$ = 0.75^{ }at.^{ }%, α^{2}σ ∼ 4.7 mW^{ }m^{−1 }K^{−2}, which is 8% lower than that obtained with equivalent $[PbBi\u2032]$, but still 60% higher than in Bi_{2}Te_{3} with equivalent $[BiTe6c\u2032].$

Increasing Ca doping to $[CaBi\u2032]$ = 2.2^{ }at.^{ }% creates bands with lower curvatures (Fig. 10) and thus higher effective-mass than at $[CaBi\u2032]$ = 0.75^{ }at.^{ }% near the VBE. In particular, the highest-effective-mass Γ-Z band advances to within 1.6k_{B}T of the band edge of intrinsic Bi_{2}Te_{3} (Fig. 10). The valence band maxima (VBM) along L-Z and a-Γ retract ≤ 0.8k_{B}T along with both the Z-F and $Z\u2212F\u030c$ bands, which transform into flatter, higher-effective-mass bands. The Z-a band flattens and advances to the VBE. The collective effect is a more than tenfold increase in DOS near the VBE (Fig. 11), created by a fivefold increase in $mDOS\u2217$, if we assume that DOS $\u221dmDOS\u22173/2$.^{15} The higher m* (lower $\Delta Ev$) is consistent with the higher ionicity^{35} of $CaBi\u2032$–Te_{6c} and $CaBi\u2032$–Te_{3a} bonds (Fig. 13) than the Bi–Te counterparts in nondefective Bi_{2}Te_{3}. The $CaBi\u2032$ defects form equal-length bonds of high ionic character with both Te_{6c} and Te_{3a} involving electron transfer from Ca. In nondefective Bi_{2}Te_{3}, however, Bi–Te_{6c} bonds formed by electron transfer from Te are 5% shorter than the covalent Bi–Te_{3a} bonds.^{28} Since Ca atomic mass is only ∼20% of that of Bi, $CaBi\u2032$ doping decreases spin–orbit coupling, which usually inhibits convergence of bands involving 5*p–*6*p* orbitals.^{32} Our observation of increased band convergence at all $[CaBi\u2032]$ suggests that the interaction of 4*s–*5*p* orbitals in Ca-doped Bi_{2}Te_{3} is fundamentally different than the 5*p*–6*p* orbital interactions for $BiTe6c\u2032$, $PbBi\u2032$, and $SnBi\u2032$ doping.

The increased m* with $[CaBi\u2032]$ = 2.2^{ }at.^{ }% manifests as a 20% higher α ∼ 200 *μ*V K^{−1} at h ∼ 3 × 10^{19} cm^{−3}, which is 30% higher than Bi_{2}Te_{3} with optimum $[BiTe6c\u2032]$ = 0.3^{ }at.^{ }%, and 65% higher than with equivalent $[BiTe6c\u2032]$ = 2.2^{ }at.^{ }%. The σ decreases by 30% due to a μ_{h} decrease to 160 cm^{2 }V^{−1 }s^{−1} (Fig. 12), leading to a net α^{2}σ ∼ 3.2 mW^{ }m^{−1 }K^{−2}. This power factor is only 40% of the maximum achieved with subatomic-percent $CaBi\u2032$ doping, indicating that enhancements in α due to high-effective-mass bands cannot fully counteract the concomitant diminution of μ_{h}. These results suggest that highly degenerate, low-effective-mass bands, generated by covalent bonds, are perhaps more important for realizing high α^{2}σ.

### D. Sn doping

Experimentally, subatomic-percent Sn doping yields 40%–75% higher α than undoped Bi_{2}Te_{3} at room temperature because of a 25%–100% increase in m*.^{10,23} However, our first-principles calculations indicate that $SnBi\u2032$ doping of Bi_{2}Te_{3} produces lower α^{2}σ than equivalent $BiTe6c\u2032$ doping of intrinsic Bi_{2}Te_{3}, despite the formation of an Sn-induced resonant state, i.e., a DOS peak within the valence band. Although $SnBi\u2032$ doping does not suppress Te *p*-orbital states, we observe that all secondary bands retract away from the VBE, reminiscent of the decreased quasidegeneracy in intrinsic Bi_{2}Te_{3} with $BiTe6c\u2032$ doping (Fig. 14). Decreased spin–orbit coupling resulting from the ∼43% lower Sn atomic mass than Bi^{32} is the likely reason for retraction of the bands. The $SnBi\u2032$–Te_{6c} bond has a higher covalent character than the Bi–Te_{6c} bond, as seen from the slightly higher electron localization half-way between the Sn–Te bond (Fig. 15). This result is consistent with the 30% lower electronegativity difference of Sn–Te (0.14) than Bi–Te (0.2).^{35}

The possibility of charge carrier mobility enhancement due to increased covalency, however, is counteracted by increases in band effective mass. All secondary VBMs in $SnBi\u2032$-doped Bi_{2}Te_{3} have lower curvatures (higher m*) than those in intrinsic Bi_{2}Te_{3} with equivalent $BiTe6c\u2032$ levels (Fig. 14). For example, at $[SnBi\u2032]$ = 0.75%, Z-a, Z-F, and $Z\u2212F\u030c$ band edges have ∼15% higher m* than in intrinsic Bi_{2}Te_{3}. At $[SnBi\u2032]$ = 2.2^{ }at.^{ }%, these three bands have ∼40% higher m*. These near-edge band topology changes are consistent with the 50%–200% higher DOS than that of nondefective Bi_{2}Te_{3}, attributable to distortion of the band structure, termed a resonant state deeper in the valence band [Fig. 16(a)]. Assuming DOS $\u221dmDOS\u22173/2$,^{15} the resonant state corresponds to 30%–110% higher $mDOS\u2217$ than intrinsic Bi_{2}Te_{3} in agreement with an experiment.^{10}

Decreasing $SnBi\u2032$ doping shifts the resonant state toward the VBE, e.g., from 200 meV below at $[SnBi\u2032]$ = 2.2 at. %, to 95 meV below the VBM at $[SnBi\u2032]$ = 0.3 at. % [Fig. 16(b)]. The reported experimental value is 15 meV at [Sn] ∼ 0.1 at. %.^{36} Our calculations reveal two α^{2}σ peaks: a primary peak around h ∼ 10^{19}–10^{20} cm^{−3} is followed by a smaller secondary peak at h ∼ 10^{20}–10^{21} cm^{−3} range [Figs. 16(c) and 17(d)]. These results mirror experimental results^{10} showing a primary α^{2}σ peak at $[SnBi\u2032]$ = 0.05 at. % (h = 1 × 10^{19} cm^{−3}), followed by a secondary peak at $[SnBi\u2032]$ = 0.3 at. % (h = 5 × 10^{19} cm^{−3}). Additionally, our calculations show that decreasing the Sn doping level increases the primary peak magnitude and shifts the secondary peak to lower hole concentrations. Thus, if we were able to model dopant concentrations more similar to experiment (0.05–0.1 at. %, requiring ∼$7\xd77\xd77$ to ∼$6\xd76\xd76$ supercell size), we expect that the secondary α^{2}σ peak would shift into the range of h ∼ 1–5 × 10^{19} cm^{−3}, and increase the primary peak magnitude, in agreement with the experimentally observed α^{2}σ enhancement.

Sn-doping-induced power factor behavior can be understood in terms of the nature and relative magnitude of the effects of Sn doping on σ and α. In particular, $SnBi\u2032$ doping increases m* due to resonant state formation, but decreases α^{2}σ due to a vastly greater decrease in σ. For example, at $[SnBi\u2032]$ = 0.3 at. %, σ ∼ 9 × 10^{4} Ω^{−1} m^{−1} [Fig. 17(a)] is 35% lower than Bi_{2}Te_{3} with equivalent $BiTe6c\u2032$ concentration, due to a decreased $\mu h$ ∼ 185 cm^{2 }V^{−1 }s^{−1} [Fig. 17(b)] at h ∼ 3 × 10^{19} cm^{−3}. A 10% increase in α ∼ 175 μV K^{−1} [Fig. 17(c)] slightly offsets the decreased mobility to produce a peak α^{2}σ ∼ 2.75 mW m^{−1 }K^{−2} [Fig. 17(d)] which is ∼30% lower than the value obtained from intrinsic Bi_{2}Te_{3} with equivalent $[BiTe6c\u2032$]. Increasing Sn doping to $[SnBi\u2032]$ = 2.2 at. % decreases α^{2}σ to ∼1.4 mW m^{−1} ^{−2}, which is ∼20% lower than the value obtained from intrinsic Bi_{2}Te_{3} with $[BiTe6c\u2032$] = 2.2 at. %.

### E. Defect configuration

The α^{2}σ of Bi_{2}Te_{3} is sensitive to both defect type and configuration. $BiTe6c\u2032$ defects show the strongest α^{2}σ disparity between clustered and distributed configurations. At $[BiTe6c\u2032]$ = 2.2 at. %, distributed defects have a separation distance of 0.75 ≤ d_{sep} ≤ 1.5 nm. However, for the same defect concentration, if we force $BiTe6c\u2032$ defects to cluster with d_{sep }∼ 0.4 nm, we get a 0.4 meV/atom lower total energy. This result indicates that clustered antisite defects, which yield a lower α^{2}σ, are energetically favored. For example, for d_{sep }∼ 0.4 nm, α^{2}σ ∼ 1.7 mW m^{−1 }K^{−2} is more than twofold lower than α^{2}σ ∼ 4 mW m^{−1 }K^{−2} obtained at higher d_{sep}.

If we decrease d_{sep} by increasing $[BiTe6c\u2032]$, we increase the Bi–$BiTe6c\u2032$ repulsion as well as Te_{6c}–$BiTe6c\u2032$ attraction for d_{sep} ≥ 0.75 nm (see Fig. 18). In contrast, if we force d_{sep }∼ 0.4 nm for constant $[BiTe6c\u2032]$ = 2.2 at. %, both bond lengths decrease$.$ Thus, $BiTe6c\u2032$–Bi and/or $BiTe6c\u2032$–$BiTe6c\u2032$ bonds are responsible for degrading α^{2}σ due to band retraction from the VBE, as shown in Fig. 19. Such differences in defect clustering (inset of Fig. 19) may explain experimental outliers, e.g., α^{2}σ > 4 mW m^{−1 }K^{−2} at h > 10^{20} cm^{−3}.^{20}

Forcing $PbBi\u2032$ defects to cluster with d_{sep} ∼ 0.4 nm at $[PbBi\u2032]$ = 2.2 at. % yields a 0.1 meV/atom lower ground state energy (Fig. 20) than sparsely distributed $PbBi\u2032$ defects with 0.75 ≤ d_{sep} ≤ 1.5 nm. This result suggests a weaker preference for $PbBi\u2032$ clustering than the above-discussed $BiTe6c\u2032$ clustering. The distributed $PbBi\u2032$ configuration yields only 5% lower α^{2}σ due to the similarity of the valence 6*p* orbitals of Pb and Bi. Given the similar energetics of defect clustering, and the opposing effects of $PbBi\u2032$and $BiTe6c\u2032$ on band structure and α^{2}σ, it is conceivable that $PbBi\u2032$ doping suppresses $BiTe6c\u2032$ formation by altering the defect energetics.

Contrary to Pb doping, $CaBi\u2032$ defects with separations in the 0.75 ≤ d_{sep} ≤ 1.5 nm range exhibit a 0.3 meV/atom lower total energy (Fig. 20) than when clustered such that d_{sep} ∼ 0.4 nm, indicating that distributed $CaBi\u2032$ defects are energetically preferred. When $CaBi\u2032$ defects are forced to cluster, α^{2}σ decreases by 20%. Similar to Ca doping, distributed $SnBi\u2032$ defects yield a ∼0.2 meV/atom lower total energy (Fig. 20) at 2.2 at. % concentration, suggesting that Sn does not prefer to cluster in Bi_{2}Te_{3}. When $SnBi\u2032$ defects are forced to cluster, α^{2}σ decreases by 25%.

### F. Maximum ZT in nanostructured Bi_{2}Te_{3}

Results of our first-principles calculations described above clearly demonstrate $PbBi\u2032$ and $CaBi\u2032$ doping of Bi_{2}Te_{3} can result in ≥30% higher α^{2}σ/τ than that in intrinsic Bi_{2}Te_{3} with $BiTe\u2032$ antisite defects (Fig. 21). In addition, we show that subatomic-percent dopant concentrations yield the highest maximum α^{2}σ, providing up to ≥40% higher α^{2}σ compared to the respective highly doped systems. Since τ increases with decreasing dopant concentration, our use of a constant τ implies that our calculated increases in α^{2}σ at subatomic-percent dopant concentration are a conservative estimate. We show below that we can potentially realize room-temperature ZT ∼ 1.7 in Bi_{2}Te_{3} by combining such doping-induced electronic property enhancements with low thermal conductivity achieved in bulk pellets with an average grain size of d_{G} ∼ 10 nm (Fig. 22). This represents a nearly twofold increase over the highest reported ZT = 0.87 (Refs. 37–39) from high-quality intrinsic single-crystal Bi_{2}Te_{3}, and a 240% increase over single-crystal p-Bi_{2}Te_{3} with a more typical ZT ∼ 0.5,^{9,18,20} which we used to benchmark our DFT-based transport calculations (Sec. IV A). Nanostructuring with 10 nm grains enables a 60% decreased lattice thermal conductivity κ_{L}, accompanied by only a 20% decrease in α^{2}σ, due to a 35% decrease in σ and a 10% increase in α (Fig. 21).

The key question raised by the above results is the extent to which the benefits of dopant induced band structure distortions can be isolated from possible changes in relaxation time. Since BoltzTraP cannot model the scattering, we analytically solved Boltzmann transport equations for α , σ, and κ_{e}, using the nonparabolic two-band Kane model.^{15} We baselined the models^{40} against experimental single-crystal transport data^{19} by implementing experimental carrier effective masses and bandgap (see Sec. III B). Energy-dependent optical, acoustic, and polar-optical phonon–electron scattering were included using standard expressions,^{15} and electron-defect scattering with the Brooks–Herring model.^{41} We have experimentally and theoretically shown that κ_{L} in nanograined Bi_{2}Te_{3} is unaffected by the subatomic-percent doping and is lower than either the single-crystal $[0001]$ or $[101\xaf0]$ direction.^{42} Thus, our κ_{L} modeling methodology^{42} using modified Debye–Callaway equations^{43} started with κ_{L }∼ 0.85 W m^{−1 }K^{−1} for polycrystalline Bi_{2}Te_{3} with randomly oriented grains, previously employed for the large-grain limit.^{42} Grain scattering was subsequently added to both electronic and thermal models as an additional scattering mechanism with mean-free-path equal to grain diameter d_{G}. At d_{G} = 10 nm, this additional scattering lowers κ_{L} more than σ, due to the grain size being higher than the mean-free-path of holes, but lower than that of phonons.^{42} The slight α enhancement is due to the energy filtering effect, i.e., increased scattering of high-energy carriers.^{44,45}

Using our analytical models, we calculated each parameter in $Z=\alpha 2\sigma (\kappa e+\kappa L)$ for 5 nm ≤ d_{G} ≤ 1 *μ*m at h = 3 × 10^{19} cm^{−3} (Fig. 22) and normalized all values to a grain size of d_{G} = 1 *μ*m, beyond which transport values are constant with grain size. Our DFT-based BoltzTraP transport values for α, σ, and κ_{e} (h = 3 × 10^{19} cm^{−3}) as well as the experimental κ_{L} = 0.85 Wm^{−1} K^{−1} were subsequently scaled using these normalized transport trends to calculate nanostructured ZT for each dopant at 0.3 at. % concentration.

## V. CONCLUSIONS

Our first-principles calculations on Bi_{2}Te_{3} reveal band structure alterations created by Pb, Ca, and Sn doping, as well as by intrinsic Bi antisite defects, and their correlation to thermoelectric power factor. Subatomic-percent Pb and Ca doping enhance the power factor by 30%–75% over Bi_{2}Te_{3} with intrinsic antisite defects, due to the retention and creation of highly degenerate, low-effective-mass bands near the valence band edge, which enhance carrier mobility. Concentrated Ca doping fosters high α, due to the ionic-bonding-induced actuation of high-effective-mass bands, despite the loss in spin–orbit coupling strength. Sn doping also increases the carrier effective mass and reduces the spin–orbit coupling but provides no enhancement in power factor compared with intrinsic p-Bi_{2}Te_{3}. Sn resonant states underlying the increased effective mass do, however, show the capability of sustaining α and μ_{h}, thereby inducing a secondary power factor peak at higher hole concentrations.

We show that higher-than-bulk ZT is best achieved in Bi_{2}Te_{3} nanocrystals doped with subatomic-percent quantities of Pb or Ca. Although a high DOS near band edges is considered beneficial for high thermoelectric efficiency, our work illustrates the importance of the underlying band topology. Only the confluence of many, highly degenerate, low-effective-mass bands, rather than high-effective-mass bands or resonant states, can create the electronic landscape necessary to realize higher-than-bulk power factor. With nanostructuring-induced suppression of κ, ZT ≥ 1.7 can be achieved. We suggest exploring additional elements outside of periodic table groups 4–7 to find dopants that have different orbital characters, can occupy different Wyckoff positions, and thereby introduce new bands along different momenta.

## SUPPLEMENTARY MATERIAL

See the supplementary material for reference data for band structure, effective mass, bonding, transport, Fermi levels and band gaps for each dopant concentration. It also includes experimental comparison data for the correlation of hole and doping concentrations..

## ACKNOWLEDGMENTS

For our calculations, we used the IBM Blue Gene/Q supercomputer at the Rensselaer Center for Computational Innovations. We gratefully acknowledge funding from the Department of Energy Basic Energy Sciences program through the S^{3}TEC EFRC, Award No. DE-SC0001299/DE-FG02-09ER4657.