Computational study of the Seebeck coefficient of one-dimensional composite nanostructures

The Seebeck coefficient (S) of composite nano-structures is theoretically explored within a self-consistent electro-thermal transport simulation framework using the non-equilibrium Green’s function method and a heat diffusion equation. Seebeck coefficients are determined using numerical techniques that mimic experimental measurements. Simulation results show that, without energy relaxing scattering, the overall S of a composite structure is determined by the highest barrier within the device. For a diffusive, composite structure with energy relaxation due to electron-phonon scattering, however, the measured S is an average of the position-dependent values with the weighting factor being the lattice temperature gradient. The results stress the importance of self-consistent solutions of phonon heat transport and the resulting lattice temperature distribution in understanding the thermoelectric properties of a composite structure. It is also clarified that the measured S of a composite structure reflects its power generation performance rather than its cooling performance. The results suggest that the lattice thermal conductivity within the composite structure might be engineered to improve the power factor over the bulk by avoiding the conventional trade-off between S and the electrical conductivity.The Seebeck coefficient (S) of composite nano-structures is theoretically explored within a self-consistent electro-thermal transport simulation framework using the non-equilibrium Green’s function method and a heat diffusion equation. Seebeck coefficients are determined using numerical techniques that mimic experimental measurements. Simulation results show that, without energy relaxing scattering, the overall S of a composite structure is determined by the highest barrier within the device. For a diffusive, composite structure with energy relaxation due to electron-phonon scattering, however, the measured S is an average of the position-dependent values with the weighting factor being the lattice temperature gradient. The results stress the importance of self-consistent solutions of phonon heat transport and the resulting lattice temperature distribution in understanding the thermoelectric properties of a composite structure. It is also clarified that the measured S of a composite structure reflects it...


I. INTRODUCTION
When a temperature difference, DT, is applied across a material sample, an open-circuit voltage, V oc , can be induced, generating an electric field that opposes the temperature gradient.This is called the Seebeck effect, 1 and the proportionality constant is the Seebeck coefficient (or thermopower) S, which gives V oc ¼ ÀSDT.The Seebeck coefficient can be looked on as the entropy transport per charged particle, 2 and S > 0 for hole conduction and S < 0 for electron conduction.The Seebeck effect is central to thermoelectric (TE) operation, and S is a key parameter in the thermoelectric figure of merit ZT ¼ S 2 rT=j that represents the efficiency of thermoelectric energy conversion, 1 where r is the electrical conductivity, T is the absolute temperature, and j is the thermal conductivity.For a homogeneous material, S can be calculated using the Boltzmann transport equation (BTE) 3,4 or Landauer formalism 5,6 in an integral form as where q is the elementary charge (q ¼ -e for electrons and q ¼ e for holes), E is the energy, E F is the Fermi level, and f 0 is the Fermi-Dirac distribution.6][7] It includes the band structure information and the effect of carrier scattering.][11] For realistic devices, however, Eq. ( 1) may not be used because the device structure may be inhomogeneous.For example, nanowire (NW) devices usually have metal contacts at the ends, [12][13][14] which introduce Schottky barriers.The potential barriers at the channel ends produce a non-uniform potential profile along the channel, which may affect the overall S of the NW device, making it different than that of a homogeneous NW calculated from Eq. (1).Understanding the Seebeck coefficient of a non-uniform, composite structure is becoming more relevant as the nano-engineered structures, such as superlattices 15 and nano-composites, 16 are attracting much attention as a promising way to further improve ZT.
In this paper, we explore the Seebeck coefficient of composite nano-structures within a self-consistent quantum transport simulation framework.The key questions to be addressed are: 1) What determines the overall S of a composite nano-structure?2) What roles do the length scales (e.g., energy relaxation length k E , the grain size d, etc.) play?3) How does the measured S relate to the cooling or power generation performance of a composite TE material?4) Is there a way to use the composite structure to modify the S versus r trade-off 1 and improve the power factor S 2 r?We use simple model structures and scattering mechanisms and restrict our attention to one-dimensional (1D) structures, but we expect that the general understanding established in this paper will be broadly applicable.
This paper is organized as follows: in Sec.II, we describe the simulation framework for the self-consistent electro-thermal transport in a 1D composite nano-structure.We also explain the techniques to numerically "measure" S a) Author to whom correspondence should be addressed.Electronic mail: kim369@purdue.edu. in our simulation.In Sec.III, simulation results are presented for ballistic and diffusive transport, and the results are related to a simple electrical-thermal circuit model.The role of k E on the measured S is also clarified.In Sec.IV, we discuss the meaning of the measured S regarding the TE performance.Possibilities to modify the S versus r trade-off and improve the power factor in a composite structure are also discussed.Conclusions follow in Sec.V. We assume a 1D NW with the doping densities varying along the transport (x) direction.In Fig. 1(a), L w is the length of the well region, L b is the barrier thickness, N w is the doping density in the well, and N b is the doping density in the barrier.The wire ends are connected to contact 1 and contact 2, which are ideal reservoirs maintained under equilibrium. 7o convert the three-dimensional (3D) doping densities to 1D and solve the Poisson equation, we assume a circular cross-section with a diameter D. In all following simulation results, we use N w ¼ 3 Â 10 19 cm À3 (heavily doped), N b ¼ 1 Â 10 16 cm À3 (essentially undoped), and D ¼ 3 nm.The model structure is intended to represent a potential barrier (grain boundary) between two doped regions (grains).In a realistic nano-composite material, the potential barrier comes from a charge at the grain boundary due to point defects, etc. 17 The simulation framework for the self-consistent electro-thermal transport is summarized in Fig. 1(b).The "electron part" treating the carrier transport and electrostatics is solved self-consistently with the "phonon part" that describes the lattice heat conduction.Within the electron part, the carrier transport is treated using the 1D non-equilibrium Green's function (NEGF) method considering three transport models, i.e., ballistic transport, elastic scattering, and inelastic scattering.Most of the previous theoretical studies on nano-composite materials adopt the BTE approach, 18,19 treating the effect of grain boundaries as another scattering mechanism with some relaxation time, s.As the grain size gets smaller and approaches the electron wavelength, 16,20 however, a quantum transport simulation framework, such as the NEGF, 21,22 is required to better understand the electron transport in nano-composites.As shown in Fig. 1 The simulation framework for the electron part described so far is widely used for conventional electronic device simulations, assuming a constant lattice temperature T L . 23In this work, we add another branch, i.e., the phonon part, to treat the power dissipation P from the electron to the phonon bath 21,22 and calculate T L by solving a heat transport equation, which again affects the electron part, as shown in Fig. 1(b).In our simulation, we solve a 1D heat diffusion equation, assuming some lattice thermal conductivity, j L .More details are discussed in the Appendix.][26] Next, we discuss how to determine the overall S. We first review the coupled current equation 1,5

II. APPROACH
where I is the electrical current, G is the electrical conductance, DV is the voltage difference, I q is the heat current, and K 0 is the electronic thermal conductance for zero DV.Alternatively, Eq. ( 2) can be expressed as where P is the Peltier coefficient and K e is the electronic thermal conductance for zero current.Note that the Kelvin relation gives P ¼ TS. Figure 2 shows the two possible configurations for the Seebeck coefficient measurement.First, S can be determined by measuring electrical currents, as shown in Fig. 2(a).From Eq. (2), I ¼ GDV for a finite and T 1 (T 2 ) are the voltage and temperature applied to contact 1 (contact 2), respectively.For DV ¼ 0 and a finite DT, we obtain I ¼ SGDT and then the Seebeck coefficient can be calculated from the ratio of the two coefficients as S ¼ SG=G.Another way to determine S is to measure the open circuit voltage DV for a DT, as shown in Fig. 2(b), and use the relation S ¼ ÀDV=DT from Eq. (3).3][14] In Sec.III, we use the two approaches to numerically "measure" the Seebeck coefficient of a composite nano-sturcture and compare the results.For the approach in

III. RESULTS
In this section, we present simulation results for ballistic transport, elastic scattering, and inelastic scattering and explore how carrier scattering affects the measured S of a composite nano-structure.Figure 3(a) shows the ballistic transport simulation results for the energy-and positionresolved electrical current and j L ¼ 150 W/m-K, which is the bulk Si value. 12The Fermi level of contact 1 (E F1 ) lies at 0 eV; E F1 ¼ E F ¼ 0 eV.There is no carrier scattering within the device, so I E; x ð Þ is uniform along the xdirection, and the average energy of the current flow E h i is constant.We numerically measure S using the two approaches described in Fig. 2, and the results are S ¼ -346 lV=K for the current measurement approach and S ¼ -348 lV=K for the voltage measurement approach.The two results are consistent and can be verified in the following way.From the constant E h i in Fig. 3(a), we can calculate the Peltier coefficient of the device as P ¼ E À E F h i=q ¼ -0.103 V, and the Kelvin relation gives S ¼ P=T ¼ -344 lV=K, which is consistent with the numerical measurement results.Also note that, for a 1D ballistic conductor, S can be calculated analytically as 5 where and E C is the conduction bandedge.For g F ¼ -3.04, which is extracted from the simulation result at the top of the barrier in Fig. 3(a), we obtain S ¼ -351 lV=K, which is again consistent with the numerically measured S values.
Figure 3(b) shows the simulation results for I E; x ð Þ for elastic scattering with a deformation potential of D 0 ¼ 0.01 eV 2 .(The meaning of D 0 within the NEGF framework 7 and its relation to the conventional expressions in the BTE approach 28 are discussed in the Appendix.)All other parameters in Fig. 3(b) are the same as those of Fig. 3(a).Still, the two approaches in Fig. 2 give consistent results as S ¼ -317 lV=K.Note that carrier scattering broadens energy levels, 7 which results in a smaller E h i than that of the ballistic case in Fig. 3(a).Although elastic scattering relaxes momentum and reduces G, it does not relax energy so that I E; x ð Þ and E h i are still uniform along the x-direction.Therefore, we can calculate the overall P of the device as P ¼ E À E F h i=q ¼ -0.0945 V, and the Kelvin relation gives S ¼ -315 lV=K, which is consistent with the numerical measurement results.
As discussed so far, when there is no energy relaxation within the device, S is determined by the potential barrier.Phase or momentum breaking scattering broadens the levels, effectively lowering the barrier height, as shown in Fig. 3(b), 7 but still the overall S is determined by the barrier region unless there is energy relaxing scattering.In these cases, P is uniform along the device, so we can define an overall P of the device, and the S obtained from the Kelvin relation is consistent with the numerically measured S.
ð Þand the average E of the current flow E h i are uniform along the x-direction.Using the approaches in Fig. 2, S is determined to be -346 lV=K (from current measurements) and -348 lV=K (from the voltage measurement).The Kelvin relation gives Next, we introduce energy relaxing scattering and see how the results change.Figure 4(a) shows the simulation result for I E; x ð Þ with DV ¼ 1 mV and T 1 ¼ T 2 ¼ 300 K for optical phonon scattering with D 0 ¼ 0.01 eV 2 and hx o ¼ 20 meV, where h is the reduced Planck constant and x o is the frequency of the optical phonon.All other parameters are the same as those in Fig. 3. Unlike the results in Fig. 3 without energy relaxation, I E; x ð Þ and E h i are non-uniform along the x-direction and we cannot define a constant P for the overall device.For example, the Peltier coefficient near contact 1 is P 1 ¼ -0.0702 V while it is higher in magnitude in the barrier region, where P b ¼ -0.0893 V in Fig. 4(a).The two approaches to measure S (recall Fig. 2) give consistent results as S ¼ -245 lV=K.Note that the measured S is neither P 1 =T ¼ -234 lV=K nor P b =T ¼ -298 lV=K, but a value somewhere in-between.
To understand the measured S of the composite nanostructure with energy relaxing scattering, we re-visit the coupled current equations in Eq. ( 2).Note that the T in Eq. ( 2) actually means the electron temperature T e as In the contacts, electrons are in equilibrium with the phonon bath, so T e1 ¼ T L1 ¼ T 1 and T e2 ¼ T L2 ¼ T 2 , where T e1 (T e2 ) and T L1 (T L2 ) are the electron and lattices temperatures for contact 1 (contact 2), respectively.When there is no energy relaxation (e.g., ballistic transport), T e is not well-defined within the device, 8 but the transport properties, such as S, are well-defined across the terminals.In Fig. 3, for example, carriers injected from the contacts with T e1 ¼ T 1 and T e2 ¼ T 2 only see the highest potential barrier and that determines the overall S of the device, regardless of the potential profile in the well region.When the carrier energy is relaxed by electron-phonon (e-ph) scattering, however, T e follows T L . 8As discussed in more detail in the Appendix, the lattice heat transport model gives T L x ð Þ, which determines the carrier energy relaxation rate due to e-ph interaction along the x-direction.This results in the carrier distribution with T e x ð Þ $ T L x ð Þ, and T e x ð Þ governs the coupled current equations in Eq. ( 5) within the device.
To understand the measured S of a diffusive composite structure in a simple way, we consider a model composite structure in Fig. 5. Two dissimilar regions with different S and G (S 1 and G 1 for region 1, S 2 and G 2 for region 2) are connected in series, and we measure the overall S using the two approaches in Fig. 2. Using a simple electrical-thermal circuit model (see Appendix for details) for the two configurations in Fig. 5, we obtain the same expression for the overall S as where DT L1 (DT L2 ) is the lattice temperature difference across region 1 (region 2) and DT L1 þ DT L2 ¼ DT.Equation (6) shows that the overall S of a diffusive composite structure is a weighted average of its component S values, and the weighting factor is the temperature difference across each component. 29,30More generally, Eq. ( 6) can be expressed in an integration form as where S x ð Þ is the x-dependent Seebeck coefficient.In Figs.

4(b)-4(c), we show simulation results for
for the model device in Fig. 4(a).From Eq. ( 7) we obtain S ¼ -243 lV=K, which is consistent with the numerical measurement result, -245 lV=K.
From Eq. ( 7), we obtain S ¼ -243 lV=K, which is consistent with the numerical measurement result.To derive Eq. ( 6) or Eq. ( 7), we assumed that T e ¼ T L , which holds in the strong scattering limit.In general, however, it may be hard to define T e x ð Þ in the device, because electron energy may not be fully relaxed by e-ph scattering, so that electrons and phonons are not fully in equilibrium. 8ur results above, however, show that the assumption works quite well when analyzing the simulation results, where we numerically "measure" S across the terminals and do not define or assume any T e x ð Þ within the device.As demonstrated in Fig. 4(c), for a constant j L , the heat diffusion equation gives a linear T L x ð Þ for a finite DT.In such cases, Eq. ( 7) can be simplified as where L is the total length of the device.Equation ( 8) implies that, for a uniform j L , the overall S of a composite structure is dominated by the region with a longer length.Figure 6(a) shows the simulation results for E h i versus x for our model device with optical phonon scattering (D 0 ¼ 0.01 eV 2 , hx o ¼ 20 meV) for various L w values with L b ¼ 10 nm (DV ¼ 10 mV, T 1 ¼ T 2 ¼ 300 K).The maximum E h i in the barrier region remains the same, but E h i in the well region decays more as L w increases.The energy relaxation length k E is estimated to be about 30 nm in our model device, and we expect that, as L w becomes much longer than k E , E h i deep in the well region will approach the value of a wire with a uniform doping density of N w .We numerically measure the overall S using the approaches in Fig. 2, and the results are shown in Fig. 6(b).Note that the results from the two approaches are still consistent.For L w <$ k E , the carrier energy is not fully relaxed in the well region, so that the overall S is more dominated by the barrier, resulting in a higher S j j.As L w increases, however, carrier energy is relaxed more in the well region, as shown in Fig. 6(a), and the overall S is more dominated by L w and approaches the measured S of a wire with a uniform, high doping density N w (no potential barrier).Note that the values of S for bulk wires in Fig. 6(b) are consistent with experimental results reported for heavily doped NWs. 13,14n summary, for a composite structure with no energy relaxation or L ( k E , the overall S is determined by the highest barrier within the device.For a diffusive composite structure with energy relaxation, however, the overall S is an average of its constituent S values, with the weighting factor being the lattice temperature drop across each region.For a simple case of a constant j L , T L x ð Þ is linear for a finite DT, and the overall S is the average of S x ð Þ over the device length.In Sec.IV, we discuss the meaning of the measured S of a composite structure in the context of TE performance, explore the effect of a x-dependent j L , and suggest possible ways to improve the TE performance by using the composite nano-structures.

A. Measured S and TE performance
As shown in the previous section, the two approaches to measure S in Fig. 2 give consistent results.Note that this measured S directly reflects the power generation performance of a TE device. 29,30In Fig. 2(a), the measured S is related to the open circuit voltage generated by a temperature difference as DV ¼ ÀSDT or, equivalently in Fig. 2(b), it represents the device's ability to drive an electrical current for a given temperature difference as I ¼ SGDT.For a homogeneous structure with a uniform S, the same S also determines the cooling performance.The cooling performance is determined by the heat current I q taken away from contact 1, i.e., I q x ¼ 0 ð Þfor DV > 0 in Fig. 1(a), and this is related to the P near the contact as I q x ¼ 0 ð Þ¼P x ¼ 0 ð ÞI, where I is uniform along the device.For a composite structure without energy relaxation, P x ð Þ is uniform, as shown in Fig. 3, and the uniform Seebeck coefficient from S ¼ P=T determines both the cooling and power generation performances.For a composite structure with energy relaxing scattering, however, P x ð Þ is non-uniform, as shown in Figs. 4 and 6, and cooling and power generation performances are represented by different S values.The measured overall S, which is related to the power generation performance, is an average of the non-uniform S x ð Þ ¼ P x ð Þ=T, and the weighting factor is the temperature gradient, as discussed in Sec.III.In the case of a constant j L (a uniform temperature gradient), the measured S can be related to the average Peltier coefficient (P avg ) as S ¼ P avg =T, where P avg ¼ Ð dxP x ð Þ=L.The cooling performance, however, is more related to the local value at the contact end, S x ¼ 0 ð Þ¼P x ¼ 0 ð Þ=T.

B. Effects of a non-uniform j L
In Sec.III, we discussed simulation results for a constant j L , which gives a linear T L x ð Þ for a finite DT.For a composite structure, j L may vary within the device, which alters T L x ð Þ and may affect the measured S. Figure 7(a) shows the simulation results for T L x ð Þ (DT ¼ 1 K, DV ¼ 0) for the two model devices, where one of them has a constant j L of 150 W/m-K (device 1) and the other one has a much lower j L of 1.5 W/m-K in the barrier region (device 2).Both devices are assumed to be diffusive with optical phonon scattering (D 0 ¼ 0.01 eV 2 , hx o ¼ 20 meV).For a constant j L (device 1), T L x ð Þ is linear, as already shown in Sec.III.When j L x ð Þ is non-uniform (device 2), however, a large portion of DT is applied across the barrier region, which has a smaller j L . Figure 7(b) shows simulation results for G versus L w for the two model devices (DV ¼ 1 mV, T 1 ¼ T 2 ¼ 300 K).Note that all the parameters are the same for devices 1 and 2 except for j L , so they give similar G values.In Fig. 7(c), however, we see a significant difference for the measured S. As shown in Eq. ( 7), the temperature gradient dT L x ð Þ=dx is the weighting factor in calculating the average of S x ð Þ.This means that the region with a larger DT dominates in determining the overall S. As shown in Fig. 7(a), for device 2 (a lower j L in the barrier region), most of DT is applied across the barrier, which has a high S x ð Þ j j, as shown in Fig. 6(a).Therefore, as shown in Fig. 7(c), the overall S j j is dominated by the barrier region and remains high even for L w > k E for device 2. For device 1 (a uniform j L ), however, the overall S is more dominated by the well region as L w increases, as discussed in Eq. ( 8), so S j j decreases and approaches the value of a uniform wire with a high doping density N w , as shown in Figs.6(c) and 7(c).These results suggest that, for a properly designed composite structure, it may be possible to alter the S versus G characteristic from its value in the bulk.

C. Improving the power factor
The effect of non-uniform j L on the measured S may open up new possibilities to improve the thermoelectric power factor S 2 r.It has been suggested that the power factor can be improved by using composite nano-structures composed of grains and grain boundaries, 17,31,32 where the grain is a doped crystalline region and grain boundaries are thought of as potential barriers that can result in the socalled "energy filtering" effect. 33The basic idea is the following: first, if the grain size d is much longer than the momentum relaxation length k p , then introducing grain boundaries within the device may not decrease r much, because the device is already in the diffusive limit.Note that the potential barriers at the grain boundaries filter out low energy carriers.If d is shorter or comparable to the energy relaxation length k E , then the carrier energy is not fully relaxed within the grain, so S j j increases, as also shown in our simulation results in Fig. 6.Note that k E is usually significantly longer than k p , 8 so, by engineering d as k p < d < k E , S may improve, while not hurting r much, which may result in a net improvement in S 2 r over a bulk material.If d is very large, as d ) k E , then one may expect that the composite structure will behave similarly to the bulk material.If the j L is non-uniform and specifically low in the barrier region, however, S j j may remain high even for d ) k E , due to the non-uniform T L -distribution, as shown in Fig. 7.This suggests that the use of polycrystalline materials may enhance S and the power factor, even for large grain sizes, if the grain boundaries highly impede phonon transport. 34,35A detailed quantitative study of the power factor enhancement in composite nano-structure is beyond the scope of this paper, but our results show that it is essential to treat the electron transport and phonon heat transport self-consistently to clearly understand the TE properties of composite structures and explore possible ways to improve the power factor performance.

V. CONCLUSIONS
In this work, we computationally explored the Seebeck coefficient of composite nano-structure within a self-consistent electro-thermal transport simulation framework.Quantum transport of electrons was treated using the non-equilibrium Green's function method coupled with the Poisson equation, and electron transport was solved self-consistently with the lattice heat diffusion equation.We numerically "measured" Seebeck coefficients using the techniques that mimic experimental methods and explored the effects of energy relaxing scattering and coupling with phonon heat transport.Simulation results show that, without energy relaxing scattering, the overall S of a composite structure is determined by the highest barrier within the device.For a diffusive composite structure with energy relaxation due to e-ph interaction, however, Peltier and Seebeck coefficients are position-dependent.The measured overall S is an average of the position-dependent values, and the weighting factor is the carrier temperature gradient, which follows the lattice temperature gradient due to e-ph interaction.Therefore, self-consistent simulations of the phonon heat transport and the resulting lattice temperature distributions are very important to understand the TE properties of a diffusive composite structure.
We also clarified the meaning of the measured S regarding the TE performance.For a diffusive composite structure, the measured S directly reflects the electrical power generation performance, but it is less related to the cooling performance.Our simulations also suggest that j L -distribution within the composite structure may be engineered to modify the S versus r trade-off and improve the power factor.For example, by making j L smaller in the region with a high local Seebeck coefficient (e.g., grain boundaries), we can apply a larger temperature gradient across that region and improve the overall S while not decreasing r.Note that this idea may work even for a large grain size of d > k E .
In this work, we used a simple 1D model and a classical heat diffusion model to treat the lattice heat transport, but we believe that the general understanding established in this paper should be useful in exploring nano-composite structure as a promising TE material.For future work, a more advanced phonon transport model 34,35 may be required to better treat non-equilibrium phonon transport in nanoscale devices and explore its effect on electron transport.While we assumed potential barriers with smooth interfaces where lateral momentum is conserved, 36 interface roughness scattering 37 may result in different thermoelectric properties.It will be also important to understand carrier transport in the network of 3D grains to address issues in realistic, bulk nano-composite materials.

APPENDIX: DEVICE MODELS 1.1D Poisson
The 1D Poisson equation along the x-direction is solved as where A is the cross-sectional area of the wire, V j is the electric potential at the jth grid point, a is the grid size, N is the total number of grid points, e S is the dielectric constant (e S ¼ 10 for our model device), e 0 is the vacuum permittivity, N 3D is the 3D doping density (N w or N b ), and n 1D is the 1D carrier density, which is obtained from the 1D NEGF, as discussed in Subsection Appendix B. For numerical stability, we adopt the non-linear Poisson technique 38 when solving Eq. (A1a).Equation (A1b) represents the boundary conditions at the device ends.

2.1D NEGF
The general model for dissipative quantum transport is 7 where G E ð Þ is the retarded Green's function, I is the identity matrix, H is the device Hamiltonian, and U ¼ ÀeV is the self-consistent potential energy obtained from the Poisson scheme in Eq. (A1).In Eq.
ð Þ is the spectral function, and C E ð Þ is the broadening, which are all N Â N matrices.For our model device, H is constructed using the effective mass approach 7 with m Ã ¼ 0.25m 0 and a ¼ 0.25 nm, where m 0 is the free electron mass.The matrices, R and R in , describe the effects of contacts and carrier scattering as , where the subscripts 1, 2, and s represent contact 1, contact 2, and scattering, respectively.
In this work, we consider three transport models, i.e., ballistic transport, elastic scattering, and inelastic scattering.For ballistic transport, R s and R in s are all zero in Eq. (A2).Elastic scattering can be treated as 7 where R s E ð Þ ¼ ÀiC s E ð Þ=2 and D 0 is the deformation potential.For acoustic phonon scattering with elastic approximation, D 0 in Eq. (A3) can be related to the acoustic phonon deformation potential D A as 28 where F is the wavefunction overlap, 8 q is the mass density, and t S is the sound velocity.In Eq. (A4), D A is given in eV and frequently appears in the conventional BTE approach. 8or Si bulk parameters, 39 Eq. (A4) gives D 0 $ 0.002 eV 2 for a cylindrical NW with D ¼ 3 nm and intraband transition within the ground state (F ' 2:66a 2 =D 2 ).In Sec.III, we use D 0 ¼ 0.01 eV 2 for elastic scattering, which gives k p $ 6.5 nm for our model device.
For inelastic scattering, we consider optical phonon with a single frequency x o as 7 where N x o is the Bose-Einstein factor and G p E ð Þ is the hole correlation function, which gives Note that, in Eq. (A5c), we assume an equilibrium occupation factor for phonons, but T L is position-dependent.We can also relate D 0 in Eq. (A5) to the optical phonon deformation potential D o as 28 where D o is given in eV=cm.For a cylindrical NW with D ¼ 3 nm and Si bulk parameters (LO mode), 39 we obtain D 0 $ 0.01 eV 2 .In Sec.III, we use this typical value of D 0 ¼ 0.01 eV 2 for optical phonon scattering, which gives k p $ 4.5 nm and k E $ 30 nm for our model device.
Once the solutions for G E ð Þ and G n E ð Þare obtained, physical quantities at the jth grid point can be calculated as 7,23 where I j!jþ1 is the electrical current flow from the jth to the (j þ 1)th grid points.Note that the spin degeneracy of 2 is included in Eq. (A7).The n 1D in Eq. (A7a) is again an input to the Poisson scheme in Eq. (A1), and the process is repeated until the self-consistent solutions for n 1D and V are obtained.

Lattice heat transport
To treat the heat conduction due to phonons, we solve a 1D heat diffusion equation as 21,25 d dx Àj L dT L dx À Á ¼ P; P ¼ ÀdI E =dx ; (A8a) where I E is the energy current by electrons, which is calculated from Eq. ( A7b) with e substituted by E. As shown in Eq. (A8b), we use fixed boundary conditions at the device ends.As described in Fig. 1(b), T L is calculated by solving Eq. (A8) for a given P from the electron part, and the updated T L again affects the electron part by changing the phonon scattering rates in Eq. (A5).The process is repeated until we obtain converged results for T L .

Figure 1 (
Figure 1(a) shows the schematic of our model device.We assume a 1D NW with the doping densities varying along the transport (x) direction.In Fig.1(a), L w is the length of the well region, L b is the barrier thickness, N w is the doping density in the well, and N b is the doping density in the barrier.The wire ends are connected to contact 1 and contact 2, which are ideal reservoirs maintained under equilibrium.7To convert the three-dimensional (3D) doping densities to 1D and solve the Poisson equation, we assume a circular cross-section with a diameter D. In all following simulation results, we use N w ¼ 3 Â 10 19 cm À3 (heavily doped), N b ¼ 1 Â 10 16 cm À3 (essentially undoped), and D ¼ 3 nm.The model structure is intended to represent a potential barrier (grain boundary) between two doped regions (grains).In a realistic nano-composite material, the potential barrier comes from a charge at the grain boundary due to point defects, etc.17 The simulation framework for the self-consistent electro-thermal transport is summarized in Fig.1(b).The "electron part" treating the carrier transport and electrostatics is solved self-consistently with the "phonon part" that describes the lattice heat conduction.Within the electron (b), the electrostatics are captured by solving the 1D Poisson equation, and the self-consistent solutions for the NEGF and the Poisson equation give the carrier density n and potential profile V along the x-direction.More details of the 1D NEGF and Poisson schemes are discussed in the Appendix.

FIG. 1 .
FIG. 1. (Color online) (a) Schematic of the 1D NW model device.L w is the length of the well region, L b is the barrier thickness, N w is the doping density in the well, N b is the doping density in the barrier (N w > N b ), D is the wire diameter, and x is the transport direction.The device is connected to ideal reservoirs: contact 1 and contact 2. (b) Simulation framework for the selfconsistent electro-thermal transport.The "electron part" calculates the self-consistent carrier density n and electric potential V, and the "phonon part" is solved for the self-consistent solutions for the power dissipation P and the lattice temperature T L .

Fig. 2 (
Fig. 2(a) (current measurement), we simply apply DV or DT and calculate the terminal electrical currents.For the approach in Fig. 2(b) (open-circuit voltage measurement), we first apply DT and then increase DV until there is no net electrical current flow.Under all DV and DT conditions, simulations are carried out self-consistently, as described in Fig. 1(b).

FIG. 2 .
FIG. 2. (Color online) Configurations for Seebeck coefficient measurement.(a) S from current measurements.I ¼ GDV for a finite DV(DT ¼ 0) and I ¼ SGDT for a finite DT (DV ¼ 0) and S ¼ SG=G.(b) S from the voltage measurement.For an open circuit voltage DV for a DT, S ¼ ÀDV=DT.

FIG. 3 .
FIG. 3. (Color online) (a) Ballistic transport simulation results for I E; x ð Þ for L b ¼ L w ¼ 10 nm, DV ¼ 1 mV, and T 1 ¼ T 2 ¼ 300 K.I E; x ð Þand the average E of the current flow E h i are uniform along the x-direction.Using the approaches in Fig. 2, S is determined to be -346 lV=K (from current measurements) and -348 lV=K (from the voltage measurement).The Kelvin relation gives S ¼ P=T ¼ E À E F h i =qT ¼ -344 lV=K, which is consistent with the numerical measurement results.(b) Simulation results for I E; x ð Þ for elastic scattering with D 0 ¼ 0.01 eV 2 .I E; x ð Þ and E h i are still uniform.The two approaches in Fig. 2 give consistent results of S ¼ -317 lV=K.The Kelvin relation gives S ¼ P=T ¼ E À E F h i =qT ¼ -315 lV=K, which is consistent with the numerical measurements.
FIG. 3. (Color online) (a) Ballistic transport simulation results for I E; x ð Þ for L b ¼ L w ¼ 10 nm, DV ¼ 1 mV, and T 1 ¼ T 2 ¼ 300 K.I E; x ð Þand the average E of the current flow E h i are uniform along the x-direction.Using the approaches in Fig. 2, S is determined to be -346 lV=K (from current measurements) and -348 lV=K (from the voltage measurement).The Kelvin relation gives S ¼ P=T ¼ E À E F h i =qT ¼ -344 lV=K, which is consistent with the numerical measurement results.(b) Simulation results for I E; x ð Þ for elastic scattering with D 0 ¼ 0.01 eV 2 .I E; x ð Þ and E h i are still uniform.The two approaches in Fig. 2 give consistent results of S ¼ -317 lV=K.The Kelvin relation gives S ¼ P=T ¼ E À E F h i =qT ¼ -315 lV=K, which is consistent with the numerical measurements.

FIG. 5 .
FIG. 5. (Color online) Seebeck coefficient measurement of a composite structure with region 1 (S 1 and G 1 ) and region 2 (S 2 and G 2 ).The device is assumed to be diffusive with energy relaxation due to e-ph scattering.(a) Open-circuit voltage measurement.(b) Current measurements.A simple electrical-thermal circuit analysis gives S ¼ S 1 DT L1 þ S 2 DT L2 ð Þ =DT for both cases, where DT L1 and DT L2 are the DT L values applied across region 1 and region 2 and DT L1 þ DT L2 ¼ DT.

FIG. 6 .
FIG. 6. (Color online) (a) Simulation results for E h i vs x with optical phonon scattering (D 0 ¼ 0.01 eV 2 , hx o ¼ 20 meV) for various L w values, with L b ¼ 10 nm and j L ¼ 150 W/m-K (DV ¼ 10 mV, T 1 ¼ T 2 ¼ 300 K).The maximum E h i in the barrier region remains the same, but E h i in the well region decays more as L w increases.(b) Simulation results for the overall S vs L w for the composite structure and bulk wire with no potential barriers.Two approaches in Fig. 2 still give consistent results (solid lines: open-circuit voltage measurements; dashed lines: current measurements).As L w increases, the overall S is more dominated by L w and approaches the value of a wire with a uniform doping density.

FIG. 7 .
FIG. 7. (Color online) (a) Simulation results for T L vs x (DT ¼ 1 K, DV ¼ 0) for the two model devices (L w ¼ 50 nm, L b ¼ 10 nm), where "device 1" has a constant j L of 150 W/m-K and "device 2" has a lower j L of 1.5 W/m-K in the barrier region.Both devices are diffusive with optical phonon scattering (D 0 ¼ 0.01 eV 2 , hx o ¼ 20 meV).For device 1, T L x ð Þ is linear, and for device 2, a large portion of DT is applied across the barrier region.(b) Simulation results for G vs L w for the two model devices (DV ¼ 1 mV, T 1 ¼ T 2 ¼ 300 K).They give similar G values, because all the parameters are the same except for j L .(c) Simulation results for the measured S vs L w (from opencircuit voltage measurements).The measured S is significantly higher for device 2 and stays high even for L w > k E , because it is dominated by the barrier region with a large DT, as shown in Fig. 7(a), and the barrier has a high S x ð Þ j j, as shown in Fig. 6(a).