## Theoretical model for the Seebeck coefficient in superlattice materials with energy relaxation

Vassilios Vargiamidis, 1, ∗ Mischa Thesberg, 2 and Neophytos Neophytou 1 1 School of Engineering, University of Warwick, Coventry, CV4 7AL, UK 2 Institute for Microelectronics, Technical University of Vienna, Vienna, A-1040, Austria

We present an analytical model for the Seebeck coefficient S of superlattice materials that explicitly takes into account the energy relaxation due to electron-optical phonon (e-ph) scattering. In such materials the Seebeck coefficient is not only determined by the bulk Seebeck values of the materials but, in addition, is dependent on the energy relaxation process of charge carriers as they propagate from the less-conductive barrier region into the more-conductive well region. We calculate S as a function of the well size d , where carrier energy becomes increasingly relaxed within the well for d &gt; λ E , where λ E is the energy relaxation length. We validate the model against more advanced quantum transport simulations based on the nonequilibrium Green's function (NEGF) method and also with experiment, and we find very good agreement. In the case in which no energy relaxation is taken into account the results deviate substantially from the NEGF results. The model also yields accurate results with only a small deviation (up to ∼ 3%) when varying the optical phonon energy glyph[planckover2pi1] ω or the e-ph coupling strength D 0 , physical parameters that would determine λ E . As a first order approximation, the model is valid for nanocomposite materials and it could prove useful in the identification of material combinations and in the estimation of ideal sizes in the design of nano-engineered thermoelectric materials with enhanced power factor performance.

PACS numbers: 73.20.-r, 73.43.-f, 72.10.-d

## I. INTRODUCTION

When a temperature gradient is applied in a solid material with free electronic carriers, a voltage gradient arises as carriers migrate from the hot side to the cold side. The strength of this thermoelectric effect is quantified by the Seebeck coefficient S , which is defined as the ratio of the voltage difference ∆ V to the temperature difference ∆ T . The absolute value of S is referred to as the thermopower.

The Seebeck coefficient is central to the performance of a thermoelectric (TE) material, which is quantified by its TE figure of merit ZT = σS T/κ 2 , where σ is the electronic conductivity, T is the temperature, and κ is the thermal conductivity. The product σS 2 is known as the power factor ( PF ). Although Bi Te 2 3 and PbTe are traditionally the most extensively studied TE materials, over the last several years various other materials have been explored with respect to their TE performance, such as transition-metal dichalcogenides [1-3], phonon-glasselectron crystals [4], half-Heuslers [5], tin selenide [6], etc. Most of these materials exhibit ZT above 1, primarily due to the reduction of their thermal conductivity [7].

Superlattices and nanocomposite materials are also currently being explored aiming to achieve even higher TE performance [8-12]. This is due to two reasons. First, they usually cause reduction of the phonon thermal conductivity to ultra-low values as a result of extensive phonon-boundary scattering [8]. In fact, this is considered as one of the most effective ways to enhance

∗ Electronic address: V.Vargiamidis@warwick.ac.uk

TE performance. Second, such nanostructures quite often also cause increase in the Seebeck coefficient [13-16], and interestingly, in some cases cause increase in the PF as well [17-19].

The design of superlattice TE materials requires extensive theoretical and computational modeling. For this purpose, several methods have been employed previously. Some of these methods adopt semi-classical approaches using the Boltzmann transport equation (BTE) where the effects of grain boundaries are treated as a scattering mechanism with some relaxation time [20-22]. However, in order to capture key aspects of the physics -such as tunneling, non-equilibrium carrier relaxation and confinement - in a single setting (especially as the material feature sizes shrink to the nanoscale), the use of a quantum transport method such as the nonequilibrium Green's function (NEGF) method [23, 24] is necessary. However, these methods are either complex or time consuming, or both. On the other hand, it is important to be able to determine fast and relatively accurately the Seebeck coefficient of superlattice (SL) materials using simpler models. This can be especially useful in experimental settings, and will guide nanostructured designs that will allow for high Seebeck coefficients, potentially high PF s as well.

Simplified models to describe the Seebeck coefficient in superlattices (and nanocomposites) exist and are widely used in the literature. The simplest way is to describe the overall Seebeck coefficient as the weighted average of the Seebeck coefficients of the well and the barrier regions, with the weighting factor being the length of each region [15]. In an additional step, in order to satisfy the continuity of heat transfer, the individual components are also weighted not only by the length of the regions,

but also by the inverse of their thermal conductivity [25]. The individual coefficients are usually obtained from the Boltzmann transport formalism separately for each region. Other works, on the other hand, use the energy dependencies of the coefficients from the BTE and by assuming thermionic emission over the potential barriers [14, 26]. A phenomenological model has also been proposed [20] for the calculation of the Seebeck coefficient of nanocomposites where interface potential barriers due to grains have been included and the effect of various scattering mechanisms was examined. In more elaborate cases, wave solutions of electronic transport are employed, which account for the formation of mini-bands as well, which are then included in transport [16],although such mini-bands would be weakened in the presence of electron-phonon interaction [10]. In principle, however, there is an intermediate region, where the Seebeck coefficient transits from the barrier into the well and viseversa, as electrons relax their energy (or gain energy) to go from one region to the other. In structures where the energy relaxation mean-free-path is comparable to the well size, this region becomes important. In fact, we have shown in the past that it is the existence of this region that allows for signicant power factor improvements in SLs and nanocomposites [10, 17, 25]. Thus, this region needs to be properly described in compact models that apply to the new generation nanocomposite TE materials, and currently, no compact model exists (despite the importance of it being evident in large scale simulations [9, 10, 18]).

In this paper, we develop a simple analytical model for the Seebeck coefficient of a channel with embedded SL barriers for energy filtering, which mimics either a SL or a nanocomposite to first order approximation, or any material in which carrier transport alternates between potential barriers and wells. Using the average energy of the current flow, and taking the energy relaxation length λ E as calculated from NEGF we derive an expression for S . The relaxation length is generally used to describe the relaxation of the carrier energy along the transport direction due to phonon emission, and is therefore a measure of the distance that the relaxation process occurs. For carriers flowing over barriers and relaxing into wells, it essentially denotes the region where the individual attributes of conductivity and Seebeck coefficient intermix. It is directly connected with the more familiar energy relaxation time , τ E , which can be calculated for different scattering mechanisms [27] and is known for many materials [28].

We present results for the case in which the Fermi level E F is ≈ k B T below the barrier height, and also for the case where E F = V B . We find that the results for S as a function of well size d are in very good agreement with the corresponding results of NEGF. Further, in NEGF we alter the optical phonon energy glyph[planckover2pi1] ω and e-ph coupling strength D 0 , which are the physical parameters that affect the energy relaxation length λ E . For all extracted λ E 's that we consider, the model predicts accurately the

FIG. 1: (Colour online) (a) Schematic representation of a superlattice nanostructure. A channel of length L is connected to ideal reservoirs (left and right contacts) with chemical potentials and temperatures µ 1 , T 1 and µ 2 , T 2 , respectively. The well size is d and the barrier thickness is b . We also define t = b + d (see Appendix A). (b) Average energy of the current flow 〈 E x ( ) 〉 as defined in Eq. (5) along the channel with SL barriers calculated with NEGF. The black lines represent the potential barriers and the blue line represents the position of the Fermi level E F . The color map indicates the current flow I ( E,x ), with yellow indicating high, and green low current density.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-d7aff31f-5989-41d3-81e5-8ce1357e13c6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-42ae55ac-3e24-45df-8c88-da58ba61d782.jpg)

dependence of S on these physical parameters.

The rest of the paper is organized as follows. In Sec. II we derive the model for the Seebeck coefficient. Then, in Sec. III, we present, analyse, and validate the results with those from NEGF. In Sec. IV we validate the results with an experiment, while in Sec. V we summarize and conclude.

## II. SEEBECK COEFFICIENT

## A. Model without energy relaxation

We consider a nanostructure composed of two different materials in which charge carriers propagate through low and high energy regions, or potential wells and barriers and abstract such a material to that of a onedimensional (1D)-like system of potential barriers. A schematic representation of such a SL nanostructure is shown in Fig. 1(a). A channel of length L is connected to two contacts (left and right), which are ideal reservoirs in equilibrium with µ 1 , T 1 and µ 2 , T 2 their chemical po-

tentials and temperatures, respectively, while b and d are the barrier thickness and well length, respectively.

Note that, although the system is conceptualized as a SL and only one direction of transport is considered, the final analytical model that we derive below depends only on bulk Seebeck coefficients and energy relaxation lengths of the constituent materials. Thus, the final model is considered to be agnostic to issues of dimensionality and valid for all dimensional structures. Furthermore, we also argue that this model should, on average, be valid for nanocomposite materials. This is because the primary conceptual difference between nanocomposites and SLs is that in a nanocomposite the barriers only have an average spacing of d drawn from a statistical distribution, rather than the fixed rigid spacing of a superlattice. However, since the central crux of the model rests on the dominant effect of carrier relaxation physics, we would expect any coherent phenomena related to quantum reflections, resonances, etc., whose existence separates the nanocomposite and SL cases, to be negligible regardless.

The most commonly used model describes the total Seebeck coefficient of such a system, S sys , as a combination of the Seebeck coefficients of the well and barrier regions, S W and S B , respectively. This is derived from [25]

$$S _ { s y s } = \frac { 1 } { \Delta T } \int _ { 0 } ^ { L } S ( x ) \left ( \frac { d T _ { L } } { d x } \right ) d x, \quad \ \ ( 1 ) \quad \ e l e c { \prod _ { \substack { \text{$\text{$\text{$per}$} } } } }$$

where S x ( ) is the local Seebeck coefficient, and ∆ T is the lattice temperature difference along the channel. Note that this expression, strictly speaking, depends on the lattice temperature, T L , however here we only consider the temperature of carrier flow, T . We take T L = T , a point justified in Refs. [25, 29] since optical phonon scattering plays the dominant role in energy relaxation (encouraging equilibrium with the phonon bath). Through Fourier's Law, we express the temperature gradient in a barrier (well) region as ( dT/dx ) B(W) = J/κ B(W) , where κ B and κ W are the thermal conductivities in the barrier and well regions, respectively, and J is the heat flux. Also, we express ∆ T as

$$\Delta T = \left ( \frac { J } { \kappa _ { \text{B} } } \right ) L _ { \text{B} } + \frac { J } { \kappa _ { \text{W} } } \left ( L _ { \text{W} } + L _ { \text{W} } ^ { \prime } \right ), \quad \left ( 2 \right )$$

where L B = nb is the sum of all barrier thicknesses with n the total number of barriers, L W = ( n -1) d is the sum of all well lengths, and L ′ W = 2 d ′ is the total length of the two wells at the ends of the channel. Note that x 0 = d ′ [see Fig. 1(a)]. The reason why the terminating regions are treated separately is to allow direct comparison with NEGF simulation later on, where this is necessary. Using Eq. (2), we can express Eq. (1) as

$$S _ { s y s } ^ { \text{no rel} } = \frac { ( S _ { B } L _ { B } / \kappa _ { B } ) + ( S _ { W } \tilde { L } _ { W } / \kappa _ { W } ) } { ( L _ { B } / \kappa _ { B } ) + ( \tilde { L } _ { W } / \kappa _ { W } ) }, \quad ( 3 ) \quad \text{In E} \\ \text{curr} \quad \text{even}$$

where ˜ L W = L W + L ′ W .

Although Eq. (3) describes well the composite Seebeck coefficient in macroscale materials, when the feature sizes of the composite phases are scaled below a few tens of nanometers, this model is inadequate. The reason is that in the vicinity of the materials' interfaces, the Seebeck coefficient does not abruptly change from S B to S W , but carriers have to gradually relax their energy (and also their momentum) to the value imposed by the equilibrium conditions of each material. In fact, this takes place within a distance determined by the energy relaxation length λ E [25, 29], see Fig. 1(b). It is important to note here that any possible PF improvement in such materials originates from the intermixing of the high Seebeck coefficient of barrier material S B with the high conductivity of the well material σ W , thus making these regions very important in composite nanostructures. We emphasize that, just by considering separately the individual (bulk) Seebeck coefficients of the two regions S B and S W , it is not easy to achieve PF improvements compared to the maximum of the two PF s of the individual barrier or well.

In Sec. IIB, we develop a simple and relatively accurate model for the Seebeck coefficient of SL materials, taking into account the energy relaxation process due to e-ph scattering. We validate the model against the results from NEGF. In the NEGF simulations we consider only electron-optical phonon scattering mediated through the e-ph coupling strength D 0 , which is the mechanism most responsible for energy relaxation.

## B. Model with energy relaxation

In order to derive an analytical model for the Seebeck coefficient of a nanocomposite system consisting of potential barriers and wells as shown in Fig. 1(a), we assume that the charge carriers are fully relaxed in the barrier regions, but in the well regions the carriers undergo a relaxation process, which is quantified by λ E . Apart from the optical phonon energy glyph[planckover2pi1] ω and the deformation potential D 0 , a key parameter that influences the relaxation process is the well size d , as discussed below.

The x dependent (local) Seebeck coefficient is given as

$$S ( x ) = \frac { \langle E ( x ) \rangle - E _ { F } } { q T }, \text{ \quad \ \ } ( 4 )$$

where q is the carrier charge ( q = -| | e for electrons and q = + | e | for holes), E F is the Fermi level, and 〈 E x ( ) 〉 is the average energy of the current flow along the x direction (propagation direction), defined as

$$\langle E ( x ) \rangle = \frac { 1 } { J } \int _ { E = E _ { c } } ^ { E = \infty } I ( E, x ) E d E. \quad \quad ( 5 )$$

In Eq. (5), I ( E,x ) is the energy and position resolved current, while J = ∫ I ( E,x dE ) = constant. Note that even though the current is constant along the channel at each cross section, its energy is not constant, i.e., the

FIG. 2: (Colour online) Average energy of the current flow 〈 E x ( ) 〉 (dashed, black lines) along the channel with (a) a single barrier, and (b) no barrier and E F = 0 05eV, calculated . with NEGF. The solid (red) line is fitting of Eq. (A5) in order to extract the energy relaxation λ E (see Sec. III). 〈 E 〉 B and 〈 E 〉 W are the average energies on the barrier, and in the well under equilibrium, respectively. The color map indicates the current flow I ( E,x ), with yellow indicating high, and green low current density.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-7d3a075f-f4d0-43d0-8fdd-c0b1fe226804.jpg)

charge carriers can gain or lose energy as they propagate. This happens in the presence of inelastic scattering (optical phonons). We emphasize that knowledge of 〈 E x ( ) 〉 in Eq. (4) allows one to determine the Seebeck coefficient S x ( ) and vice versa, regardless the complexities of the nanostructure.

We substitute now S x ( ) from Eq. (4) into Eq. (1) and we obtain

$$S _ { s y s } = \frac { 1 } { q T \Delta T } \int _ { 0 } ^ { L } \left ( \langle E ( x ) \rangle - E _ { F } \right ) \left ( \frac { d T _ { L } } { d x } \right ) d x. \quad ( 6 ) \quad \text{$\dots$}$$

In Eq. (6), in order to make contact with NEGF simulation, it is necessary to consider 〈 E x ( ) 〉 in the well regions as being different from that in the wells at the channel ends (i.e., those that are close to the contacts) as shown in Appendix A. Further, in the barrier regions 〈 E x ( ) 〉 is taken to be constant (see Eq. (A1) of Appendix A), which turns out to be a good approximation. Accordingly, Eq. (6) can be expressed as a sum of four terms, each of which pertains to different region, i.e., well regions, wells at the left and right channel ends, and barrier regions, ˜ s W , ˜ , ˜ s L s R , and ˜ s B , respectively. The total Seebeck coefficient can then be expressed as

$$S _ { s y s } = \tilde { s } _ { L } + n \tilde { s } _ { B } + ( n - 1 ) \tilde { s } _ { W } + \tilde { s } _ { R }, \quad \ ( 7 )$$

where n is the total number of barriers in the channel. Each one of the four terms in Eq. (7) is of the same form as Eq. (6), but each one with different form of average energy 〈 E x ( ) 〉 and with different integration region.

Starting with ˜ s B , we use Eq. (A1) for a single barrier, i.e.,

$$\langle E ( x ) \rangle = \langle E \rangle _ { \text{B} },$$

where 〈 E 〉 B is the average energy on the top of a barrier [see Fig. 2(a)]. Thus,

$$\tilde { s } _ { B } = \left ( \frac { b } { \Delta T } \right ) \frac { J } { \kappa _ { B } } S _ { B },$$

where

$$S _ { \text{B} } = \frac { \langle E \rangle _ { \text{B} } - E _ { \text{F} } } { q T }. \text{ \quad \quad } ( 1 0 )$$

In order to evaluate S B in Eq. (10), we use NEGF to evaluate first the average energy 〈 E 〉 B on top of a single barrier placed in the channel, as shown in Fig. 2(a). However, in general the value of S B could be extracted from bulk experimental values. As it turns out S B is actually one of three parameters needed to determine the total Seebeck coefficient S sys .

For the calculation of ˜ s W we use Eq. (A3) for a single well (see rationale in the Appendix), as also seen in Ref. [29],

$$\tilde { \ } g e n \quad & \langle E ( x ) \rangle = ( \langle E \rangle _ { B } - \langle E \rangle _ { W } ) \\ & \times \left ( e ^ { - ( x - x _ { 0 } - b ) / \lambda _ { E } } + e ^ { ( x - x _ { 0 } - t ) / \lambda _ { E } } - e ^ { - d / \lambda _ { E } } \right ) + \langle E \rangle _ { W }, ( 1 1 )$$

where 〈 E 〉 W is the average energy within the well under equilibrium and x 0 + b &lt; x &lt; x 0 + . The value of t 〈 E 〉 W can be extracted from NEGF as the average energy of the current flow for a pristine channel (as is done here), or taken from experimental measurement. The result is then that of an 'effective' well, as shown in Fig. 2(b). Using Eq. (11), we can express ˜ s W in the form

$$\tilde { s } _ { \text{W} } = \left ( \frac { d } { \Delta T } \right ) \frac { J } { \kappa _ { \text{W} } } S _ { \text{W-relax} }, \text{ \quad \ \ } ( 1 2 )$$

where S W-relax is obtained using Eq. (6) and the form of 〈 E x ( ) 〉 given in Eq. (11). It is given by

$$\text{gous} \\ \text{annel} \\ \text{$\mathcal{E}(x)$} \\ x \text{ A} ),$$

In Eq. (13) we made use of Eq. (10) for S B and the corresponding relation for S W . The value of S W is obtained

from 〈 E 〉 W of the pristine channel shown in Fig. 2(b). Note that in the two limits λ E → 0 and λ E glyph[greatermuch] d , Eq. (13) yields S W-relax → S W and S W-relax → S B , respectively, i.e., the corresponding bulk values of the Seebeck coefficient in the well and barrier regions, as expected. A similar calculation for ˜ s R using Eq. (A5) yields

$$\tilde { s } _ { \text{R} } = \left ( \frac { d ^ { \prime } } { \Delta T } \right ) \frac { J } { \kappa _ { \text{W} } } S _ { \text{W-relax} } ^ { \prime }, \quad \quad ( 1 4 ) \quad \text{par} \\ \quad \text{is e}$$

where

$$S _ { W - \text{relax} } ^ { \prime } = S _ { W } + ( S _ { B } - S _ { W } ) \left ( \frac { \lambda _ { E } } { d ^ { \prime } } \right ) \left ( 1 - e ^ { - d ^ { \prime } / \lambda _ { E } } \right ). \ ( 1 5 ) \quad \text{between} \\ \quad \text{therm}$$

Note that d ′ = x 0 (see Fig. 1(a)). The contribution of ˜ s L is identical to that of ˜ s R , i.e., ˜ s L = ˜ s R . Thus, finally, from Eq. (7) the total Seebeck coefficient of the system takes the form

$$S _ { s y s } & = & \frac { \int _ { \substack { \text{$L_{B}S_{B}/\kappa_{B})+(L_{W}S_{W-relax}/\kappa_{W})+(L_{W}S_{W-relax}/\kappa_{W})$} \\ (L_{B}/\kappa_{B})+(\tilde{L}_{W}/\kappa_{W})$} \\ (L_{B}/\kappa_{B})+(\tilde{L}_{W}/\kappa_{W})$} \\ ( 1 6 )$$

In Eq. (16) we notice that -aside from the geometric factors of L W , L B , etc. - the Seebeck coefficient of the SL material is determined in terms of only three parameters; namely, the bulk Seebeck coefficients of the barrier and well materials S B and S W , respectively, and the energy relaxation length λ E . The difference between Eq. (3) and Eq. (16) is the presence of energy relaxation in the latter, which is a result of electron-optical phonon scattering. In fact, the energy relaxation process leads to the partial extension of S B into the well region up to a distance ≈ λ E from the barrier. We remark that in a long channel we can neglect the two wells close to the contacts to a good approximation, in which case Eq. (16) can be expressed as

$$S _ { s y s } \simeq \frac { ( L _ { \text{B} } S _ { \text{B} } / \kappa _ { \text{B} } ) + ( L _ { \text{W} } S _ { \text{W-relax} } / \kappa _ { \text{W} } ) } { ( L _ { \text{B} } / \kappa _ { \text{B} } ) + ( L _ { \text{W} } / \kappa _ { \text{W} } ) }. \quad ( 1 7 ) \quad \text{enter} _ { \text{$\text{sol}$} }.$$

Note that the presence of barriers causes increase of the average energy of the current and a consequent increase in the Seebeck coefficient, which is due to the energy filtering provided by the barriers [19]. On the other hand, the energy relaxation process of charge carriers in the well region causes electrons to propagate at lower energy states, which leads to reduction of the Seebeck coefficient. However, at the edge of a well and immediately after the barrier the Seebeck coefficient remains close to its highest value, i.e., close to S B . The conductivity mean-free-path of the carriers of these spatial regions can still be high as carriers propagate at higher velocity states compared to the relaxed well states, resulting in an increase of the PF for suitable well sizes [17, 31, 32]. This increase originates from these energy non-relaxed regions, which emphasizes the need for them to be captured accurately.

## III. VALIDATION OF MODEL WITH NEGF

We compare now the analytical model for the Seebeck coefficient with simulation results from 1D NEGF. However, we re-iterate that NEGF simulation plays the role here of validation, and the reduced dimensionality considered is the result of computational necessity on the part of NEGF simulation, but the analytical model itself is expected to be dimensionally agnostic. We consider a channel of length L = 250nm, an initial number of twenty four rectangular barriers ( n = 24) with spacing d = 4nm between them, and each one of thickness b = 5nm and height V B = 0 05eV. For simplicity . we consider equal thermal conductivities for the barrier and well materials, i.e., κ B = κ W . For the NEGF simulations we assume the same parameter values and a channel with arbitrarily small width W = 3nm to help with convergence which can be difficult in a truly 1D structure with such intense optical phonon scattering. The small width of the channel gives rise to an upward shift of the subband energies by an amount 0 025eV, resulting in an effective barrier . height of V B = 0 075eV. The channel with the effective . barrier height is as shown in Fig. 3(a).

## A. S vs d

In the following we use an e-ph coupling strength D 0 = 0 0026eV . 2 and optical phonon energy of glyph[planckover2pi1] ω = 0 06eV . (which is close to the Si value), while we place the Fermi level at E F = 0 05eV. This corresponds to degenerately . doped channel, where high PF s were observed [17]. The relaxation length is extracted by fitting Eq. (A5) [solid red line in Fig. 2(a)] on the NEGF simulation for 〈 E x ( ) 〉 yielding λ E = 16 5nm. . For these parameter values we also find 〈 E 〉 B = 0 112eV and . 〈 E 〉 W = 0 075eV, as de-. scribed in the context of Fig. 2. Using Eq. (10) and a similar one for S W we find S B = 2 06 . × 10 -4 V/K and S W = 0 833 . × 10 -4 V/K. In Fig. 3(a) we show the average energy of the current flow 〈 E x ( ) 〉 calculated from NEGF (solid red line) and that from the model, i.e., Eqs. (A1) and (A3) (dashed blue line) in the case of five barriers for which d = 50nm and x 0 = 12 5nm. . We notice that 〈 E x ( ) 〉 plotted from the model describes very well the simulation results and captures all essential features of the current flow, including the relaxation process. We remark here that the decay rate of 〈 E x ( ) 〉 within each well depends on the size of the well compared to the energy relaxation length, as discussed below [see also Eq. (11)].

In Fig. 3(b) we show the Seebeck coefficient S NEGF calculated from NEGF (solid, red line) as a function of well size d , and the Seebeck coefficient S sys calculated from the model with energy relaxation Eq. (16) (dashed, blue line). The dashed-dotted magenta line shows the Seebeck coefficient S no rel sys from the model without energy relaxation Eq. (3). S NEGF is calculated by integrating the average energy of the current flow with respect to the Fermi level when a voltage difference ∆ V is applied

FIG. 3: (Colour online) (a) Average energy of the current flow 〈 E x ( ) 〉 along the channel with five barriers calculated from Eq. (5) using NEGF (solid, red line) and from the analytical result Eq. (A3) (dashed, blue line) for which d = 50nm and λ E = 16 5nm.The color map indicates the current flow . I ( E,x ), with yellow indicating high, and green low current density. (b) Seebeck coefficient vs well size d calculated from (i) NEGF (solid, red line), (ii) model with energy relaxation (dashed, blue line), and (iii) model without energy relaxation (dashed-dotted, magenta line).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-cfefdafa-b987-4be7-9813-20f3de1ebb38.jpg)

at the channel contacts, as shown in Eqs. (4) - (6) (see also Ref. [33]). The bulk values of the Seebeck coefficients S B and S W are also shown for reference (dashed-dotted, black likes). In these calculations d increases by removing barriers sequentially one at a time while keeping L fixed. Notice that in our model we take into account the finite thickness of the barriers via Eq. (A1). Notice that S sys decreases with increasing d as a consequence of increasing energy relaxation in the well regions, and agrees very well with S NEGF . As d increases, 〈 E x ( ) 〉 gradually relaxes more in the well regions and it approaches 〈 E 〉 W in the middle of each well. However, d should be significantly larger than λ E in order to have full energy relaxation and to achieve the limits 〈 E x ( ) 〉 → 〈 E 〉 W and S sys → S W .

FIG. 4: (Colour online) (a) Average energy of the current flow 〈 E x ( ) 〉 along a channel with five barriers calculated from Eq. (5) using NEGF for increasing values of λ E (i.e., decreasing values of the deformation potential D 0 ). The color map indicates the current flow I ( E,x ), with yellow indicating high, and green low current density. (b) Seebeck coefficient vs D 0 calculated with NEGF (solid, red line) and with the model (dashed, blue line). Inset: the same but vs λ E .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-1ec0a370-80ba-4d82-afa1-b53b489e2255.jpg)

## B. S vs D 0

We compare now the model with the NEGF simulations in the case in which the e-ph coupling strength D 0 is varied. Here, again glyph[planckover2pi1] ω = 0 06eV. In order to illus-. trate our results we use a channel with five SL barriers. Since the energy relaxation length λ E decreases as D 0 increases, for each value of D 0 separately we made fitting of Eq. (A5) to the NEGF 〈 E x ( ) 〉 as in Fig. 2(a) and extracted the corresponding values of λ E . The values of 〈 E 〉 B and 〈 E 〉 W are determined as before. However, as D 0 varies, 〈 E 〉 W does not remain constant and varies slightly. Accordingly, S W is determined separately for each value of D 0 .

In Fig. 4(a) we show the average energy of the current flow 〈 E x ( ) 〉 calculated with NEGF for increasing values of λ E . Notice the gradually faster decay rate of 〈 E x ( ) 〉 within each well as λ E decreases (i.e., as D 0 increases). The physical origin of this behavior is the enhanced scat-

FIG. 5: (Colour online) (a) Average energy of the current flow 〈 E x ( ) 〉 along a channel with five barriers calculated from Eq. (5) using NEGF for increasing values of λ E (i.e., increasing values of the optical phonon energy glyph[planckover2pi1] ω ). The color map indicates the current flow I ( E,x ), with yellow indicating high, and green low current density. (b) Seebeck coefficient vs glyph[planckover2pi1] ω calculated with NEGF (solid, red line) and with the model (dashed, blue line). Inset: the same but vs λ E .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-0715f8cd-32b8-4bb2-9d93-df748ae65c7f.jpg)

tering of electrons with optical phonons as D 0 increases, resulting in gradually stronger carrier energy relaxation. Consequently the Seebeck coefficient decreases with increasing D 0 , as shown in Fig. 4(b). The solid red line is the NEGF simulation results while the dashed blue line is the result of the model. It can be seen that the agreement between the two results is very good. The inset shows the Seebeck coefficient vs the values of the energy relaxation length λ E that correspond to the values of D 0 that were used.

## C. S vs glyph[planckover2pi1] ω

We illustrate now the case in which the optical phonon energy glyph[planckover2pi1] ω is varied. Here, we fix D 0 = 0 0026eV . 2 despite the fact that the e-ph coupling strength is ∼ 1 /ω [27], because we intent to investigate the effects of the phonon energy alone, independent of the e-ph coupling. Again we use a channel with five SL barriers. In order to deter-

FIG. 6: (Colour online) Optical phonon energies glyph[planckover2pi1] ω for common thermoelectric materials.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-bf3933bb-ba75-46d1-91a1-e910e69b6a82.jpg)

mine the values of λ E we made fitting of Eq. (A5) to the NEGF 〈 E x ( ) 〉 separately for each value of glyph[planckover2pi1] ω . Also, for each value of glyph[planckover2pi1] ω we find the corresponding value of 〈 E 〉 W as described in the context of Fig. 2(b). In Fig. 5(a) we show the average energy of the current flow 〈 E x ( ) 〉 calculated with NEGF for increasing values of λ E . We notice that for longer energy relaxation lengths, i.e., λ E = 16nm (where glyph[planckover2pi1] ω &gt; V B -E F ) charge carriers, which travel at energies ≈ k B T above the barrier height, cannot now easily emit phonons as the final scattering states reside below E F , and they are almost filled, thus the relaxation rate is lower [see black line in Fig. 5(a)]. In fact, we have performed simulations with even higher glyph[planckover2pi1] ω , and found that the relaxation rate is suppressed even more. This could be a generic filtering design direction to suppress relaxation in nanostructured materials, by choosing V B -E F smaller compared to the material's glyph[planckover2pi1] ω . Thus, the average energy of the current flow increases and as a consequence the Seebeck coefficient also increases. This is shown in Fig. 5(b) where it can also be seen that S sys agrees very well with S NEGF to an accuracy of 1-2%. We also note that for glyph[planckover2pi1] ω = 0 02eV the Seebeck coefficient increases . slightly. In fact, as glyph[planckover2pi1] ω becomes even smaller the energy relaxation gradually diminishes and, in the limit glyph[planckover2pi1] ω → 0, the energy 〈 E x ( ) 〉 should become constant reflecting the absence of optical phonon processes (i.e., resembling the elastic acoustic phonon case or even the ballistic one). The inset shows the Seebeck coefficient vs the values of λ E that correspond to the values of glyph[planckover2pi1] ω on the x axis.

The values of glyph[planckover2pi1] ω for which we plotted the Seebeck coefficient correspond to those of the most common TE materials and semiconductors that are being explored experimentally. The optical phonon energies for some of these materials is shown in Fig. 6, where it is seen that glyph[planckover2pi1] ω ranges from ≈ 13meV for PbTe up to ≈ 63meV for Si. In passing, we remark that the slightly lower values of S sys throughout the range of glyph[planckover2pi1] ω compared to S NEGF in Fig. 5(b) are probably due to the fact that 〈 E 〉 B in the model is constant on the top of each barrier and slightly

FIG. 7: (Colour online) (a) Average energy of the current flow 〈 E x ( ) 〉 along the channel with five barriers calculated from Eq. (5) using NEGF (solid, red line) and from the analytical result Eq. (A3) (dashed, blue line) for which d = 50nm and λ E = 17 5nm. Here . E F = V B = 0 075eV. The color map . indicates the current flow I ( E,x ), with yellow indicating high, and green low current density. (b) Seebeck coefficient vs well size d calculated from (i) NEGF (solid, red line), (ii) model with energy relaxation (dashed, blue line), and (iii) model without energy relaxation (dashed-dotted, magenta line).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-83a7ce01-a33c-4a69-91b4-a6f478af008c.jpg)

lower than 〈 E x ( ) 〉 of NEGF.

## D. S vs d ( E F = V B )

In the above analysis the Fermi level was taken at E F = 0 05eV, i.e., . ≈ k B T below the barrier height. We explore now the regime for which the Fermi level is aligned with the barrier height V B , i.e., E F = V B = 0 075eV, while the . rest of the parameters are the same as previously. This case is interesting because it has been shown previously [18] that this is the optimal case for PF improvement if relaxation is suppressed.

The relaxation length is again extracted graphically by fitting Eq. (A5) on the NEGF simulation result, which now yields λ E = 17 5nm. Also, in the same manner as in . the previous case, we find 〈 E 〉 B = 0 117eV and . 〈 E 〉 W = 0 096eV, which yield . S B = 1 4 . × 10 -4 V/K and S W =

FIG. 8: (Colour online) (a) Average energy of the current flow 〈 E x ( ) 〉 along the channel with five barriers calculated from Eq. (5) using NEGF (solid, red line) and from the analytical result Eq. (A3) (dashed, blue line) for which d = 50nm and λ E = 17 5nm. Here . E F = V B = 0 075eV. The color map . indicates the current flow I ( E,x ), with yellow indicating high, and green low current density. (b) Seebeck coefficient vs grain size d calculated from (i) NEGF (solid, red line), (ii) model with energy relaxation (dashed, blue line), and (iii) model without energy relaxation (dashed-dotted, magenta line).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-528dca0f-7320-4106-b774-cce7fd3029c8.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-ce992780-577c-43dd-bac7-53afd6722865.jpg)

0 7 . × 10 -4 V/K. In Fig. 7(a) we show the average energy of the current flow 〈 E x ( ) 〉 calculated from NEGF (solid, red line) and that from the model, i.e., Eqs. (A1) and (A3) (dashed, blue line) in the case of five barriers. The rate of energy relaxation in this case is significantly smaller than in the previous case of Sec. IIIA where the Fermi level was at E F = V B -k B T [compare with Fig. 3(a)]. Indeed, 〈 E 〉 B - 〈 E 〉 W = 0 023eV when . E F = V B , while 〈 E 〉 B - 〈 E 〉 W = 0 039eV when . E F = V B -k B T , i.e., more than 40% decrease. The reason for the suppressed relaxation in this case where V B = E F compared to the previous one where V B &gt; E F , is simply because electrons in the wells now tend to relax at the Fermi level (which is at the barrier level) and not below the barrier level (although the Seebeck coefficients are lower now due to the higher E F . The variation of the Seebeck coefficient with glyph[planckover2pi1] ω is in this case smaller as well, since the current flow is closer to E F , and out-scattering to filled lower

energies away from E F is more difficult.

In Fig. 7(b) we show the Seebeck coefficient S NEGF calculated from NEGF (solid, red line) as a function of well size d , and the Seebeck coefficient S sys calculated from the model with energy relaxation Eq. (16) (dashed, blue line). The dashed-dotted magenta line shows the Seebeck coefficient S no rel sys from the model without energy relaxation Eq. (3). The bulk values of the Seebeck coefficients S B and S W are also shown for reference (dashed-dotted, black likes). Again, as in the previous case of Fig. 3, in these calculations d increases by removing barriers sequentially one at a time. We note the very good agreement of S sys with S NEGF to an accuracy of up to ≈ 5% for large d . We also note that, as a consequence of the slower rate of energy relaxation, the Seebeck coefficient also decreases at a slower rate with increasing d than in the previous case illustrated in Fig. 3(b).

## E. S vs glyph[planckover2pi1] ω ( E F = V B )

We vary now the optical phonon energy glyph[planckover2pi1] ω , while keeping the deformation potential fixed at D 0 = 0 0026eV . . 2 We expect that the average energy of the current flow and the Seebeck coefficient exhibit the same behaviour as in the previous case. In Fig. 8(a) we show the average energy of the current flow 〈 E x ( ) 〉 calculated with NEGF for increasing values of λ E . These values were extracted from fitting of Eq. (A5) to the NEGF result for each value of glyph[planckover2pi1] ω as we have done in Sec. IIIC and, in addition, each value of 〈 E 〉 W was determined as in Sec. IIIC. We notice the small effect of glyph[planckover2pi1] ω , which is also reflected in the Seebeck coefficient. This is shown in Fig. 8(b) where it is also seen that S sys agrees very well with S NEGF to an accuracy of 1 -2% for smaller values of glyph[planckover2pi1] ω . We also notice that the Seebeck coefficient exhibits identical behaviour as that shown in Sec. IIIC [see Fig. 5(b)] except the small variation with glyph[planckover2pi1] ω .

## IV. COMPARISON TO EXPERIMENT

To partially validate the Seebeck model with energy relaxation we constructed, we use the measured data in the experiment of Ref. [15, 16], for the case of the Seebeck coefficients in SLs based on ErAs doped InGaAs wells and InGaAlAs barriers. The papers provide the estimated band offsets of the wells and barriers compared to the position of the Fermi level, as well as the measured Seebeck coefficients for the in-plane and cross-plane directions. Although we do not have access to other necessary parameters to compute electronic transport reliably in correlation with the experiment, at first order we can still approximate the Seebeck coefficient using Boltzmann Transport theory under the relaxation time approximation as:

$$S = \frac { q k _ { \text{B} } } { \sigma } \int _ { E _ { 0 } } ^ { \infty } d E \left ( - \frac { \partial f _ { 0 } } { \partial E } \right ) \Xi ( E ) \left ( \frac { E - E _ { \text{F} } } { k _ { \text{B} } T } \right ), \quad ( 1 8 ) \quad \text{to be} \quad \\ \text{comp}$$

FIG. 9: (Colour online) Comparison between experimental and theoretical Seebeck coefficients for ErAs:InGaAs/InGaAlAs SLs with different doping concentrations from Ref. [15]. The black-solid and black-dashed lines show measured data for cross-plane and in-plane Seebeck coefficients in the SLs. The red-solid and red-dashed lines show the theoretically estimated Seebeck coefficients of the uniform well S W and barrier S B , respectively. The bluesolid and blue-dahed lines show the theoretical calculations for the model that accounts for the relaxations physics, and the model that considers independent regions, respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01090v1-figure-46325c93-78c6-47c2-a926-a5f4393dd310.jpg)

where

$$\sigma = q ^ { 2 } \int _ { E _ { 0 } } ^ { \infty } d E \left ( - \frac { \partial f _ { 0 } } { \partial E } \right ) \Xi ( E ), \text{ \quad \ \ } ( 1 9 )$$

with the transport distribution function Ξ( E ) defined as:

$$\Xi ( E ) = \tau ( E ) \upsilon ( E ) ^ { 2 } g ( E ), \text{ \quad \ \ } ( 2 0 )$$

where υ E ( ) is the carriers' velocity, g E ( ) is the density of states, and τ ( E ) is the relaxation time. Notice that υ E ( ) 2 ∼ E , while in 3D g E ( ) ∼ E 1 / 2 , and it is common to express the energy dependence of the relaxation times for acoustic phonons and ionized impurity scattering as ∼ E -1 / 2 and ∼ E 3 / 2 , respectively. In the experiment a series of doping values in the 10 18 - 10 19 cm -3 range were used, thus, we employ a mixed scattering relaxation time exponent in BTE as r = 1 2 [34]. / Thus, the Seebeck coefficients for the well S W and the barrier S B can be approximated using the band edges provided in the experimental paper as:

$$S = \frac { k _ { \text{B} } } { q } \frac { \int _ { E _ { 0 } } ^ { \infty } E ^ { 2 } \left ( - \frac { \partial f _ { 0 } } { \partial E } \right ) \left ( \frac { E - E _ { F } } { k _ { B } T } \right ) d E } { \int _ { E _ { 0 } } ^ { \infty } E ^ { 2 } \left ( - \frac { \partial f _ { 0 } } { \partial E } \right ) d E }, \quad ( 2 1 )$$

where now at first order it does not depend on material parameters. In this structure the barrier thickness is b = 10nm, while the well length is d = 20nm, and the only parameter needed is the relaxation length, which is taken to be λ E = 30nm, to reflect the higher mobility of InGaAs compared to Si.

Figure 9 shows the measured data for the in-plane and cross plane Seebeck coefficients by the black-dashed and black solid lines, respectively versus carrier density. The red-dashed and red-solid lines show the uniform channel calculated S W and S B (upper and lower limits of our calculations). The blue-dashed line shows the calculated Seebeck coefficients in the case where each region in considered to be independent (no relaxation physics considered), whereas the blue-solid line when the relaxation physics is considered by the model developed. Despite the large uncertainties of this evaluation, the model (blue-solid line) is in the range of the measured crossplane data (black-solid line). We find that slightly better fit can be obtained by adjusting the scattering time exponents, however, we do not attempt better fit as any exponent we use will be purely speculative. On the other hand, when the barrier and well are considered to be independent, the independent region model resides lower and coincides with the in-plane measured Seebeck data (although this could just be accidental). Despite the large uncertainties, this analysis shows the validity of the developed model in describing the Seebeck coefficient of SLs and nanocomposite systems.

## V. SUMMARY AND CONCLUSIONS

In summary, we presented a simple analytical model for the Seebeck coefficient S of superlattice materials (or nanocomposites to first order approcimation) in the presence of energy relaxation due to electron-optical phonon scattering. This model casts the complex and crucial physics of semi-relaxation and its role in Seebeck enhancement in terms of only three material parameters: the bulk Seebeck coefficients of the constituent materials and the energy relaxation length λ E , which is related to the energy relaxation time τ E , of the more conductive 'well' material. Thus, it is our hope that the model can help guide future nano-engineering efforts aimed at Seebeck coefficient and power factor enhancement.

To validate this model, numerical simulations were performed using the fully quantum mechanical nonequilibrium Green's function method and very good agreement was found. We also compared the model with experimental values for ErAs:InGaAs/InGaAlAs superlattice, and good agreement was found. The variation of the calculated Seebeck coefficient with increasing well size d and Fermi level E F , as well as with increasing e-ph coupling strength D 0 , and with increasing optical phonon energy glyph[planckover2pi1] ω , the physical parameters that determine λ E , were also studied. We also provided an expression for the average energy of the current flow 〈 E x ( ) 〉 , which agrees very well with the NEGF result and captures accurately the behavior of 〈 E x ( ) 〉 in the well regions. We expect that these results will be helpful and useful to experimentalists in their determination of the Seebeck coefficients of nanocomposite structures and superlattices. Note that the paper deals exclusively with the Seebeck coefficient because it can be trivially mapped to the average energy of the current flow. Similar considerations for the electrical conductance in the semi-energy relaxing regions between the barriers and the wells are more complicated, as there is no direct map to the average energy of the current flow.

Finally, we note that in all simulations of this work and in the construction of the model, we considered periodic superlattice structures. However, we argue that the model is at first order applicable to nanocomposite/nanocrystalline materials as well. Nanocomposites are described by a 3D aperiodic geometry, and strictly speaking the complexity of the transport paths is such that would not allow us to map the 3D onto 1D paths beyond a first order estimation. Superlattice geometries thought, can be considered as a limiting case for a nanocomposite system, becoming more accurate as the variance of barrier spacing and size decreases and the structure becomes more uniform in shape and distribution. Indeed, in previous works of ours [9, 35] we pointed out that in the presence of statistical variability in the sizes of the domains, the overall Seebeck coefficient is rather robust. Thus, we still believe that the model developed provides a first order estimate to the Seebeck coefficient of 3D nanocomposites and nanocrystalline structures as well.

Acknowledgments This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 678763). We also thank Samuel Foster and Dhritiman Chakraborty for helpful discussions.

## Appendix A: Average energy of the current flow 〈 E x ( ) 〉 for SL structures

We consider a 1D SL structure of length L , as shown in Fig. 1(a). We make the approximation of constant 〈 E x ( ) 〉 over each barrier, i.e.,

$$\text{ment} \\ \text{and} \\ \text{alc} \text{in} -$$

where 〈 E 〉 B is the average energy on the barrier, Θ( ... ) is the Heaviside function, n is the total number of barriers, and t = b + d . However, within a well, 〈 E x ( ) 〉 decays exponentially [30]. For simplicity, we consider one well with length d , surrounded by two barriers at x = 0 and x = d . In addition to the exponential decay, 〈 E x ( ) 〉 should satisfy two boundary conditions, i.e., 〈 E x ( ) 〉 = 〈 E 〉 B at x = 0 and x = d . The expression for 〈 E x ( ) 〉 should also satisfy the equilibrium condition for large d , i.e., 〈 E x ( ) 〉 = 〈 E 〉 W at x = d/ 2 as d →∞ , where 〈 E 〉 W is the average energy in the well region under equilibrium. This means that as the well size becomes large (i.e., d glyph[greatermuch]

λ E ), the energy of the charge carriers is fully relaxed within the well and 〈 E x ( ) 〉 reaches the bulk limit. The differential equation that yields solutions which satisfy these conditions is of second order and given as

$$\frac { d ^ { 2 } \langle E \rangle } { d x ^ { 2 } } - \frac { \langle E \rangle } { \lambda _ { E } ^ { 2 } } = \frac { 1 } { \lambda _ { E } ^ { 2 } } \left ( \langle E \rangle _ { B } - \langle E \rangle _ { W } \right ) e ^ { - d / \lambda _ { E } } - \frac { \langle E \rangle _ { W } } { \lambda _ { E } ^ { 2 } }. \quad \text{differ} \, \\ \intertext { T h e s s t i n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e n d e }$$

The solution of Eq. (A2) proceeds in a straightforward manner and is given as

$$\langle E ( x ) \rangle = & \sum _ { \ell = 0 } ^ { n - 2 } ( \langle E \rangle _ { B } - \langle E \rangle _ { \text{W} } ) [ e ^ { - ( x - x _ { 0 } - b - \ell t ) / \lambda _ { E } } \\ & + e ^ { ( x - x _ { 0 } - ( \ell + 1 ) t ) / \lambda _ { E } } - e ^ { - d / \lambda _ { E } } ] + \langle E \rangle _ { \text{W} }. \ ( A 3 )$$

In Eq. (A3) we notice that as d becomes large, 〈 E x ( ) 〉 gradually relaxes, and 〈 E d/ ( 2) 〉 → 〈 E 〉 W as d →∞ .

In the wells at the channel ends, i.e., close to the contacts, 〈 E x ( ) 〉 satisfies different boundary conditions. For

- [1] G. Ding, G. Y. Gao, Z. Huang, W. Zhang and K. Yao, Nanotechnology 27 , 375703 (2016).
- [2] H. Huang, Y. Cui, Q. Li, C. Dun, W. Zhou, W. Huang, L. Chen, C. A. Hewitt, and D. L. Carroll, Nano Energy 26 , 172 (2016).
- [3] W. Huang, X. Luo, C. K. Gan, S. Y. Quek, and G. Liang, Phys. Chem. Chem. Phys. 16 , 10866 (2014).
- [4] M. Beekman, D. T. Morelli, and G. S. Nolas, Nat. Mater. 14 , 1182 (2015).
- [5] C. Fu, S. Bai, Y. Liu, Y. Tang, L. Chen, X. Zhao, and T. Zhu, Nat. Commun. 6 , 8144 (2015).
- [6] L.-D. Zhao, S.-H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. P. Dravid, and M. G. Kanatzidis, Nature (London) 508 , 373 (201).
- [7] D. Beretta, N. Neophytou, J. M. Hodges, M. G. Kanatzidis, D. Narducci, M. M.Gonzalez, M. Beekman, B. Balke, G. Cerrettii, W. Tremel et al., 'Thermoelectrics: From history, a window to the future,' Mater. Sci. Eng. R Rep. (to be published).
- [8] H. Mizuno, S. Mossa, and J.-L. Barrat, Sci.Rep. 5 , 14116 (2015).
- [9] M. Thesberg, M. Pourfath, H. Kosina, and N. Neophytou, J. Appl. Phys. 118 , 224301 (2015).
- [10] M. Thesberg, H. Kosina, and N. Neophytou, J. Appl. Phys. 120 , 234302 (2016).
- [11] Y. M. Zhou and L.-D. Zhao, Adv. Mater. 29 , 1702676 (2017).
- [12] P. Priyadarshi, A. Sharma, S. Mukherjee, and B. Muralidharan, J. Phys. D: Appl. Phys. 51 , 185301 (2018).
- [13] R. Venkatasubramanian, T. Colpitts, E. Watko, and J. Hutchby, in Proceedings IEEE 15th International Conference on Thermoelectrics (IEEE Service Centre, 1996).
- [14] A. Ishida, T. Yamada, D. Cao, Y. Inoue, M. Veis, and T. Kita, J. Appl. Phys. 106 , 023718 (2009).
- [15] G. Zeng, J. M. O. Zide, W. Kim, J. E. Bowers, A. C. Gossard, Z. Bian, Y. Zhang, A. Shakouri, S. L. Singer, and A. Majumdar, J. Appl. Phys. 101 , 034502 (2007).
- [16] Z. Bian, M. Zebarjadi, R. Singh, Y. Ezzahri, A. Shakouri,

simplicity, we consider the well at the right channel end and we make the transformation x → -x x 0 -( n -1) + . t b Then, 〈 E x ( ) 〉 should only satisfy the boundary conditions 〈 E x ( ) 〉 = 〈 E 〉 B at x = 0 and, in addition, the equilibrium condition 〈 E x ( ) 〉 = 〈 E 〉 W at x = d / ′ 2 as d ′ → ∞ . The differential equation that yields solutions which satisfy these conditions is of first order and given as

$$\frac { d \langle E \rangle } { d x } + \frac { \langle E \rangle } { \lambda _ { \text{E} } } = \frac { \langle E \rangle _ { \text{W} } } { \lambda _ { \text{E} } }. \quad \quad \quad ( \text{A4} )$$

The presence of the term on the right hand side of Eq. (A4) guarantees that the solution satisfies the equilibrium condition. The solution of Eq. (A4) is given as

$$\langle E ( x ) \rangle = \left ( \langle E \rangle _ { \mathrm B } - \langle E \rangle _ { \mathrm W } \right ) e ^ { - x / \lambda _ { \mathrm E } } + \langle E \rangle _ { \mathrm W }. \quad ( A 5 )$$

The solution Eq. (A5) is relevant to the well in the right channel end. For the well in the left channel end, we just make the replacement x → d ′ -x in Eq. (A5).

- G. Zeng, J.-H. Bahk, J. E. Bowers, J. M. O. Zide, and A. C. Gossard, Phys. Rev. B 76 , 205311 (2007).
- [17] N. Neophytou, X. Zianni, H. Kosina, S. Frabboni, B. Lorenzi, and D. Narducci, Nanotechnology 24 , 205402 (2013).
- [18] V. Vargiamidis and N. Neophytou, Phys. Rev. B 99 , 045405 (2019).
- [19] N. Neophytou and H. Kosina, J. Appl. Phys. 114 , 044315 (2013).
- [20] A. Popescu, L. Woods, J. Martin, and G. Nolas, Phys. Rev. B 79 , 205302 (2009).
- [21] C. Bera, M. Soulier, C. Navone, G. Roux, J. Simon, S. Volt, and N. Mingo, J. Appl. Phys. 108 , 124306 (2010).
- [22] M. Bartkowiak, G. D. Mahan, and M. T. Terry, in Semiconductors and Semimetals (Elsevier, 2001), Vol. 70, pp. 245-271.
- [23] S. O. Koswatta, S. Hasan, M. S. Lundstrom, M. P. Anantram, and D. E. Nikonov, IEEE Trans. Electron Devices 54 , 2339 (2007).
- [24] M. P. Anantram, M. Lundstrom, and D. Nikonov, Proc. IEEE 96 , 1511 (2008).
- [25] R. Kim and M. Lundstrom, J. Appl. Phys. 110 , 034511 (2011).
- [26] X. Zianni and D. Narducci, J. Appl. Phys. 117 , 035102 (2015).
- [27] M. Lundstrom, Fundamentals of Carrier Transport , Cambridge, UK: Cambridge Univ. Press, 2000.
- [28] B. Gonzalez, V. Palankovski, H. Kosina, A. Hernandez, and S. Selberherr, Solid-State Electron. 43 , 1791 (1999).
- [29] R. Kim and M. Lundstrom, J. Appl. Phys. 111 , 024508 (2012).
- [30] B. Moyzhes and V. Nemchinsky, Appl. Phys. Lett. 73 , 1895 (1998).
- [31] N. Neophytou, X. Zianni, H. Kosina, S. Frabboni, B. Lorenzi, and D. Narducci, J. Elec. Mater. 43 , 1896 (2014).
- [32] M. Thesberg, H. Kosina, and N. Neophytou, Phys. Rev. B 95 , 125206 (2017).

glyph[negationslash]

[33] In experimental settings, one extracts the Seebeck coefficient from the open circuit voltage upon the application of a thermal gradient along the channel, as S = ∆ V/ ∆ , which equivalently T can also be computed by S = I (∆ V =0) /G ∆ . In [25], it was validated that the two T methods of extracting the Seebeck coefficient are equivalent, which makes it easier in time consuming simulations (as the ones we undertake) to only run the ∆ V = 0 case and still be able to extract the Seebeck coefficient by in-

- tegrating the energy of the current flow over the length of the channel.
- [34] J. Mao, J. Shuai, S. Song, Y. Wu, R. Dally, J. Zhou, Z. Liu, J. Sun, Q. Zhang, C. dela Cruz, S. Wilson, Y. Pei, D. J. Singh, G. Chen, C.-W. Chu, and Z. Ren, PNAS 114 , 40, 10548 (2017).
- [35] M. Thesberg, M. Pourfath, N. Neophytou and H. Kosina, J. Electron. Mater. 45 , 1584 (2016).
Vassilios Vargiamidis Mischa Thesberg Neophytos Neophytou
Theoretical model for the Seebeck coefficient in superlattice materials with energy relaxation
We present an analytical model for the Seebeck coefficient S of superlattice
materials that explicitly takes into account the energy relaxation due to
electron-optical phonon (e-ph) scattering. In such materials, the Seebeck
coefficient is not only determined by the bulk Seebeck values of the materials
but, in addition, is dependent on the energy relaxation process of charge
carriers as they propagate from the less-conductive barrier region into the
more-conductive well region. We calculate S as a function of the well size d,
where carrier energy becomes increasingly relaxed within the well for d greater
than l, where l is the energy relaxation length. We validate the model against
more advanced quantum transport simulations based on the nonequilibrium Green
function (NEGF) method and also with an experiment, and we find very good
agreement. In the case in which no energy relaxation is taken into account, the
results deviate substantially from the NEGF results. The model also yields
accurate results with only a small deviation (up to ~3%) when varying the
optical phonon energy hw or the e-ph coupling strength D0, physical parameters
that would determine l. As a first order approximation, the model is valid for
nanocomposite materials, and it could prove useful in the identification of
material combinations and in the estimation of ideal sizes in the design of
nanoengineered thermoelectric materials with enhanced power factor performance.