Divalent doping-induced thermoelectric power factor increase in p-type Bi₂Te₃ via electronic structure tuning

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=α2σTκ$⁠, 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 α²σ.^1,4,5 However, it is widely recognized that further ZT increases require realizing α²σ 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 α²σ values multifold higher than those of intrinsic single crystals.⁶ For example, we showed that manipulating lattice strain in Bi_2–xSb_xTe_3–ySe_y alloys can increase α²σ 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₂Te₃. 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 α²σ, 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,¹¹ but the consequent effects on band structure, α, and σ are not known. Accordingly, for each dopant addition, we computed the Bi₂Te₃ 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₂Te₃ to separate the effects of antisite point defects.

Our results show that subatomic-percent Pb doping leads to a 75% α²σ/τ enhancement because Pb occupancy of Bi sites (⁠$PbBi′$⁠) preserves the high-mobility bands that are otherwise suppressed in intrinsic Bi₂Te₃ due to Bi occupancy of Te sites (i.e., $BiTe6c′$ antisite defects). Ca occupancy of Bi sites (⁠$CaBi′$⁠) creates up to 60% higher α²σ/τ by forming a corrugated band structure from a mix of high and low-effective-mass bands. In contrast, Sn occupancy of Bi sites (⁠$SnBi′$⁠) results in a 30% lower α²σ/τ 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₂Te₃ suggest that ZT ≥ 1.7 could be attained by Pb or Ca doping of nanostructured Bi₂Te₃.

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.¹⁴ 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⁻⁵ 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₂Te₃ 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⁶ 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_BT of E_F dominate transport.¹⁵ 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,¹⁶ 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¯0]$ 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₂Te₃ to match the experimental value E_G = 0.15 eV.¹⁷ 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 στ⁻¹ from the band structures of Bi₂Te₃ with $BiTe6c′$ defects (see Sec. IV A). We used a charge scattering relaxation time τ = 3.4 × 10⁻¹⁴ s to match our calculated α²σ 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¹⁹ ≤ h ≤ 10²⁰ cm⁻³, 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.

Bi₂Te₃ has a rhombohedral structure [Fig. 1(a)] specified by the $R3¯m$ 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 $θR$ = 24.322° are within 1% of experimentally determined values.²¹ Since Pb,^9,22 Ca,^11,23 and Sn^10,23 are known to occupy pnictogen sites as acceptors in V₂VI₃ materials, we placed the dopants at Bi sites in the supercell.

In Bi₂Te₃, experimental charge carrier effective mass ranges between 0.06m₀ ≤ m* ≤ 0.08m₀ for electrons and 0.08m₀ ≤ m* ≤ 0.11m₀ for holes,^24,25 where m₀ is the free electron mass and $m∗=(mx∗my∗mz∗)1/3$ is the geometric mean of the Cartesian surface curvatures of energy surfaces in the IBZ. Bohr ionization energies $EI=e4m∗2(4πϵϵ0ℏ)2$ are, therefore, ≤ 0.44 meV for donors and ≤ 0.6 meV for acceptors,²⁶ both of which are <0.5% of the bandgap 0.13 ≤ E_g ≤ 0.15 eV.¹⁷ 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₂Te₃.¹⁹ The Bohr radii $rB=4πϵϵ0ℏ2m∗e2$ ²⁶ are >33 nm for donors and >24 nm for acceptors, indicating delocalized effects for defect concentrations ≤ 10⁻⁴ at. %. Despite a factor-of-two uncertainty,²⁶ our estimates of E_I and r_B indicate that Bi₂Te₃ should have nearly completely ionized dopants at all concentrations investigated here.

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

The relatively low value of our calculated 200 ≤ $μh$ ≤ 300 cm²V⁻¹s⁻¹ 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.¹⁵ Increasing $[BiTe6c′]$⁠, therefore, tends to retract valence band states farther from the edge in nondefective Bi₂Te₃ (Fig. 3). For example, $[BiTe6c′]$ = 2.2 at. % causes electronic bands along Z-a, a-Γ, Z-F, and $Z−F̌$ to retract by 30–70 meV (1.2k_BT to 2.7k_BT at 300 K), leading to a 50%–75% lower density of states (DOS) near the VBE compared to that in nondefective Bi₂Te₃ (Fig. 4).

Increasing $[BiTe6c′]$ from 0.3 to 2.2 at. % decreases σ by 30% [Fig. 5(a)] due to an equivalent decrease in $μ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³² and recognizing multiple inequivalent bands with 2 ≤ g ≤ 6 near the VBE of Bi₂Te₃, we note that movement of the bands alters the quasidegeneracy, i.e., the total number of bands within ∼3k_BT of the VBE. In particular, doping-induced retraction of bands decreases quasidegeneracy, which in turn decreases both α and σ because $α=2kB23eℏ2mDOS∗Tπ3n23EF$⁠,³³ where $mDOS∗∝(g2mx∗my∗mz∗)1/3$⁠,⁶ and $σ≅∑i=1gσi,$³⁴ where σ_i is the contribution from each band. We obtain the highest α ²σ ∼ 3.9 mW m⁻¹K⁻² at h = 3 × 10¹⁹ cm⁻³ with $[BiTe6c′]$ = 0.3 at. % [Fig. 5(d)]. Increasing [$BiTe6c′]$ to 2.2 at. % reduces α²σ to 1.7 mW m⁻¹K⁻².

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

Contrary to the effect of $BiTe6c′$ as discussed above, $PbBi′$ doping does not suppress Te p-orbital states. The similar atomic masses of Pb and Bi ensure that $PbBi′$ doping does not significantly alter the spin–orbit coupling effect observed in nondefective Bi₂Te₃.³² However, the higher electronegativity difference between Pb and Te (0.3) compared to Bi and Te (0.2)³⁵ 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₂Te₃. This Pb-doping-induced ionicity increase thereby decreases electron wave function overlap, which decreases intraband energy level splitting and thus, $ΔEv$⁠. The consequent increase in m* for $[PbBi′]$ = 2.2 at. % manifests as a 0.4k_BT retraction of the low-$m∗$ 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′$-doped Bi₂Te₃ increases α²σ by enhancing $μh$ relative to intrinsic Bi₂Te₃ with equivalent $BiTe6c′$ doping. For example, at $[PbBi′]$ = 0.3at.%, this band contributes to a 25% higher σ ∼ 1.9 × 10⁵ Ω⁻¹ m⁻¹ [Fig. 9(a)] due to an equivalent increase in mobility to $μh$ ∼ 400 cm²V⁻¹s⁻¹ [Fig. 9(b)], at h ∼ 3 × 10¹⁹ cm⁻³. Since band contributions are averaged as $α=∑i=1gαiσi∑i=1gσi$⁠,⁶ α only increases by ∼5% [Fig. 9(c)], suggesting that the contribution of the high-$m∗$ bands to α dominates over that of the low-$m∗$ Z-F band. Therefore, $PbBi′$ doping of 0.3 at. % provides a ∼30% higher α²σ ∼ 5.1 mW m⁻¹K⁻² [Fig. 9(d)]. Increasing Pb doping to $[PbBi′]$ = 0.75 at. % does not further alter α²σ, in contrast to a 25% α²σ decrease seen upon a similar increase in $[BiTe6c′]$ in intrinsic Bi₂Te₃.

At the highest studied doping level of [$PbBi′$] = 2.2 at. %, the charge carrier mobility decreases by 30% to $μh$ ∼ 275 cm²V⁻¹s⁻¹ [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 α²σ ∼ 3.3 mW m⁻¹K⁻² [Fig. 9(d)]. Our results agree well with experimentally measured α²σ ∼ 4.5 mW m⁻¹K⁻² in Bi₂Te₃ with $[PbBi′]$ = 0.03at.%, which is 50% higher than that measured in Bi₂Te₃ with $BiTe6c′$⁠.⁹ No power factor enhancements over intrinsic Bi₂Te₃ with $BiTe6c′$ defects were experimentally reported for $[PbBi′]$ > 0.2at.%.⁹ 

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

Such a convergence of low-m* bands at $[CaBi′]$ = 0.3at.% yields a 20% higher σ ∼ 1.8 × 10⁵Ω⁻¹m⁻¹ at h ∼ 3 × 10¹⁹ cm⁻³ [Fig. 12(a)] compared to that with equivalent $[BiTe6c′]$⁠, due to a high $μh$ ∼ 375cm²V⁻¹s⁻¹ [Fig. 12(b)]. The 10% increase in α [Fig. 12(c)] is twofold greater than that observed for equivalent $PbBi′$ doping. The net effect is a ∼33% increase in the power factor to α²σ ∼ 5.3 mW m⁻¹K⁻² [Fig. 12(d)]. At $[CaBi′]$ = 0.75at.%, α²σ ∼ 4.7 mWm⁻¹K⁻², which is 8% lower than that obtained with equivalent $[PbBi′]$⁠, but still 60% higher than in Bi₂Te₃ with equivalent $[BiTe6c′].$

Increasing Ca doping to $[CaBi′]$ = 2.2at.% creates bands with lower curvatures (Fig. 10) and thus higher effective-mass than at $[CaBi′]$ = 0.75at.% near the VBE. In particular, the highest-effective-mass Γ-Z band advances to within 1.6k_BT of the band edge of intrinsic Bi₂Te₃ (Fig. 10). The valence band maxima (VBM) along L-Z and a-Γ retract ≤ 0.8k_BT along with both the Z-F and $Z−F̌$ 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∗$⁠, if we assume that DOS $∝mDOS∗3/2$⁠.¹⁵ The higher m* (lower $ΔEv$⁠) is consistent with the higher ionicity³⁵ of $CaBi′$–Te_6c and $CaBi′$–Te_3a bonds (Fig. 13) than the Bi–Te counterparts in nondefective Bi₂Te₃. The $CaBi′$ defects form equal-length bonds of high ionic character with both Te_6c and Te_3a involving electron transfer from Ca. In nondefective Bi₂Te₃, however, Bi–Te_6c bonds formed by electron transfer from Te are 5% shorter than the covalent Bi–Te_3a bonds.²⁸ Since Ca atomic mass is only ∼20% of that of Bi, $CaBi′$ doping decreases spin–orbit coupling, which usually inhibits convergence of bands involving 5p–6p orbitals.³² Our observation of increased band convergence at all $[CaBi′]$ suggests that the interaction of 4s–5p orbitals in Ca-doped Bi₂Te₃ is fundamentally different than the 5p–6p orbital interactions for $BiTe6c′$⁠, $PbBi′$⁠, and $SnBi′$ doping.

The increased m* with $[CaBi′]$ = 2.2at.% manifests as a 20% higher α ∼ 200 μV K⁻¹ at h ∼ 3 × 10¹⁹ cm⁻³, which is 30% higher than Bi₂Te₃ with optimum $[BiTe6c′]$ = 0.3at.%, and 65% higher than with equivalent $[BiTe6c′]$ = 2.2at.%. The σ decreases by 30% due to a μ_h decrease to 160 cm²V⁻¹s⁻¹ (Fig. 12), leading to a net α²σ ∼ 3.2 mWm⁻¹K⁻². This power factor is only 40% of the maximum achieved with subatomic-percent $CaBi′$ 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 α²σ.

Experimentally, subatomic-percent Sn doping yields 40%–75% higher α than undoped Bi₂Te₃ at room temperature because of a 25%–100% increase in m*.^10,23 However, our first-principles calculations indicate that $SnBi′$ doping of Bi₂Te₃ produces lower α²σ than equivalent $BiTe6c′$ doping of intrinsic Bi₂Te₃, despite the formation of an Sn-induced resonant state, i.e., a DOS peak within the valence band. Although $SnBi′$ 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₂Te₃ with $BiTe6c′$ doping (Fig. 14). Decreased spin–orbit coupling resulting from the ∼43% lower Sn atomic mass than Bi³² is the likely reason for retraction of the bands. The $SnBi′$–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).³⁵ 

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′$-doped Bi₂Te₃ have lower curvatures (higher m*) than those in intrinsic Bi₂Te₃ with equivalent $BiTe6c′$ levels (Fig. 14). For example, at $[SnBi′]$ = 0.75%, Z-a, Z-F, and $Z−F̌$ band edges have ∼15% higher m* than in intrinsic Bi₂Te₃. At $[SnBi′]$ = 2.2at.%, 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₂Te₃, attributable to distortion of the band structure, termed a resonant state deeper in the valence band [Fig. 16(a)]. Assuming DOS $∝mDOS∗3/2$⁠,¹⁵ the resonant state corresponds to 30%–110% higher $mDOS∗$ than intrinsic Bi₂Te₃ in agreement with an experiment.¹⁰ 

Decreasing $SnBi′$ doping shifts the resonant state toward the VBE, e.g., from 200 meV below at $[SnBi′]$ = 2.2 at. %, to 95 meV below the VBM at $[SnBi′]$ = 0.3 at. % [Fig. 16(b)]. The reported experimental value is 15 meV at [Sn] ∼ 0.1 at. %.³⁶ Our calculations reveal two α²σ peaks: a primary peak around h ∼ 10¹⁹–10²⁰ cm⁻³ is followed by a smaller secondary peak at h ∼ 10²⁰–10²¹ cm⁻³ range [Figs. 16(c) and 17(d)]. These results mirror experimental results¹⁰ showing a primary α²σ peak at $[SnBi′]$ = 0.05 at. % (h = 1 × 10¹⁹ cm⁻³), followed by a secondary peak at $[SnBi′]$ = 0.3 at. % (h = 5 × 10¹⁹ cm⁻³). 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×7×7$ to ∼$6×6×6$ supercell size), we expect that the secondary α²σ peak would shift into the range of h ∼ 1–5 × 10¹⁹ cm⁻³, and increase the primary peak magnitude, in agreement with the experimentally observed α²σ 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′$ doping increases m* due to resonant state formation, but decreases α²σ due to a vastly greater decrease in σ. For example, at $[SnBi′]$ = 0.3 at. %, σ ∼ 9 × 10⁴ Ω⁻¹ m⁻¹ [Fig. 17(a)] is 35% lower than Bi₂Te₃ with equivalent $BiTe6c′$ concentration, due to a decreased $μh$ ∼ 185 cm²V⁻¹s⁻¹ [Fig. 17(b)] at h ∼ 3 × 10¹⁹ cm⁻³. A 10% increase in α ∼ 175 μV K⁻¹ [Fig. 17(c)] slightly offsets the decreased mobility to produce a peak α²σ ∼ 2.75 mW m⁻¹K⁻² [Fig. 17(d)] which is ∼30% lower than the value obtained from intrinsic Bi₂Te₃ with equivalent $[BiTe6c′$]. Increasing Sn doping to $[SnBi′]$ = 2.2 at. % decreases α²σ to ∼1.4 mW m^−1 −2, which is ∼20% lower than the value obtained from intrinsic Bi₂Te₃ with $[BiTe6c′$] = 2.2 at. %.

The α²σ of Bi₂Te₃ is sensitive to both defect type and configuration. $BiTe6c′$ defects show the strongest α²σ disparity between clustered and distributed configurations. At $[BiTe6c′]$ = 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′$ 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 α²σ, are energetically favored. For example, for d_sep∼ 0.4 nm, α²σ ∼ 1.7 mW m⁻¹K⁻² is more than twofold lower than α²σ ∼ 4 mW m⁻¹K⁻² obtained at higher d_sep.

If we decrease d_sep by increasing $[BiTe6c′]$⁠, we increase the Bi–$BiTe6c′$ repulsion as well as Te_6c–$BiTe6c′$ attraction for d_sep ≥ 0.75 nm (see Fig. 18). In contrast, if we force d_sep∼ 0.4 nm for constant $[BiTe6c′]$ = 2.2 at. %, both bond lengths decrease$.$ Thus, $BiTe6c′$–Bi and/or $BiTe6c′$–$BiTe6c′$ bonds are responsible for degrading α²σ 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., α²σ > 4 mW m⁻¹K⁻² at h > 10²⁰ cm⁻³.²⁰ 

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

Contrary to Pb doping, $CaBi′$ 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′$ defects are energetically preferred. When $CaBi′$ defects are forced to cluster, α²σ decreases by 20%. Similar to Ca doping, distributed $SnBi′$ 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₂Te₃. When $SnBi′$ defects are forced to cluster, α²σ decreases by 25%.

Results of our first-principles calculations described above clearly demonstrate $PbBi′$ and $CaBi′$ doping of Bi₂Te₃ can result in ≥30% higher α²σ/τ than that in intrinsic Bi₂Te₃ with $BiTe′$ antisite defects (Fig. 21). In addition, we show that subatomic-percent dopant concentrations yield the highest maximum α²σ, providing up to ≥40% higher α²σ compared to the respective highly doped systems. Since τ increases with decreasing dopant concentration, our use of a constant τ implies that our calculated increases in α²σ at subatomic-percent dopant concentration are a conservative estimate. We show below that we can potentially realize room-temperature ZT ∼ 1.7 in Bi₂Te₃ 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₂Te₃, and a 240% increase over single-crystal p-Bi₂Te₃ 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 α²σ, 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.¹⁵ We baselined the models⁴⁰ against experimental single-crystal transport data¹⁹ 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,¹⁵ and electron-defect scattering with the Brooks–Herring model.⁴¹ We have experimentally and theoretically shown that κ_L in nanograined Bi₂Te₃ is unaffected by the subatomic-percent doping and is lower than either the single-crystal $[0001]$ or $[101¯0]$ direction.⁴² Thus, our κ_L modeling methodology⁴² using modified Debye–Callaway equations⁴³ started with κ_L∼ 0.85 W m⁻¹K⁻¹ for polycrystalline Bi₂Te₃ with randomly oriented grains, previously employed for the large-grain limit.⁴² 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.⁴² 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=α2σ(κe+κL)$ for 5 nm ≤ d_G ≤ 1 μm at h = 3 × 10¹⁹ cm⁻³ (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¹⁹ cm⁻³) as well as the experimental κ_L = 0.85 Wm⁻¹ K⁻¹ were subsequently scaled using these normalized transport trends to calculate nanostructured ZT for each dopant at 0.3 at. % concentration.

Our first-principles calculations on Bi₂Te₃ 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₂Te₃ 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₂Te₃. 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₂Te₃ 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.

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..

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³TEC EFRC, Award No. DE-SC0001299/DE-FG02-09ER4657.