Network theory for inhomogeneous thermoelectrics

The Onsager-de Groot-Callen transport theory, implemented as a network model, is used to simulate the transient Harman method, which is widely used experimentally to determine all thermoelectric transport coefficients in a single measurement setup. It is shown that this method systematically overestimates the Seebeck coefficient for samples composed of two different materials. As a consequence, the figure of merit is also overestimated, if the thermal coupling of the measurement setup to the environment is weak. For a mixture of metal and semiconductor particles near metal percolation the figure of merit obtained by the Harman method is more than 100 \% too large. For a correct interpretation of the experimental data, information on composition and microstructure of the sample are indispensable.


I. INTRODUCTION
Thermoelectric materials are important for energy harvesting, especially from waste heat 1,2 . In order to optimize the conversion into electricity, it is desirable to predict device properties by efficient computer simulations. This is one goal of the present paper. More fundamentally, we are going to point out that memory effects render the global response of composite materials, which are common among modern nanostructured thermoelectrics 3 , highly complex.
The theoretical description of transport processes goes back to Onsager 4,5 . Applied to thermoelectrics, this became known as the Onsager-de Groot-Callen theory 6 . A special case is the so-called constant property model (CPM), where the electric and heat conductivity, σ and κ, and the Seebeck coefficient α are assumed to be constant. This model has been studied in detail analytically in one dimension 7,8 and will serve as a reference system for validation in this paper.
Since analytic calculations are restricted to simple compounds, numerical models have been developed in order to describe inhomogeneous materials. Although these models are also based on the Onsager-de Groot-Callen theory, most of them do not fully describe the thermoelectric effects, since they do not include Joule heat and/or Peltier heat 9,10 . They were used to calculate either the heat and electrical conductances or the Seebeck coefficient [11][12][13] . More complex models 14 exist in the framework of drift-diffusion models, and they are applied e.g. for the simulation of generators in complex geometries 15 .
In this paper a simple and coherent way to discretize the Onsager transport theory for thermoelectric materials will be used, which includes all relevant effects and time-dependencies. It can be seen as a version of the finite difference method reviewed in 16 . Originally it was designed for the investigation of transport processes in nano-particle configurations ( 17 , 13 , 18 ), which were mapped to a network model. In this paper it will be applied to disordered bulk systems rather than particle agglomerates.
The model, derived in Sec. II, enables us to study thermoelectric effects in geometries, which can not be solved analytically. It is validated by comparing simulations and analytical results for one-dimensional segmented CPM thermoelectrics. As an application, the transient Harman method is simulated focusing on inhomogeneous systems like segmented thermoelectrics, superlattices and composite media. Furthermore, we support and extend our findings concerning segmented structures by an analytical treatment.

II. MODEL
The Onsager-de Groot-Callen theory 8 describes the flow of electrical current density j and heat current density j q j = −σ (∇µ/q + α∇T ) (1) j q = −κ∇T + Π j (2) driven by the gradients of temperature T and electrochemical potential µ (q is the charge of the mobile particles). According to the Kelvin relation the Peltier coefficient Π is related to the Seebeck coefficient α by Π = αT . If the Seebeck coefficient is zero, heat and electric transport decouple resulting in Ohm's and Fourier's laws with electrical conductivity σ and heat conductivity κ. Volumetric heat production due to Joule heating taken into account by with heat capacity per volume, c. Discretization of these equations in form of a network model is achieved by assigning variable temperatures T i and electrochemical potentials µ i to each lattice site i. The bonds between neighboring sites are characterized by an electric conductance G ij , a heat conductance K ij , a Seebeck and a Peltier coefficient, α ij and Π ij . In the network model the electric current I ij and the heat current I q,ij between sites i and j read The coefficients G ij , K ij , α ij , Π ij depend on the material properties at sites i and j. Let G i denote the conductance of the material site i. The electrical resistance 1/G ij of the bond between sites i and j, half of which is material i or j, respectively, is modeled as with an extra interface contribution R 12 . Likewise, the thermal conductance K ij is calculated from Being contact rather than material properties, the interface contributions R and R q are particularly important, if one has only one crystalline material that is divided into grains by grain boundaries. In this paper, however, we want to focus on compound materials, where the different bulk properties have a strong influence on the thermoelectric response. In order not to obscure this bulk effect by the additional influence of the interfaces, we set R = 0 and R q = 0 in the following. Next, we discuss the Seebeck coefficients α ij , if the materials at the neighboring sites i and j are different. The Seebeck voltage α ij (T i − T j ) is the sum of two contributions: first, the voltage between site i and the interface in the middle of the bond ij, which is at temperature T c,ij , and second, the voltage between the interface and site j: whereᾱ ij = (α i + α j )/2,T ij = (T i + T j )/2 and ∆α ij = (α i − α j ). ApproximatingT ij = T c,ij , eq. (8) simplifies to For the Peltier coefficient Π ij , finally, one has to take the Peltier heat properly into account, which is delivered to or taken from the adjacent materials i and j at the interface. The electrical current I ij carries the heat Π i I ij away from site i and delivers Π j I ij to site j. We assume that the energy difference (Π i − Π j )I ij , which is set free (or consumed) at the interface, is given to (or taken from) both adjacent sites in equal parts. The net heat current induced by I ij is then with Due to the Kelvin relation, Π i = T i α i for material i.
For the temporal evolution of the local temperatures, which is essential for the Harman method, the heat capacities C i must be given for the material belonging to site i. Without specifying the material further, the ratio C/K is the typical time scale, on which temperature differences between neighboring sites are leveled out in the network by heat conduction. On the other hand, the electrical site capacitances C el can be neglected, because the equilibration time C el /G of the electrochemical potential is much shorter than C/K. Therefore it is a reasonable approximation that the electrostatic potential adjusts instantaneously compared to the slow thermal evolution of the system.
The time evolution of the temperature T i is described byṪ The second term is the Joule heat (eq. (3)) set free on bond ij. The factor 1/2 accounts for its delivery in equal parts to sites i and j. For the numerical integration the time step is about one-hundredth of the time scale of the fastest heat exchanging mechanism C/K. The simulation procedure works as follows: An initial temperature distribution and a total current I entering one electrode and leaving the other are provided (see Fig. 1). The electrochemical potentials are calculated from Kirchhoff's first law and then fed into eq. (13) in order to calculate the temperatures for the next time step. The temporal evolution of the temperatures implies continuous adjustments of the local currents (according to Eq. (4)) and hence the electrochemical potentials. Iterating these steps gives the transients and finally the stationary state.
Electrodes are attached to the opposite surfaces normal to the x-axis. They are assumed to have uniform electrochemical potential and temperature. In this paper the electrical current I flowing through the material (and the electrodes) is fixed. Heat losses through the surfaces  parallel to the overall current direction are not taken into account. Two different boundary conditions for the temperature at the electrodes are considered. In Sections III A and III B the temperatures T 0 > T L = 300 K of the electrodes at x = 0 and x = L are fixed. We denote this as "setup C". In the Sections starting with III C, a different "setup H" is used (see Fig. 1): The electrodes are in contact with the environment, which acts as a heat bath of fixed temperature T HB = 300 K. The temperatures of the electrodes are determined by Joule heating, Peltier effect and heat exchange with the environment. In general they differ from T HB and depend on the current direction.

A. Constant property model (CPM)
Here we briefly recall the analytical solution of the onedimensional CPM 7,8 . Equations (1) -(3) yield in steady state where j represents the x-component of the electrical current density. For fixed temperatures T 0 and T L the solution reads The symmetry of the j 2 -term with respect to x = L/2 implies that the Joule heat produced in the sample is flowing equally to both electrodes. Similarly, the electrochemical potential profile µ(x) follows from eq. (1) and eq. (15) Without loss of generality, the integration constant in (16) is µ(0) = µ 0 = 0.

B. Segmented thermoelectrics
Segmented thermoelectrics consist of two layers of different material. If represented by a one-dimensional piecewise constant property model, analytic expressions for the temperature and potential profile, as well as for the effective transport coefficients can be obtained 20 . These expressions shall be compared to the corresponding simulation results in order to validate the model presented in Sec. II.
The analytical calculation makes use of the fact that for each segment Eq. (14) holds, so that the temperature profile is piecewise parabolic with the respective parameter sets. The six integration constants (for each material two from Eq. (14) and one from ∂µ/∂x) are fixed by the boundary conditions T 0 , T L and µ 0 = 0 and by requiring continuity of T (x), µ(x) and j q,x (x) at the interface. Figure 2 shows the temperature profile for setup C with the parameters given in table I. The electrodes are characterized by the same parameters as the adjacent material, which prevents additional Peltier heating/cooling at the electrode-sample interfaces. The size of the sample is L = L x = L y = L z = 10 −2 m, discretized by N y = N z = 1, N x = 100 lattice sites. The current |I| = 50 A flows from left to right (squares) or in reverse direction (blue points). The effect of Peltier heating/cooling at the interface can be clearly seen: The temperature profiles differ significantly for the two current directions. It is this internal temperature profile which affects the results of the Harman method, as will be discussed in Sec. III D.
Considering the bond at the interface where ∆α ij = 0, we note that eq. (9) is an approximationl, if a non-linear temperature distributions is present. It modifies the total voltage by an interface term, the relative contribution of which vanishes if N x becomes large. As we discarded interface resistances previously, this interface effect may be neglected as well.

C. Harman method
The Harman method 21,22 is a measurement technique, which allows to determine the three transport parameters and consequently the figure of merit within one single measurement procedure. A setup H as shown in fig. 1 is used, and a known dc current I is applied. Due to the Peltier effect one electrode/sample boundary heats up and the other one cools down. At the same time, Joule heating leads to an additional change of the sample temperature. This continues until a steady state is reached, where heat conduction compensates further temperature changes. The time evolution of the electrode temperatures and the voltage across the sample (see fig. 3) are measured. After reaching the steady state the current is switched off. Starting with a homogeneous temperature the electric conductance can be determined from a measurement of the voltage V σ (see fig. 3) via Primed quantities denote Harman method measurement results. In experiments this measurement might be difficult, as the Peltier effect starts immediately creating a temperature difference between the electrodes, when the current is switched on 22 . However, determining G tot during simulations is simpler, since V σ is recorded before the temperature change affects the voltage. The Seebeck coefficient is calculated from the temperature difference ∆T and the Seebeck voltage V α measured immediately after switching off the current I.
The heat conductance of the sample, is obtained from energy balance in the steady state 22 .
Here, K env denotes the heat conductance between the heat bath and the electrodes. The bars indicate averaging of the values at x = 0 and x = L, e.g.T = (T 0 + T L )/2. µ env,0 and µ env,L are the electrochemical potentials of the leads next to the left and right electrode, respectively. Let us briefly recall the derivation of (20) according to 22 . The energy current I e = I q + µ/qI arriving at the left electrode from the environment (heat bath or lead) is In steady state it must be equal to the energy current leaving the sample on the right hand side. In these expressions, the last terms are the fractions of Joule heat produced in the leads, that flow into the sample. They are assumed to be equal for simplicity. The Peltier coefficient of the (metallic) leads is assumed to be negligibly small. Taking the average of these expressions eliminates the parameters T HB and G env : On the other hand, the energy current entering the sample from the left must be equal to (21). In the framework of the constant property model it is given by where the subscripts refer to the values at x = 0 or x = L, respectively, as before. This in turn must be equal to the energy current arriving at the right electrode from the sample, Again taking the average of these two expressions eliminates the parameter G tot : Identifying this with (23) leads to (20). (μ −μ env )/q is the average contact potential between the electrodes and the leads. Here we do not take contact resistances into account. Hence it does not enter the evaluation of the simulation data. In experiments additional losses due • -12  to convection and heat radiation occur, which can be taken into account via correction terms added to eq. (20) (e.g. 23 ).
The Harman method gives correct results for homogeneous samples. However, as we are going to explain in the next two sections, when applied to heterogeneous systems, it can give values for the Seebeck coefficient, which can be more than 100% off the true value. Consequently, also the heat conductance derived from (20) will be misleading.

D. The Harman method applied to segmented thermoelectrics
The basic problem encountered, when the Harman method is applied to inhomogeneous systems, will be explained in this section for the example of a bilayer thermoelectric. We simulated a setup H ( fig. 1) consisting of materials A and B with parameters given in table I. The volume fraction of the A-layer is denoted by f . The elec-trodes are characterized by the same parameters as the adjacent layers, except of the Seebeck coefficient, which is set to zero. The dimensions are chosen as in sec. III B.
The transport coefficients are measured with the Harman method for different ratios f and currents I. The results are then compared to the analytical expressions, which will be given first.
Denoting the electrochemical potential at the interface betweem the A-and B-segment by µ AB , the electrical conductivity of material A in the absence of a temperature gradient is The heat conductivity, on the other hand, is defined under open circuit conditions, j = 0, as with the temperature T AB at the interface. The expressions for σ B and κ B read accordingly. The effective elec- trical and heat conductivities are then derived as a connection in series: The effective Seebeck coefficient, also defined for open circuit conditions, corresponds to a series connection of the Seebeck voltages created from material A and B: Fig. 4 shows a comparison between eqs. (29), (30), (31) and the results of a simulation of the Harman method. As expected, the electric conductivity obtained from the Harman method coincides perfectly with the analytic results ( fig. 4, upper left). But the Seebeck coefficient (upper right) and the heat conductivity measured via the Harman method depend on the current, which was applied before the measurement of V α and deviate strongly from the open circuit values, eqs. (30) and (31). Only in the limiting cases f → 1 and f → 0 the differences vanish. Consequently, the same is true for z eff .
Considering the Seebeck coefficient, the deviations arise from the influence of the Peltier heating/cooling at the interfaces on the temperature profile (see fig. 5). It causes the current dependence of the temperatures T 0 , T L , T AB and hence of the thermopower If one compares the analytic temperature profiles for f = 1/2 created by the external current I = −12 A (black, solid line) and I = 12 A (gray, solid line) to the Although the Seebeck coefficient is measured at I = 0 when using the Harman method, it depends strongly on the previous current (see fig. 6), because the temperature profile in inhomogeneous samples acts as a memory. Remarkably, even for arbitrarily small currents |I| the Harman measurement gives a Seebeck coefficient which deviates from the true open circuit value. In this limit Joule heating is negligible and the temperature at all three interfaces depend linearly on the current due to Peltier interface effect. Hence, the relation between the temperature differences across both materials becomes independent of I for I → 0.
We elaborate this argument by deriving the Seebeck coefficient of a segmented thermoelectric measured by the Harman method analytically. For this purpose we consider the AB structure depicted in fig. 7, which is connected to heat bathes with fixed temperature T HB via heat conductance K env . For each segment i the temperature follows eq. (32), which coincides with the simulations (lines in upper right of fig. 4). Usually, small currents are applied allowing for a expansion of T (x) around j = 0, which leads to with the cross section area of the sample S. It is important to note that the difference always has the same sign as α eff meaning that the Harman method systematically overestimates the absolute value of the Seebeck coefficient.
Measuring the heat conductivity according to eq. (20) leads to flawed results in a segmented thermoelectric. However, we present a possibility to infer the heat conductivity eq. (30) from the Harman method requiring α eff , which can be gained from Harman measurments using eq. (34). Basically, we repeat the derivation for K tot as presented in sec. III C. We determine the energy currents at x = 0 and x = L using the temperature distribution and its derivative as derived for the strucure depicted in fig. 7. The resulting energy currents are expanded to first order around j = 0 and their average reads Following the arguments of sec. III C this must be equal to eq. (23), which leads to where we have set µ 0 − µ env,0 = µ L − µ env,L = 0. Using with α c = α eff −α eff (see eq. (34) and κ c = α c T HB j/(T L − T 0 ). According to eq. (37) the figure of merit may be over or underestimated by the Harman method. Just for K env → 0 we yield Hence, for small currents and weak heat coupling the figure of merit is always overestimated by the Harman technique. Considering the relative deviation ∆ z = (z eff − z eff )/z eff in the limit K env → 0 we find that its maximum occuring at the most unfavorable f = f max is which is solely affected by the difference of both Seebeck coefficients and can be very large. E.g. for α A = 0.0004 V/K and α B = 0.0001 V/K we get ∆ z (f max ) = 0.5625 and the Harman method overestimates the z by 56.25%. We conclude that applying the Harman method to segmented structures as done in 24 is tricky and may lead to results that are too optimistic. Now we show that the systematic error, which we found, if the Harman method is applied to a doublelayer system, remains unchanged for periodic superlattices in the limit of weak current. This is particularly easy to see for even numbers of layers. First, for open circuit conditions, a given temperature difference T L − T 0 will be evenly distributed among all N double layers, T ν+1 − T ν = (T L − T 0 )/N (see fig. 8). In first order this is also the case, if a current is imposed. In particular, the temperature dependence of the Peltier effect may then be neglected. Hence, each double layer has  the same systematic error of the effective Seebeck coefficient, eq. (34). Consequently the homogeneous array of double layers gives rise to the same error. If the number of layers is odd, this conclusion remains true up to a correction proportional to 1/N which vanishes for large superlattices.
Superlattices attracted a lot of interest and a record zT ≈ 2.4 was measured using the Harman method in Bi 2 Te 3 /Sb 2 Te 3 superlattices 25 , which has not been reproduced so far. The individual segments are very thin: Bi 2 Te 3 segments have a thickness of 1 nm and the Sb 2 Te 3 segments have a thickness of 5 nm implying f = 1/6. Neglecting quantum effects and using α BiTe ≈ 2.2 · 10 −4 V/K 26 , α SbTe ≈ 0.9 · 10 −4 V/K 27 , κ BiTe ≈ 2 W/(mK) 28  In contrast to a segmented thermoelectric, the number of interfaces between different materials is large and the error produced by the approximation eq. (9) might become relevant. We checked that this is not the case by comparing with a simulation, in which the Seebeck coefficient α ij for a bond between the two materials was determined from α i and α j weighted with K i and K j as in eq. (31). No significant differences could be detected. Therefore we keep eq. (9) for the following simulations.
For such disordered two-component systems, an effective medium theory was developed 29-31 with the following expression for the Seebeck coefficient: The effective medium electrical and heat conductivities are given by and an analogous equation for κ eff . The parameter A = (1 − f c )/f c is connected to the percolation threshold f c . The Harman method gives an accurate effective electrical conductivity for the composite material, as explained before. By fitting these simulation results with eq. 41, we determine the effective medium parameters f c = 0.594 (2) and t = 1.315 (8). Note that f c is close to the expected percolation threshold 0.592746 for site percolation on a square lattice. The effective medium theory expressions for σ eff and κ eff with the above fitting parameters are fed into eq. (40) to obtain the effective medium value for the Seebeck coefficient.
The Seebeck coefficient α eff obtained by the Harman method deviates strongly from the results for open circuit conditions and from the effective medium theory, especially slightly below the percolation threshold (see fig. 9). This phenomenon can be understood by looking at the temperature and current distributions at the moment the Seebeck coefficient is measured. For f = 0.8, far above the percolation threshold for the metal-like material, the temperature and current distribution appear homogeneous (see fig. 10a) leading to a good match of the Harman and the open circuit measurements. In another sample at f = 0.5 (see fig. 10b and 10c), however, a path of well conducting material A almost percolates and thus carries a majority of the current. The percolation is interrupted by material B and a strong Peltier heating/cooling appears. This together with differences in heat conductivity of material A and B results in a strongly inhomogeneous potential and temperature profile. A big part of the temperature drop is located across material B, which is characterized by α B > α A , hence the effective Seebeck coefficient is enhanced close to the percolation threshold (see fig. 9).
The overestimated Seebeck coefficient α eff in turn affects the heat conductivity shown in fig. 9 (right) leading to a an overestimation of κ eff by a factor larger than 3. Taking α eff determined by eq. (40) to derive κ eff with eq. (20) (yellow squares) a significantly better agreement with effective medium theory (black line) is achieved.

IV. SUMMARY
A simple but powerful simulation model has been derived that describes all thermoelectric responses according to the Onsager-de Groot-Callen transport theory. Although not discussed in the present work we would like to emphasize, that the generalizations to arbitrary (including three dimensional) lattices and the inclusion of charge dynamics 17 are straight forward.
By means of this model we pointed out that the Harman method to measure the transport coefficients experimentally has substantial systematic errors, when applied to composite materials. In its usual operation mode (weak imposed current, weak thermal coupling to the environment) it always overestimates the absolute value of the Seebeck coefficient and the figure of merit. We gave a numerical example, where the Seebeck coefficient turned out to be wrong by more than a factor of 2 ( fig. 9, f ≈ 0.5).
In order to explain this effect, we calculated it analytically for a superlattice of alternating layers of two different materials. The reason is that the Peltier heating/cooling at the interfaces enlarges the temperature gradient in the layers with the larger absolute value of the Seebeck coefficient, and reduces it in the other layers. This self organized correlation between Seebeck coefficient and temperature gradient persists, when the current is switched off. The temperature profile acts as a memory of the previous current and affects the thermal voltage, from which the Seebeck coefficient is inferred in the Harman method. The correct Seebeck coefficient would be measured under open circuit conditions for a different internal temperature profile that does not reflect any previous Peltier heating/cooling. The difference between the two measurements of the Seebeck coefficient could be calculated analytically for a superlattice and is given in eq. (34). It can be used to correct the Harman measurement.
Overestimating the Seebeck coefficient implies an overestimation of the heat conductivity in the Harman method, as well. However, for small currents we determined an expression, eq. (36), which represents the correct open circuit heat conductivity using quantities available during a Harman measurement.
Random composites are only accessible by simulations. In this paper we discussed a material composed of metal and semiconductor particles or domains ( fig. 10). Below the percolation threshold of the metal-like phase we confirmed that the Harman method overestimates α and κ.
In summary, the application of the Harman method for inhomogeneous media is tricky. However, for segmented thermoelectrics including superlattices results from the Harman method can be corrected.

V. ACKNOWLEDGMENT
We thank R. Chavez, G. Schierning and R. Schmechel for valuable discussions. Financial support by the German Research Foundation (DFG) within the Priority Program on nanoscaled thermoelectric materials, SPP 1386 is gratefully acknowledged.