A review of recent progress in thermoelectric materials through computational methods

Abstract

Reducing our overwhelming dependence on fossil fuels requires groundbreaking innovations in increasing our efficiency in energy consumption for current technologies and moving towards renewable energy sources. Thermoelectric materials can help in achieving both goals. Moreover, because of recent advances in high-performance computing, researchers more increasingly rely on computational methods in discovering new thermoelectric materials with economically feasible performance. In this article, significant thermoelectric materials discovered through these computational methods are systematically reviewed. Furthermore, the primary computational tools that aid the design of the next-generation thermoelectric materials are introduced and discussed. These techniques include various levels of density functional theory, electronic transport simulations, and phonon calculations.

Introduction

More than 70% of the global energy consumption is lost during conversion as waste heat [1], released to the environment at a wide range of temperatures (300–1000 K). This wasted energy originates from numerous sources such as mining and smelting plants, factories, electricity generators, air-conditioning and lighting appliances, pipelines, vehicles’ engines, brakes and exhausts, and even the human body. Harvesting this waste heat is of paramount importance for improving energy consumption efficiency. One promising solution is the thermoelectric (TE) effect [2], which governs the direct conversion of heat to electricity in a TE material. Additionally, TE materials can further reduce the CO₂ emissions, and enable extensive deployment of battery- and grid-free devices, wearable electronics, the Internet of Things (IoT), and much more. The performance of a thermoelectric material at a given temperature T is quantified by its figure of merit (ZT). ZT is defined as $ZT = S_{{{\text{Coeff}}}}^{2} \sigma T/\kappa$, where S_Coeff is the Seebeck coefficient, σ is the electrical conductivity, and κ is the thermal conductivity. κ is made of two contributions: the electronic thermal conductivity (κ_e) and the lattice or phonon thermal conductivity (κ_L). Maximizing ZT in thermoelectric materials requires three criteria satisfied: first, a large S_Coeff, to obtain the highest electric potential difference for a given temperature gradient; second, a low κ, to delay the thermal equilibration between the heat source and the environment; last, a high σ, to generate the highest current for the thermo-voltage.

To construct a TE harvesting unit, one must interface an n-type TE material to a p-type TE material, forming an n–p TE junction. Then, by putting one side of the TE junction in contact with a hot surface (waste heat source) and the other side with a cold surface (ambient), electricity is generated on the cold side. Several couples of the TE junctions are usually coupled together, both serially and in parallel, forming a TE module, shown in Fig. 1, to increase the amount of harvested energy. Such a configuration is commonly achieved by cutting n- and p-type bars from sintered pellets [3], or preparing n- and p-type thin films epitaxially deposited on an inflexible single-crystal substrate [4]. The primary appeal of TE generators is the fact that they are capable of converting waste heat directly into electrical energy without any moving parts. As a result, research thrived and plenty of TE materials (alloys, oxides, and nitrides) were discovered. After nearly a century of research and development, however, the practical use of TE devices is still limited to niche applications [5,6,7]. There are three significant obstacles against the widespread application of TE harvesters: (i) the low conversion efficiency and low chemical stability of TE materials; (ii) the toxicity and criticality of elements used in TE materials; and (iii) the mechanical rigidity of TE modules that impedes application in wearables.

Fig. 1

Adapted with permission from Ref.[12]. Copyright 2018 IEEE

Structural design of wearable thermoelectric harvesting device: a schematic view; b close-up view of one thermoelectric module.

In the search for high-performance thermoelectric materials, it has now become evident that the traditional materials’ engineering approach based on judicious choice of materials and synthesis conditions via trial-and-error would only result in little improvement of the TE performance [8, 9]. As a result, a theory-guided approach in the search for suitable high-performance thermoelectric materials with favorable electronic and thermal transport is becoming common [10,11,12]. In this review, we, therefore, systematically examine and present the recent developments in thermoelectric research propelled by theoretical and computational methods. The theoretical tools reviewed here can be divided into two broad categories: (a) electronic-structure calculations and that deal with S, σ, and κ_e; and (b) lattice dynamics or phonon calculations that determine κ_L. Although these two types of calculations are interrelated, they, nonetheless, describe fundamentally different aspects of the materials’ properties. As such, different computational tools have been developed for electronic and lattice dynamics calculations, of which both are covered here.

The need for computational insight

The second law of thermodynamics fundamentally dictates the upper limit efficiency of any thermodynamic engine, including the thermoelectric power generators. The efficiency (η) of a solid-state power generator, as a function of ZT, operating between hot and cold reservoirs with temperatures T_H and T_C is given by the following nonlinear expression [13]:

$$\eta = \frac{{T_{{\text{H}}} - T_{{\text{C}}} }}{{T_{{\text{H}}} }}\left[ {\frac{{\left( {1 + {ZT}_{{\text{M}}} } \right)^{1/2} - 1}}{{\left( {1 + {ZT}_{{\text{M}}} } \right)^{1/2} + \left( {T_{{\text{C}}} /T_{{\text{H}}} } \right)}}} \right] , \quad T_{{\text{M}}} = \left( {T_{{\text{H}}} + T_{{\text{C}}} } \right)/2.$$

(1)

The challenge in achieving high ZT ($S_{{{\text{Coeff}}}}^{2} \sigma T/\kappa$) materials lies in the interdependency of the κ, σ, and S_Coeff in any given material system, as all are functions of the carrier concentration (n). For instance, in an isotropic n-type semiconductor, σ can be expressed as:

$$\sigma \left( {T,n} \right) = e \cdot n\left( T \right) \cdot \mu_{e} \left( {T,n} \right),$$

(2)

in which e is the unit charge and μ_e (T,n) is the electron mobility. μ_e (T,n) is governed by a multitude of scattering mechanisms that often vary with T and carrier concentration, in a complex manner. n(T), in turn, is commonly controlled by extrinsic doping. S_Coeff can be generally expressed by the formula [14]:

$$S_{{{\text{Coeff}}}} \left( {T,n} \right) = - \frac{{k_{{\text{B}}} }}{e} \cdot \left( {\frac{5}{2} - {\ln}\frac{n\left( T \right)}{{N_{{\text{D}}} }}} \right),$$

(3)

in which k_B is the Boltzmann constant, and N_D is the concentration of n-type (donor) dopants. One should note that Eq. 3 only holds at relatively large n(T) as S_Coeff diverges when n(T) approaches zero [15]. As can be easily deciphered, σ in thermoelectric semiconductors can be raised by carrier doping. As an undesired side effect, however, the additional carriers decrease S_Coeff, resulting in a meager or no net benefit to the ZT.

Furthermore, the relation between κ and n is even more complicated. κ_e is related to σ by the Wiedemann–Franz law [16]. As a result, any carrier doping unfavorably boosts the electronic share of the thermal conductivity, while κ_L, which is the dominant thermal conduction mechanism over a wide range of carrier concentrations in semiconductors, depends on the density, heat capacity, and thermal diffusivity of which all are influenced by doping. For all these complications, κ is often phenomenologically modeled as a temperature- and concentration-dependent power law [17]:

$$\kappa = \kappa_{300} \left( {\frac{T}{{300 \,{\text{K}}}}} \right)^{{\sigma_{\kappa } }} \cdot \left( {1 + \left( {\frac{n\left( T \right)}{{n_{{{\text{ref}}}} }}} \right)^{{\beta_{\kappa } }} } \right).$$

(4)

Here, κ₃₀₀ is the thermal conductivity at 300 K, while n_ref, σ_κ, and β_κ are constants to be determined by fitting experimental measurements. In the intrinsic range, both electrons and holes contribute toward κ, constituting a bipolar contribution towards κ. The heavier an n-type semiconductor is doped, however, the more generous the κ_e contribution is towards κ.

A reflection upon the complex interdependence of S_Coeff, σ, and κ, one recognizes the challenges of designing high-efficiency thermoelectric materials. The challenges arising from carrier doping are schematically presented in Fig. 2. In addition to the electronic requirements, one must also consider other constraints such as the non-toxicity and the abundance of the elements used (usually associated with the price and sustainability), mechanical robustness, moist resistance, flexibility (for wearables), and the low dopant solubility in some cases, to mention a few. As a result, the typical experimental approach in finding new thermoelectric materials inspired by search within a composition or structural space similar to known materials is very restrictive by nature. Consequently, by resorting to computationally guided materials search, design, and discovery, one may accelerate the pace of the progress in TE research beyond the limit imposed by the judicious and serendipitous experiments.

Fig. 2

adapted from Refs. [18] and [19], while κ was adapted from Ref. [20]

The dependence of S_Coeff, σ, and κ on carrier concentration (n) for polycrystalline Pb_1−xSn_xTe at 300 K. ZT reaches its maximum only for a very narrow window of n. A change in temperature shifts all these curves in a non-trivial manner, resulting in a maximum ZT at a different n. Furthermore, synthesis protocols may drastically change the presented values. S_coeff and σ values were

Electronic calculations

The electronic description

The efficiency of the thermoelectric energy conversion is directly related to the materials’ non-dimensional figure of merit, which grows linearly with the electrical conductivity [21]. The essential quantities are then σ, together with the S_Coeff, which are determined by the electronic structures of the materials and by carrier-scattering processes. Currently, there are two main approaches for calculating the electron transport properties and S_Coeff: the semiclassical formulas based on Boltzmann transport theory; and the fully quantum–mechanical approach described in the Kubo–Greenwood theory [22,23,24]. An alternative, fully quantum–mechanical approach has also been developed by Landauer [25]. Here, the scattering processes are modeled as if they build charges and fields inside the materials, which represents an alternative view to that of Kubo’s method in which currents are formed in response to a given external electric field. The Landauer methodology intimately resembles the experimental viewpoint, where one usually imposes an external current and measures the resulting potential drop due to the electron scattering, which is an important tool in guiding intuition when studying mesoscopic transport properties.

The semiclassical Boltzmann transport equation, however, has gained appreciable popularity in recent years due to the ease of use and its integration with a diverse range of density functional theory (DFT) codes. The semiclassical Boltzmann transport equation, as implemented in the prominent BoltzTrap [24] and BoltzTrap2 [26] packages, relies only on the DFT calculated band and k-dependent quasiparticle energies (ε) as input. It then extrapolates the transport properties within the rigid band approximation [24]. Furthermore, this method assumes a T and ε independent relaxation time (τ) for electronic scatterings which results in a tractable and straightforward form of the equations for S_Coeff, σ, and κ_e. The constant relaxation-time approximation, despite its simplicity, predicts S_Coeff values that match well with experiments and is widely adopted in the high-throughput theoretical search of novel TE materials. Table 1 summarizes some of the recent predictions using the BoltzTrap package.

In brief, the Boltzmann transport equation describes the general statistical behavior of a system out of thermodynamic equilibrium in terms of local equilibrium distributions. The Boltzmann transport equation is derived by assuming random motion of particles—in this case, electrons—within a medium, leading to the formulation of a current (j) in terms of conductivity tensors:

$$j_{i} = \sigma_{ij} E_{j} + {N}_{ij} {\nabla}_{j} T,$$

(5)

in which σ_ij and Ν_ij are the elements of the electrical and thermoelectric conductivity tensors. From the band structure, within the rigid band approximation, the conductivity tensors along the reciprocal lattice direction α and β can be obtained as follows:

$$\sigma_{\alpha \beta } \left( {i,{\mathbf{k}}} \right) = e^{2} \tau_{{i,{\mathbf{k}}}} v_{\alpha } \left( {i,{\mathbf{k}}} \right)v_{\beta } \left( {i,{\mathbf{k}}} \right),$$

(6)

where e is the electron charge, τ_i,k is the relaxation time, v(i,k) is the group velocity, i is the band index, and k is the reciprocal (crystal momentum) vector. Here, v(i,k) can be calculated from the band structure by taking the derivative of the band energy function (f) over k. Here, the relaxation time (T_i,k) describes how quickly the system returns to thermodynamic equilibrium after scattering events. T_i,k plays an important role in determining the transport properties such as electron mobility, electrical conductivity, thermal conductivity, and Seebeck coefficient. T_i,k principally depends on both the band index (i) and the k vector direction. However, several studies of the direction dependence have shown that T_i,k can be considered, with a good approximation, direction-independent in almost all known crystals [27]. A further approximation of T_i,k to a constant value leads to an expression for the Seeback coefficient that is also independent of the relaxation time. However, an important restriction for the relaxation-time approximation to be valid is that the time the distribution function takes to relax to its equilibrium is independent of the distribution function itself and the externally applied electric field. By integrating σ_αβ over the energy (ε) space, the transport tensor can be obtained as:

$$\sigma_{\alpha \beta } \left( {T,\mu } \right) = \frac{1}{{\Omega }}\int {\sigma_{\alpha \beta } \left( \varepsilon \right)\left[ { - \frac{{\partial f_{\mu } \left( {T,\varepsilon } \right)}}{\partial \varepsilon }} \right]{\text{d}}\varepsilon } ,$$

(7)

$${N}_{\alpha \beta } \left( {T,\mu } \right) = \frac{1}{{eT{\Omega }}}\int {\sigma_{\alpha \beta } \left( {\varepsilon - \mu } \right)\left[ { - \frac{{\partial f_{\mu } \left( {T,\varepsilon } \right)}}{\partial \varepsilon }} \right]{\text{d}}\varepsilon } ,$$

(8)

$$\kappa_{\alpha \beta } \left( {T,\mu } \right) = \frac{1}{{e^{2} T{\Omega }}}\int {\sigma_{\alpha \beta } \left( {\varepsilon - \mu } \right)^{2} \left[ { - \frac{{\partial f_{\mu } \left( {T,\varepsilon } \right)}}{\partial \varepsilon }} \right]{\text{d}}\varepsilon ,}$$

(9)

in which μ is the electronic chemical potential, and Ω is the unit cell volume. Finally, the Seebeck coefficient can be easily computed as follows:

$$S_{ij} = (\sigma^{ - 1} )_{\alpha i} {N}_{\alpha j} .$$

(10)

The physical limitation to the semiclassical Boltzmann theory is that it is only valid for the first-order terms, while it more or less fails at higher orders. This shortcoming is due to the assumption of chaotic molecular motion in Boltzmann’s collision integral—an approximation that is only accurate up to first order [28]. However, the equivalence of the two most critical theoretical methods for the calculation of electronic transport within the linear regime has been explicitly demonstrated, and both of them allow the calculation of electric transport and S_Coeff from first-principles (ab initio) calculations, namely DFT. The excellent agreement between the calculated results and the experimental data reported in the literature validates the validity of the transport theories when coupled with ab initio computational methods [29,30,31].

Wannier functions

The maximally localized Wannier function representation allows us to naturally introduce the ground-state electronic structure into the lattice Green’s function approach that will be the basis for the evaluation of the Landauer quantum conductance [42, 43]. This approach is computationally efficient, can be easily implemented as a postprocessing step in standard electronic-structure DFT calculations, and allows us to link the electronic transport properties of a device to the nature of the chemical bonds directly, providing an insight into the mechanisms that govern electron flow at the nanoscale. Owing to the exponential localization of the Wannier functions obtained from planes waves, accurate integrations over the Brillouin zone can be achieved at a very moderate computational cost. Moreover, the analytical expression for the band derivatives in the Wannier basis resolves any issues that may occur when evaluating derivatives near band crossings. The accuracy of the results on a coarse k-point reciprocal-space grid using Wannier functions was confirmed by Pizzi et al., who performed calculations of the electrical conductivity on the binary and ternary skutterudites CoSb₃ and CoGe_3/2S_3/2 [44]. Moreover, the semiclassical Boltzmann transport theory based on maximally localized Wannier functions with constant relaxation-time approximation was applied to examine the transport properties of SrTiO₃ [45]. The results for the transport properties show the prospect of thermoelectric presentation for SrTiO₃ at different temperatures from 300 to 800 K as a function of chemical potential, which plays an important role. Varying the temperature in terms of chemical potential shows a remarkable enhancement in the TE performance for the properties of thermoelectric materials. In the considered cubic phase, the chemical potential of SrTiO₃ is equal to the sum of all atomic component chemical potentials in the crystal [45].

Table 1 The predicted S_Coeff values for the TE materials, calculated with the BoltzTrap package that have recently reported in the literature

In search of flat bands

The band dispersion, in other words, the slope of the bands near band edges, is of great importance for predicting materials’ TE performance. Steeper slopes indicate stronger orbital interactions and faster carrier mobility. On one hand, for stronger orbital interactions, the band becomes more disperse in energy. On the other hand, weaker orbital interactions correspond to flatter bands. This simple model elucidates the fact that dispersed bands lead to higher carrier mobility, while flat bands indicate more localized carriers. However, there is a class of fermion systems with a diverging density of states (DOS). These type of materials, such as CuAlO₂ [46], exhibit a dispersionless spectrum that have precisely zero dispersion, the so-called flat band [47]. Flat bands can emerge in strongly interacting condensed matter systems and layered compounds with integer-valued pseudo-spin [48,49,50,51,52].

Ishii et al. found that the electric current flowing along the chain suddenly increases the magnitude when a minimal electric field is applied to the system. This large current originates from the symmetry breaking of the large energy degeneracy of flat-band states in the Kagome lattice, which is nearly independent of the size of the Kagome‐lattice chain itself [53]. Under flat-band conditions, Ishii et al. showed that carriers in high-lying minibands could propagate for several superlattice periods without being scattered. These carriers follow the miniband description of transport. In a uniform electric field, the minibands break down into localized states immediately. These states participate in conduction. One strategy to achieve large S_Coeff and large electron transport is via low-energy electronic bands containing both, flat and dispersive parts, in different regions of crystal momentum space. This approach is known as the pudding-mold band structure, of which examples are shown in FIG. 3 [54, 55].

Fig. 3

Adapted with permission from Refs. [46, 54, 55]. Copyrights 2013, and 2019 American Physical Society and 2007 Physical Society of Japan

The red curves in a, b show pudding-mold band structure in CuAlO₂ [46] and BaPdS₂ [54], respectively. Schematic figure of the electronic structure within the pudding-mold structure shown in (c) and compare to the one of a typical metal (d). Moreover, e, f represent the band dispersion [55]. Here, t₁ is the hopping integral between neighboring sites, and ε is the quasiparticle energy. The dotted gray lines mark the position of the chemical potential for a given temperature T.

When we consider the problem of electron transport from a more general viewpoint, the literature suggests that most layered compounds are particularly interesting because of their unusual electronic and physical properties, in addition to their desirable thermoelectric performance [56]. For example, through electronic-structure calculations using DFT, a large and anomalous anisotropic thermoelectric transport properties were predicted in d⁰-electron layered complex nitrides AMN₂ (where A and M are metal cations). Subsequently, it was proposed that these compounds may constitute a novel and promising class of TE materials [57,58,59,60]. On one hand, because of the dimensionality of the electronic structures, anisotropic electronic transport properties have also been observed in layered oxides La_xSr_1−xCuO₄ and Sr₂RuO₄ [61,62,63]. On the other hand, complex layered nitrides in the form AMN₂ (A = Sr or Na; M = Zr, Hf, Nb, Ta), were analyzed. The Sr containing compounds SrZrN₂ and SrHfN₂ were found to have three-dimensional electronic structures and isotropic thermoelectric transport properties, despite their α-NaFeO₂-type layered crystal structure [64]. This apparent contradiction can be explained by the fact that the bottom of the conduction bands in these compounds is composed of Zr $4{d}_{{z}^{2}}+4{d}_{xz}+4{d}_{yz}$ and Hf $5{d}_{{z}^{2}}+5{d}_{xz}+5{d}_{yz}$, along with Sr $4{d}_{{x}^{2}-{y}^{2}}+4{d}_{xy}$ orbitals. Here, the A and M ions jointly provide a full set of the d orbitals. However, although the electronic structures of NaNbN₂ and NaTaN₂ are also three-dimensional, NaNbN₂ and NaTaN₂ were found to have weak anisotropic thermoelectric transport properties. This weak anisotropy arises from the NbN₂ and TaN₂ layers whose conduction band minimums are composed of Nb $4{d}_{{z}^{2}}+4{d}_{xz}+4{d}_{yz}$ and Ta $5{d}_{{z}^{2}}+5{d}_{xz}+5{d}_{yz}$ along with N 2p_z orbitals with no contribution from the Na orbitals. Understanding and elucidating the origin of anisotropic transport properties in layered compounds are of great significance for cultivating this new class of TE compounds, not only for thermoelectric materials but also for realizing novel physical properties [65].

Strongly correlated thermoelectric materials

Spin-driven Seebeck effect

Most conventional thermoelectric materials, of which we reviewed some in the previous section, are semiconductors, rather than metals, and, consequently, have low intrinsic electrical conductivity. That is because although metals are the ideal conductors, they, nonetheless, have a diminishingly small S_Coeff. The carrier symmetry at the Fermi level (E_F) can explain metals’ small Seebeck coefficient, as it implies the number of hot electrons above E_F and cold empty states below E_F are roughly equal. When the number of diffused holes and electrons under a temperature difference is roughly the same, the net Seebeck effect is minimal due to cancelation. Semiconductors circumvent this problem by having only one type of majority carriers. There is, however, one class of thermoelectric materials that noticeably deviates from the conventional paradigm, i.e., the strongly correlated thermoelectrics. Thermoelectric materials such as Ca₃Co₄O_9+δ [66], Fe_1+yTe_0.6Se_0.4 [67], and Na_xCoO₂ [68] are excellent electrical conductors—Na_xCoO₂ has a high carrier concentrations (n) which is in the order of ~ 10²¹ to 10²² cm⁻³ [69]. In these materials, spin entropy flow drives the Seebeck effect. For instance, in Na_xCoO₂, the susceptible spin population is approximately equal to the hole population and remains constant with decreasing T. Na_xCoO₂, therefore, has a conductivity similar to a metal and a magnetic behavior similar to a frustrated insulator. Such compounds are known as Curie–Weiss metals [70, 71]. Wang et al. [72] demonstrated that at T = 2 K, a magnetic field completely suppresses S_Coeff to nil. The disappearance of S_Coeff in Na_xCoO₂, when an external magnetic field aligns all the spins, implies that the Seebeck effect originates from the spin entropy flow carried by the holes in the Curie–Weiss phase. In this section, therefore, we review the structural and electronic characteristics of Na_xCoO₂—the most investigated such compound—as a representative of the strongly correlated thermoelectric material.

Na _x CoO ₂

Sodium cobaltate (Na_xCoO₂) is a strongly correlated oxide [73, 74] with a fairly high thermoelectric figure of merit (ZT) at temperatures around 800 K for x > 0.5 [75]. As shown in Fig. 4, the Na_xCoO₂ lattice is made of alternating Na and edge-sharing CoO₂ octahedral layers stacked along the c direction constituting a quasi-two-dimensional system. In Na deficient (x < 1) Na_xCoO₂, each Na vacancy creates a hole converting a Co³⁺ ion into Co⁴⁺. As a result, Na_xCoO₂ exhibits p-type conductivity with n inversely proportional to x. Moreover, the coexistence of Co³⁺ and Co⁴⁺ ions, or charge disproportionation, in the CoO₂ layers generates a competitive Seebeck potential through spin entropy flow [68, 72]. Furthermore, the mobile Na⁺ ions diffusively scatter the phonons [70, 76], decreasing κ_L to ~ 0.01 W cm⁻¹ K⁻¹ at ~ 1000 K [77]—this κ_L value is ~ 100 times smaller than that of ZnO at the same temperature [78].

Fig. 4

A schematic presentation of the Na_0.75CoO₂ structure constructed out of a desodiated 4a × 2a × 1c supercell of the P6₃/mmc Na₁CoO₂ primitive cell. Co and O ions occupy the Wyckoff 2a and 4f sites of the hexagonal lattice structure, respectively. In Na_xCoO₂, some of the Na ions occupy 2b (Na1) sites which share basal coordinates with Co, while the rest of of the Na ions occupy 2c (Na2) sites which share basal coordinates with O. There are two layers of Na ions (and two layers of CoO₂) in the unit cell. Na ions at Z = 0 are presented in red, while Na ions at Z = 0.5 are presented in blue

The electrical transport properties of Na_xCoO₂ have also been reported for different concentrations of Na in several experimental studies [70, 79, 80]; however, the results are inconsistent at times and even contradictory. For example, Kawata et al. [79] found that for x > 0.55, the electrical resistivity (ρ) of Na_xCoO₂ increases along with increasing x. Motohashi et al. [80], however, found the opposite trend for ρ as a function of x. More recently, Foo et al. used the chemical de-intercalation from Na_0.77CoO₂ to change the Na concentration and found a metal–insulator transition of the charge ordering in Na_0.5CoO₂ single crystal. A similar effect, however, could not be detected in the polycrystalline Na_xCoO₂ prepared by solid-state reaction. In another work, Cui et al. [81] reported a comprehensive study of the structure and electrical transport properties of Na_xCoO₂ (x = 0.52, 0.56, 0.60, 0.68, 0.72, 0.76, 0.8). Here, all the samples had metallic conduction, with ρ reaching a minimum at around x = 0.64. Moreover, a jump in ∂ρ/∂T was observed at the temperature corresponding to the dimensional crossover of the electronic transport properties. Some of the discrepancies can be explained by the structural differences between Na_xCoO₂ prepared by solid-state reaction and by chemical de-intercalation from Na_0.77CoO₂. For instance, it was shown that the distance between Co in the CoO₂ layers, which may depend on the synthesis protocol, is critical in determining the electrical transport property of Na_xCoO₂—and more generally, of all the Na_xTMO₂ (TM = Cr, Mn, Fe, Co, Ni)—where this critical distance R_c that determines whether the compound is insulating or metallic [82]. As seen here, the origin of this difference in observations is still an open question [83, 84] and sets a challenge for theoretical analysis.

Spin entropy flow

The high-temperature S_Coeff in Na_xCoO₂ can be expressed by Koshibae’s equation [85, 86] for strongly correlated materials, which itself is a modified form of Heikes formula [87, 88]:

$$S_{{{\text{Coeff}}}} \left( {T \to \infty } \right) = - \frac{{k_{{\text{B}}} }}{e}{\ln}\left[ {\frac{{g\left( {{\text{Co}}^{3 + } } \right)}}{{g\left( {{\text{Co}}^{4 + } } \right)}}\frac{x}{1 - x}} \right].$$

(11)

Here, g refers to the degeneracy of the electronic configuration in ions A and B (here A = Co³⁺ and B = Co⁴⁺), and it is equal to all the distinct possible ways which the electrons can be arranged in the orbitals of the ions they occupy. x is the ratio of Co⁴⁺ over all Co ions, which is proportional to the carrier (hole) concentration, and finally, k_B and e are the Boltzmann’s constant and the electron’s charge, respectively. This concept is schematically demonstrated in Fig. 5. In Na_xCoO₂, the six electrons of the low-spin Co³⁺ have a $t_{{2{\text{g}}}}^{6} e_{{\text{g}}}^{0}$ configuration for which there is no other degenerate state that is entropy S = ln[g(Co³⁺)] = ln1. In contrast, for the low-spin Co⁴⁺ ion, since one electron is removed, there are six degenerate configurations. Consequently, S = ln[g(Co⁴⁺)] = ln6. As an electron hops from a Co³⁺ ion to a Co⁴⁺ ion, an entropy flux moves opposite to the electric current flow (hole); this phenomenon is referred to as spin entropy flow [89]. According to Eq. 11, the S_Coeff for Na_0.75CoO₂ (a common representative of the Na_xCoO₂ family with metallic conduction) is 249 μV K⁻¹, which is quite comparable with the experimental measurements of ~ 200 μV K⁻¹ at 800 K [90].

Fig. 5

Reproduced with permission from Ref. [91]. Copyright 2020 Institute of Physics

The schematics of spin entropy flow in Na_xCoO₂. S and g are the entropy and the electronic degeneracy per site, respectively.

Xiang et al. attempted to explain Na_xCoO₂’s high S_Coeff value using standard Boltzmann transport theory combined with spin-polarized DFT [92]. They found that S_Coeff was smaller when the system was set to be magnetically polarized, and interpreted the results as an alternative explanation for the suppression of S_Coeff in a magnetic field. According to their interpretation, the high S_Coeff and its suppression in high magnetic fields could be understood using the conventional transport theory applicable to other thermoelectric materials and that the unique thermoelectric properties of Na_xCoO₂ are due to its unusual band structure as manifested in the narrow manifold of t_2g bands in this system. However, the dependence of S_Coeff on Co’s spin state (high spin vs. low spin) [93] still points to the spin entropy flow as the primary cause for high S_Coeff values.

What level of DFT?

Density functional theory (DFT) is one of the most common theoretical tools for probing strongly correlated TE materials such as Na_xCoO₂. Within the DFT framework, the most basic approximation for the exchange–correlation functional, i.e., the local density approximation (LDA), in which the exchange and correlation energy terms are fitted to those of a homogenous electron gas [94], fails to describe the electronic structures in highly correlated systems realistically [95]. The next most popular functional based on the general gradient approximation (GGA) [96], often, offers only a minor improvement [97]. The shortcoming stems from the fact the both LDA and GGA functionals tend to delocalize electrons over the lattice space. As such, each electron feels an average of the Coulombic potential. Consequently, for highly correlated materials, such as Na_0.75CoO₂, the significant Coulombic repulsion between localized electrons that results in charge disproportionation, are not correctly represented by LDA and GGA [83]. The thermoelectric effect, nonetheless, stems from charge disproportionation and needs an accurate theoretical description.

A standard solution for avoiding the delocalization problem of LDA and GGA functionals is to add an ad hoc Hubbard-like localized term to the LDA and GGA functionals. This approach is known as DFT + U (commonly referred to as LDA + U and GGA + U). Practically, in the DFT + U method [98, 99], one applies a Hubbard (U) and an exchange (J) terms to the localized d (or f) electrons. The delocalized electrons (s and p) are left unaltered as LDA and GGA correctly describe these delocalized electrons. To avoid the double-counting of the correlation part for localized electrons, another term—called the double-counting correction—is subtracted from the Hamiltonian. U − J, known as U_eff, determines the level of localization within the DFT + U implementation.

The U_eff value is not usually known a priori. Therefore, it is often chosen to replicate some experimental measurements [100]. For this reason, multiple values are usually reported in the literature for a given compound depending on the experimental values used as a benchmark, such as the band gap, formation enthalpy, or magnetic phase. In the case of Na_xCoO₂, among the reported U and J values in the literature [101,102,103], all reproduce the charge disproportionation [104, 105] of the Co ions more or less accurately. That is to clearly distinguish Co²⁺, Co³⁺, and Co⁴⁺ from one another in terms of magnetization and partial density of states (PDOS). For instance, the more common value of U_eff = 4 eV is at a midpoint between the lower values of ~ 3 eV implemented Co³⁺ [106] and the higher values of ~ 5 eV implemented Co⁴⁺ [107,108,109].

The criticality of U _eff

As presented in the previous section, the exact choice of U_eff for Na_xCoO₂ remains somehow ambiguous. Here, we examine how the results of the DFT calculations with various U_eff value for Na_0.75CoO₂ (as a representative for the Na_xCoO₂ compound family) differ from one another. The problem to address is the arrangement of Na ions (or equally Na vacancies) in the unit cell. Figure 6a shows the Na_0.75CoO₂ obtained with the LDA level of the theory [110]. In this structure, the unit cell contains a $\sqrt 3 a \times 4a \times 1c$ supercell in which the Na ions were arranged in rows of 3 × Na2, progressively followed by a Na vacancy, a 3 × Na1, and another Na vacancy. This structure has equal numbers of Na2 and Na1 ions (Na1/Na2 = 1). Figure 6(b) shows the Na arrangement obtained with U_eff = 4 eV (GGA + U formalism) [111]. Here, the irreducible unit cell has dimensions of $2\sqrt 3 a \times 2a \times 1c$ in which one-third of the Na ions occupied Na1 sites while two-third occupied Na2 sites (Na1/Na2 = 0.5) as the Na vacancies create a multi-vacancy cluster as marked by green enclosures. Moreover, Na ions at the Na2 sites were ever slightly off the oxygen axes conforming to the H1 structure [112]. Figure 6c shows the Na arrangement in Na_0.75CoO₂ calculated with U_eff = 5 eV (GGA + U formalism) [113]. Here, the ratio of Na1/Na2 is ~ 0.44, where the Na ions in Na1 sites form ordered tri-Na1 droplets. The margin of the stability of the structure at Fig. 6c was, however, only ~ 2 meV lower than the competing structures which the authors had considered, implying that phase separation might occur at room temperature.

Fig. 6

Adapted after Refs. [110, 111, 113]. Copyrights 2005 American Physical Society, 2015 World Scientific, and 2008 American Institute of Physics

The schematic presentation of Na ion arrangements in Na_0.75CoO₂ calculated with different U_eff values. a LDA calculations based on Zhang et al. [110]. The open circles denote those Na ions at Z = 0, while solid circles denote those Na ions at Z = 0.5 within the unit cell. b GGA + U (U_eff = 4 eV) calculations based on Assadi et al. [111]. Here, the color scheme is similar to that in Fig. 4. c GGA + U calculation based on Hinuma et al. [113] with U_eff = 5 eV. Here, the Na ions at Z = 0 are at mirror locations of those at Z = 0.5. The blue and gold spheres rather denote Na1 and Na2.

We can infer that the choice of U_eff is demonstrably critical in describing the arrangement of Na ions in Na_0.75CoO₂. Consequently, to validate the theoretical results, an accurate comparison with the experiment is required. This task is, however, easier said than done. Na is volatile, and its concentration in the sample may be smaller than the nominal value [114]. Furthermore, deciphering the real atomic environment at a given lattice point in Na_xCoO₂ experimentally, which is strongly dominated by the coordination environment, can be determined using the analysis of X-ray absorption near-edge spectra. However, such a measurement often is limited by the energy range of the instruments, rendering the interpretation of the results less than deterministic [115]. Obtaining a well-calibrated U_eff value for Na_xCoO₂, therefore, remains a difficult task. Utilizing higher level theory such as Hartree–Fock/DFT hybrid functionals [116] or GW approximation [117] may offer more robust theoretical tools for probing Na_xCoO₂ and similar materials. Unfortunately, these methods are more computationally demanding than the standard or + U corrected DFT calculations.

Phonon calculations

As mentioned above, the efficiency of a TE generator depends on the ZT of the TE material, where κ can be further separated into lattice (κ_L) and electronic (κ_e) thermal conductivities. The latter is correlated with σ as shown by the Wiedemann–Franz law:

$${{\kappa_{{\text{e}}} } \mathord{\left/ {\vphantom {{\kappa_{{\text{e}}} } \sigma }} \right. \kern-\nulldelimiterspace} \sigma } = LT,$$

(12)

where L is the Lorenz number. Since κ_L is the only independent variable in ZT, lowering κ_L has been the most effective way to increase ZT in the past decades. An ideal TE material has the so-called “phonon-glass electron-crystal” (PGEC) property [118], with a high power factor ($PF=\sigma {S}_{\text{Coeff}}^{2}$) and a low κ_L. In this section, we will briefly summarize recent advances and challenges related to the ab initio description of κ_L.

Phonons and thermal conductivity

Debye first used quantized normal mode of the collective atomic vibrations to describe the specific heat of crystalline solids, giving birth to the concept of the phonon [119]. The phonon description of lattice vibrations has been spectacularly successful and has become a pillar of solid-state physics. In a harmonic oscillating lattice, phonons never interact, there is no lattice expansion, and κ_L approaches infinity. In a real lattice, the interatomic potential has anharmonic terms, which leads to interactions between phonon quasiparticles, thus resulting in a finite κ_L and lattice thermal expansion (some intriguing materials exhibit negative thermal expansion within a temperature range [120,121,122,123]). Peierls adapted the Boltzmann theory for the transport of gases to phonon “gas” to describe their population distributions driven out of equilibrium by a temperature gradient, formulating the Peierls–Boltzmann transport (PBTE) equation for determining phonon lifetimes and κ_L [124, 125]. The solution of the PBTE, a complex set of coupled integrodifferential equations, is a complicated task by itself. A variety of approximations have thus been introduced to develop an understanding of phonon driven thermal properties. For example, one can mention the Callaway and Klemens models of intrinsic phonon scattering and conductivity [126, 127]. Furthermore, Slack et al. developed a simple model for κ_L and a simple set of rules for its interpretation based on crystal complexity, average atomic mass, bonding strength, and anharmonicity [128, 129]. These are quantified in terms of a material’s Debye temperature and average Grüneisen parameter, which gives a measure of the speed of sound of the heat-carrying phonons and anharmonicity, respectively.

With the recent advances in theoretical and computational techniques and the exponentially increasing computational power, phonon band structure, group velocities, lifetimes, and thermal conductivity can be calculated from ab initio without ad hoc adjustable parameters. The solution of the PBTE has been achieved so substantially beyond the typical relaxation-time approximation, from variational methods [130] to self-consistent iterative procedures [131] and compressive sensing [132]. Moreover, realistic interatomic potentials for a variety of materials can be obtained from DFT [133]. The harmonic and anharmonic interatomic force constants (IFCs), as in Fig. 7 and Eq. 13, can be determined numerically either from DFT electronic-structure calculations based on supercell perturbations method [134] or from the linear response methods within density functional perturbation theory (DFPT) [135]. In the anharmonic case, as shown in Fig. 7, the interatomic force is not proportional to the displacement, due to the third, fourth, and even higher order terms in the Taylor expansion of the potential energy V:

$$\begin{aligned} V = & \sum\nolimits_{ij,\alpha \beta } {\frac{{\partial^{2} V}}{{\partial u_{\alpha }^{i} \partial u_{\beta }^{j} }}\Delta u_{\alpha }^{i} \Delta u_{\beta }^{j} } + \sum\nolimits_{ijk,\alpha \beta \gamma } {\frac{{\partial^{3} V}}{{\partial u_{\alpha }^{i} \partial u_{\beta }^{j} \partial u_{\gamma }^{k} }}\Delta u_{\alpha }^{i} \Delta u_{\beta }^{j} } \Delta u_{\gamma }^{k} \\ & \quad + \sum\nolimits_{ijkl,\alpha \beta \gamma \sigma } {\frac{{\partial^{4} V}}{{\partial u_{\alpha }^{i} \partial u_{\beta }^{j} \partial u_{\gamma }^{k} \partial u_{\sigma }^{l} }}\Delta u_{\alpha }^{i} \Delta u_{\beta }^{j} } \Delta u_{\gamma }^{k} \Delta u_{\sigma }^{l} . \\ \end{aligned}$$

(13)

Fig. 7

Comparison of a harmonic and an anharmonic oscillator

Nowadays, fully ab initio calculations of lattice thermal conductivity have become standard practice, owing to the availability of computational power and open-source code packages including but not limited to ShengBTE [136], Phono3py [137], PhonTS [131], and almaBTE [138]. However, these standard computations represent perfect and infinite single crystals. Real materials are finite and with imperfections, point defects like vacancies and interstitials, dislocations, and grain boundaries. In TE materials, defects are typically desired to tune the electronic band structure via doping and provide important scattering centers for phonons to lower κ_L. Defects add considerable complexity to the phonon calculations [139] and the prediction of κ_L in realistic materials. Ab initio calculations of phonon–phonon and phonon–defect interactions on the same footing are just becoming possible recently.

Realistic phonon calculations

A single site point defect changes the local structure and chemical bonding around it, providing a scattering center for phonons decreasing the thermal conductivity. This thermal resistance mechanism can be accounted, to some extent, within the above-described software packages for calculating κ_L, often demonstrating excellent agreement with the measured data for high-quality crystals. However, the theory assumes that the point defects are randomly distributed in the lattice within the diluted limit, meaning that defects can scatter phonons without any correlation between them. As seen in Fig. 8, in the engineered materials with high point-defect concentration, this independent scattering model is questionable [140,141,142]. The inadequacy of this first-order perturbative description of point-defect scattering was shown in the recent work when comparing results from the first-order approximation with those of fully ab initio Green’s function calculations of phonon scattering by vacancies and nitrogen substitution in diamond [141]. This work showed that including the variation of IFC induced by point defects and going beyond the first-order perturbations give much more reliable phonon-defect scattering prediction than the simple mass perturbation model.

Fig. 8

Adapted with permission from Refs. [141] and [142]. Copyrights 2014 and 2018 American Society of Physics

a κ_L of diamond as a function of defect concentration. N and V stand for nitrogen and vacancy, respectively. The insets show the [100] projection view of the N impurity and vacancy defects after relaxation. Green: N atom. Blue: first coordination shell. Yellow: second and third coordination shells. b κ_L of In_1−xGa_xAs calculated at 300 K with the VCA and with the special quasirandom structure approach.

Further works have combined first-principles Green’s function and PBTE methods to describe and predict κ_L in a variety of materials with different defect types (intrinsic vacancies and antisite, and substitutional dopants) [143, 144]. Among the intrinsic defects, vacancies are generally found to have the most potent effects on phonons scattering, due to the magnitude of the missing mass from vacancy centers and structural distortions. However, the “dilute limit” approximation is still assumed in most of the typical state-of-the-art phonon-defect calculations. Various methods have been employed for high concentrations of defects and coherent scattering mechanisms, such as the virtual crystal approximation (VCA) [145], supercell unfolding [142], and ab initio molecular dynamics (AIMD) simulations [146, 147].

One avenue in the search for low κ_L and high ZT is to examine very complex materials [148]; the ones that possess a variety of bonding patterns, rattling atoms in covalent cage structures, and many low-frequency optic phonon branches for the scattering of heat-carrying acoustic phonons. Computationally, such crystal complexity brings about significant challenges. In fact, ab initio PBTE methods have been recently used to explore complex unit cells of, for example, PbRb₄Br₆ or Cu₁₂Sb₄S₁₃ [132, 149]. The numerous vibrational degrees of freedom in the complex unit cells make the DFT calculations of the IFCs computationally tricky and increase the phase space for phonon–phonon interactions. These numerical bottlenecks would become less problematic as the computational power increases, offering a better understanding of the relevant physical mechanisms governing κ_L, and allowing better approximations to be made without much loss of precision in the near future.

A part of the state-of-the-art ab initio thermal transport methodologies is based on the phonon quasiparticle scattering picture within the quantum perturbation theory. Intrinsic phonon–phonon interactions are derived from the perturbations coming from higher order anharmonicity of the interatomic potential (derivatives of the interatomic potential with respect to each atomic displacement, see Eq. 13). For materials with extreme behaviors (such as strongly anharmonic) or in extreme conditions [like high T or pressure (P)], the ab initio methods derived from lowest order perturbation theory could show discrepancies between measured and calculated phonon transport properties. Materials that are about to or have undergone phase transition around or below the temperature of interest pose a crucial challenge as their lattice is highly anharmonic. Explicitly incorporating the thermal effects onto the anharmonic IFCs, beyond the quasiharmonic approximation, therefore, becomes necessary. Furthermore, including the fourth-order anharmonicity has been shown to be essential for κ_L in Si [150] at high temperature, and in PbTe, even at room temperature [151]. Molecular dynamics (MD) simulations of phonon transport include all orders of the anharmonicity. MD, however, treats the atomic interactions classically; and AIMD remains prohibitively costly to tackle realistic system sizes.

Describing the lattice vibration and thermal transport at extreme T and P and near phase transitions driven by T or P will be the challenge for future ab initio phonon transport calculation and even to DFT in general. Phonons in materials near phase transitions can have localization as the mean-free path becomes comparable to their wavelength, constituting the Ioffe–Regel regime [152]. In this circumstance, a wave can no longer be readily defined. A localized vibrational state can become hybridized with propagating phonon states, resulting in resonant-like extended states with diffusive motion, not propagating, in nature. Perhaps, new quasiparticle picture like diffusons/locons [153, 154] in disordered materials will be needed to describe the lattice vibration phenomena in these extreme cases.

Recently, Hellman et al. have developed an accurate method, called the temperature-dependent effective potential (TDEP) technique, to determine the temperature-dependent anharmonic free energy, based on ab initio molecular dynamics followed by a mapping onto an effective model Hamiltonian describing the lattice dynamics. This method has been applied to various materials with dynamic instability, such as bcc Zr [155], and some strongly anharmonic materials, such as PbTe [156] and SnSe [157], which happen to be good thermoelectrics due to the very low thermal conductivity. Moreover, Zhou et al. calculated the electron–phonon interaction in the presence of soft-modes in SrTiO₃ perovskite from ab initio using (anharmonic) TDEP phonons [158]. The stochastic self-consistent harmonic approximation (SSCHA) is another novel development in the theoretical framework to study from ab initio the anharmonic properties of solids. SSCHA has been successfully applied to study the second-order structural phase transition and temperature-dependent anharmonic phonons primarily in various materials, such as the ferroelectric SnTe and GeTe [159], and the thermoelectric SnSe [159], as shown in Fig. 9.

Fig. 9

Adapted with permission from Ref. [159]. Copyright 2019 American Society of Physics

Harmonic and anharmonic phonons in the Lorentzian approximation for SnSe. The length of the blue bars corresponds to the linewidth (full length of the line is the full width at half maximum). The calculations were performed with LDA on the experimental structure.

For advanced TE materials, high-performance improvement needs to combine optimization of both the electronic and vibrational properties of the materials as these are tightly coupled. Thus, state-of-the-art numerical simulations of the TE materials need to incorporate electron–phonon interactions to improve our understanding of κ_L S_Coeff, and σ as coupled properties. Numerical techniques have been built and developed recently for efficient calculation of electron–phonon interaction matrix elements from DFPT and quantum perturbation theory [160, 161]. These ab initio electron–phonon coupling has been examined and used to study transport phenomena in doped TE semiconductors [158, 162,163,164,165] and metals [166,167,168]. Rather than solving the transport equation for electron and the phonon separately, the solution of their coupled distribution equations would be the focus of future calculations. These calculations shall be critically important in understanding the phonon-drag effects in TE materials, and the highly out-of-equilibrium electron and phonon dynamics in the ultrafast laser-based thermal spectroscopies.

Interfacial systems

Effective control of interfaces in thermoelectric materials is critical for improving their efficiency. In general, interfaces in thermoelectric devices are divided into two main types: (a) grain boundaries within single-phase materials and solid solutions; and (b) interfacial systems formed between different materials, e.g., superlattices or epitaxial layers. Figure 10 shows a schematic representation of grain boundaries in polycrystalline materials and multilayered heterostructures. Grain boundaries constitute the interfaces between grains in polycrystalline bulk materials. These grains may exhibit different atomic structures, crystal orientations, or different compositions. The role of grain boundaries is especially crucial in nanostructured materials, which are composed of small crystal grain sizes and exhibit a high density of grain boundaries. However, despite the numerous experimental observations of grain boundaries, their implications on transport properties are not entirely understood. Theoretical simulations are, therefore, a powerful tool to study the role of grain boundaries and interfacial systems in thermoelectric materials at the nano- and sub-nano-scale. Computer simulations can be used to discuss the stability of specific interfacial systems, in which vacancies or defects are commonly present (and serve as scattering centers, as explained in the previous section), the materials’ behavior under mechanical strain or shear deformations, and their corresponding TE performance. Such theoretical investigations are, however, computationally expensive and have just recently become popular within the thermoelectric community. In this section, we will review some of the landmark studies in the field, focusing on the computational techniques employed.

Fig. 10

Schematic representation of grain boundaries in a polycrystalline materials and b multilayered heterostructures. Variations in colors represent different crystallographic orientations in grains or different material compositions in solid solutions or heterostructures. Grain boundaries are denoted with black lines

Nanostructuring and grain boundaries

From the theoretical perspective, electron and phonon transport in polycrystalline one-dimensional graphene nanoribbons and nanotubes was investigated by density functional tight-binding (DFTB) theory combined with non-equilibrium Green functions by Lehmann et al. [169]. The introduction of lattice imperfections such as dislocations was modeled by pentagons, heptagons, and octagons embedded within the hexagonal carbon nanostructures (Fig. 11a) which is suggested as an efficient method for enhancing the TE effect. The computed S_Coeff exhibit values ~  − 0.4 mV/K and ~  − 0.6 mV/K in polycrystalline nanoribbons and nanotubes, respectively, in contrast to S_Coeff > 0 in the ideal graphene structures. Variations of the electronic and thermal properties in polycrystalline carbon-based systems give rise to an improved ZT by several orders of magnitude at room temperature.

Fig. 11

Adapted with permission from Refs. [169, 173, 180]. Copyrights 2015 and 2017 American Physical Society, and 2010 American Institute of Physics

Structures of a polycrystalline graphene nanoribbon with pentagon (blue shaded area), heptagon (orange), and octagon (yellow) defects obtained by joining two layers with an inclination angle of 24° [169], b Bi₂Te₃ one quintuple layer-thick slab (Te and Bi atoms are presented with orange and purple spheres, respectively) [173], and c p-type LaNiO₃/SrTiO₃(001) superlattice with two p-type interfaces (La, Sr, and O atoms are purple, green, and red, respectively. Light and dark blue octahedra are centered around Ti and Ni ions, respectively) [180].

Phosphorene (composed of a few-layer-thick black phosphorous) is another bidimensional material that exhibits highly anisotropic electrical and thermal conductances in single-crystal form. A combined DFT and self-energy correction (GW) study [170] demonstrated that the preferred directions for κ and σ in a phosphorene monolayer are orthogonal to one another. This results in an anisotropic ZT that is large along the armchair direction, with ZT above 1 at 300 K and above 2 at 500 K, both for 60 ps relaxation time. However, in the experimentally processed material, the presence of defects may reduce the conductance and diminish the anisotropies.

In Si nanowires, simulations carried out by the ensemble Monte Carlo technique (EMC) [171] showed that the electrical conductivity decreases along with the cross section of the wire, while, on the other hand, S_Coeff increases as the wire cross section decreases. The calculated ZT in nanowires is 20–40 times more substantial compared to the bulk Si. The decrease in the κ_L due to strong phonon-boundary scattering can explain the remarkable enhancement in ZT. The above-described studies show that the formation of grain boundaries can be used to improve the thermoelectric efficiency by maximizing S_Coeff while decreasing the resistivity (ρ) and κ. On the other hand, grain boundaries can also affect the conductance negatively, reducing the anisotropy observed in single-crystalline counterparts.

Interfaces and heterostructures

Besides the recent progress in improving the performance of thermoelectric materials by modifying grain boundaries in nanostructured compounds, superlattices (formed by alternate stacking of two materials) have also shown potential to enhance the efficiency of TE energy conversion devices [172]. Ab initio DFT calculations showed that the slab formed by one atomic quintuple layer of Bi₂Te₃ (Fig. 11b) has the capacity to increase ZT to 7.15, which is ten times larger than for the bulk structure at room temperature [173]. The significant enhancement in ZT is caused by the changes in the valence band density from quantum confinement within the thin film.

Saha et al. studied HfN/ScN [174] and ZrN/ScN [175] metal/semiconductor superlattices constructed by stacking layers of both materials alternatively. The structure of these model interfaces and the computed electrostatic potentials as a function of the distance are depicted in Fig. 12. The p-type Schottky barrier height for the HfN/ScN of 0.77 eV is larger than the 0.34 eV in ZrN/ScN superlattices. The n-type barrier of the HfN/ScN superlattices is estimated to be 0.13 eV (versus 0.56 eV for ZrN/ScN), which suggests that these superlattice models are degenerate for electron transport along the cross plane. The lower barrier in HfN/ScN indicates the metallic conduction and lower S_Coeff compared to the ZrN counterpart. On the other hand, HfN/ScN superlattices present better κ_L than the ZrN-based models. These results point out the trade-off between S_Coeff and σ, suggesting that simultaneously increasing both parameters is still challenging in two-dimensional geometries.

Fig. 12

Adapted with permission from Refs.[174, 175]. Copyrights 2012 Institute of Physics and 2011 American Institute of Physics

a Metal–nitride interface structure formed between b ScN and ZrN [175] or c ScN and HfN [174]. b and c represent the planar average electrostatic potential (green dashed line) as a function of perpendicular distance from the ScN/ZrN and ScN/HfN interfaces, respectively. The red dashed line indicates the lattice-plane oscillations, which are filtered with the macroscopic averaging technique [182]. The vertical dashed black line represents the lattice planes.

Ab initio calculations have also been used to discuss the stability of solid-state TE with high cooling power density, such as epitaxial Ni/Bi₂Te₃, NiTe/Bi₂Te₃, Co/Bi₂Te_3, and CoTe₂/Bi₂Te₃ interfaces [176]. In these models, the metallic part is compressed/expanded by up to ± 10% from its equilibrium structure to match the Bi₂Te₃ lattice parameters. All these interfaces were found to form stable Ohmic contacts. Other ab initio studies of epitaxial CoSb₃/TiCoSb [177] and CoSb₃/Ti [178] interfaces discussed the stability, electronic structure, and mechanical properties of these model structures. These publications [177, 178] examined the atomic-level mechanical strain failure mechanisms in CoSb₃/Ti-based interfacial systems and found that the fracture takes place in the CoSb₃ region, in line with experimental studies. The brittleness of the covalent Sb–Sb bond was compared with the more stable (ionic) Co–Sb bond in interfacial [177] and pure CoSb₃ models [179].

The study of epitaxial SrTiO₃/SrTi_1−xNb_xO₃ (0 ≤ x ≤ 0.5) modeled by all-electron DFT showed that Nb doping induces metallicity into the originally insulating SrTiO₃ [181]. The Nb doping yields an increase in the carrier concentration while decreasing the thermoelectric response of the heterostructure. The exceptional agreement between the theory and experiment presented in this study demonstrates the high reliability of first-principles calculations for predicting S_Coeff within a wide-doping range. In SrTiO₃ interfaces forming epitaxial superlattices with LaNiO₃ (Fig. 11c), DFT + U calculations described how coupling interfaces of opposite polarity generate an electric field in the superlattice [180]. The analyses of the electronic transport and thermoelectric properties showed the anisotropic behavior of this interfacial model, which presents a large positive cross plane S_Coeff ~ 135 μV K⁻¹ and a ZT of 0.35 through the p-type superlattice at room temperature.

Understanding the atomic-level interfacial structure, mechanical stability, and electronic features are essential for the design of novel multi-layered materials with an enhanced thermoelectric ZT. These studies showed that the TE efficiency of superlattices could be largely improved compared to their bulk counterparts. The critical parameters for increasing the ZT in superlattices are optimizing the Schottky barrier for enhancing S_Coeff while retaining a high σ and reducing the cross plane κ at the interface. The mechanical properties are also an essential factor to consider for engineering applications, since the different thermal expansion coefficient at both sides of the interface may give rise to accumulated thermo-mechanical stresses that detrimentally form interfacial cracks.

Outlook

So far, we have discussed how Boltzmann transport equation, as implemented in BoltzTrap [24, 26] can offer solid predictions of the electronic transport properties of the thermoelectric materials. Furthermore, Peierls–Boltzmann transport (PBTE) equation, as implemented in ShengBTE code [136], can determine the phonon lifetimes and κ_L. These tools, coupled with the existing material structure databases, either based on computational [183,184,185] or experimental [186, 187] compilations, can automate the search for brand new thermoelectric materials through high-throughput calculations within a particular composition or structure space [188, 189]. Furthermore, interfaces based on the compounds in these databases can be automatically generated and screened using the MPInterfaces tool [190]. The DFT and transport calculations can be automated using Atomate software package [191] to run automatically, repair, and report the outcome. The outcome then can be searched for suitable intrinsic descriptors such as stability, band gap, flat-band feature, and S_Coeff to identify the most promising materials. Recently compounds such as Cu₁₂Sb₄S₁₃, Cu₂₆V₂Ge₆S₃₂ [192], ZnSnSb₂, AgInS₂, AgGaSe₂, AgInSe₂, and LiInTe₂ [193] have been predicted to be feasible thermoelectric materials via high-throughput calculations.

High-throughput calculations, although very powerful and predictive, come with some disadvantages. First, the uniformity of the computational settings among many chemical compositions under investigation may cause convergence issues in some compounds. Consequently, the generated data set may not be as accurate for all compounds or even miss some data points for which the DFT calculations catastrophically fail [194]. Moreover, regarding the transport properties, the assumption of a constant relaxation time across a set of compositions may not be realistic and thus is qualitative at best [195]. These problems restrict the capability of most high-throughput calculations in identifying the best thermoelectric materials. Consequently, the promising compounds discovered through such calculations should be further scrutinized with a higher level of theory. For instance, the use of the AMSET package [196] for the electronic transport calculations may remedy some of the shortcomings with the BoltzTrap method as it does not assume a constant relaxation time. AMSET, however, currently relies on a one-dimensional model of the averaged band structure, which limits its predictive power for materials with highly anisotropic band features.

An alternative to high-throughput calculations is to restrict the theoretical search to a handful number of compounds at each investigation—like the studies in Refs [36, 39] of Table 1, and Refs. [197, 198]—and fine-tune and customize the theoretical settings for the compounds under investigation. As one can see, the balance between the number of compounds investigated and the level of the accuracy is a trade-off to be determined by the constraints on the resources and the design of the search for new thermoelectric materials [199].

Conclusions

This review has highlighted the critical role that advanced computational techniques play in the understanding and design of thermoelectric materials. For instance, band structure calculations based on density functional theory with appropriate functional can aid in predicting the conventional thermoelectric materials operating based on a pudding-mold-shaped band near the Fermi level, such as CuAlO₂, or strongly correlated thermoelectric materials that operate based on Heikes formula for materials with mixed valency such as Na_xCoO₂. Similarly, density functional theory can be employed to study the more complex interface effects at grain boundaries or heterostructures, for instance, in ScN/ZrN. Furthermore, the wavefunctions calculated by density functional theory can be used to predict the electronic transport (electrical conductivity, Seebeck coefficient, and electronic contribution to thermal conductivity) using Wannier Functions at minimum cost. Phonon dispersion calculations that can elucidate the lattice contribution to the thermal conductivity can be calculated with the Peierls–Boltzmann transport equation. As such, theoretical calculations can predict all properties affecting the thermoelectric performance of a given material or system. Consequently, the computational thermoelectric investigation can complement and assist in the analysis of experimental measurements; for instance, defect concentration, crystal structure, and grain boundary orientations. Moreover, such computational studies can explain the vital atomic-scale features and provide a fundamental understanding of the thermoelectric effect in novel materials for which detailed experimental characterization had not yet completed. Most importantly, such computational investigation can play a predictive role in the design of robust new thermoelectric materials.