Improved Electrical, Thermal, and Thermoelectric Properties Through Sample‐to‐Sample Fluctuations in Near‐Percolation Threshold Composite Materials

Effective medium theories (EMT) are powerful tools to calculate sample averaged thermoelectric material properties of composite materials. However, averaging over the heterogeneous spatial distribution of the phases can lead to incorrect estimates of the thermoelectric transport properties and the figure of merit ZT in compositions close to the percolation threshold. This is particularly true when the phases’ electronic properties are rather distinct leading to pronounced percolation effects. The authors propose an alternative model to calculate the thermoelectric properties of multi‐phased materials that are based on an expanded nodal analysis of random resistor networks (RRN). This method conserves the information about the morphology of the individual phases, allowing the study of the current paths through the phases and the influence of heterogeneous charge transport and cluster formation on the effective material properties of the composite. The authors show that in composites with strongly differing phases close to the percolation threshold the thermoelectric properties and the ZT value are always dominated exclusively by one phase or the other and never by an average of both. For these compositions, the individual samples display properties vastly different from EMT predictions and can be exploited for an increased thermoelectric performance.


Introduction
Printable highly efficient thermoelectric materials [1,2] offer a promising approach towards low-cost thermoelectric generators (TEG). [3] Such devices offer an elegant way to energy-harvest for autonomous sensor nodes, wearable electronics, and large-scale waste heat recovery. [4][5][6][7] TEGs operate by transforming thermal energy, usually from ambient or waste heat, directly into electrical energy. This energy conversion is based on the Seebeck effect, which produces a voltage V = S × ΔΤ when a temperature difference ΔΤ is applied to the opposite sides of a thermoelectric material. The Seebeck coefficient S is the proportionality factor of this conversion. Considering the electrical conductivity and the thermal conductivity of the material, one obtains its thermoelectric figure of merit: with T being the absolute average temperature of the sample. The ZT value is a direct measure of how efficiently a thermoelectric material can convert thermal into electrical energy. [8] In the field of printed thermoelectrics, composite materials with separated phases must be used as the printability of the materials needs a resin in which the semiconducting moiety is mixed. The focus in the field lies on composite systems of an efficient thermoelectric constituent and a resin in which a compromise has to be found that combines good thermoelectric efficiency with manufacturing properties such as printability, a low sintering temperature, and mechanical flexibility. [9][10][11][12][13][14][15] The real limitation for the ZT value of a composite is, in most practical cases, not the one given by the highest ZT value of its components, [16] but by their appropriately averaged ZT values. The latter results from the large contrast in the transport properties of the components. These mean values can be quantitatively determined using analytical approaches such as the generalized effective medium theory (GEMT, see below). [17] They easily give rise to figures of merit orders of magnitudes below the above bound. To identify an approach that can bypass the predictions of GEMT would give rise to a much wider spectrum of allowed transport properties and circumvent principles that usually reduce the figure of merit. As we will show, large sample-to-sample fluctuations near the percolation threshold allow creating samples with a sizable yield that have much more formidable thermoelectric performance than EMTs would imply.
In composite materials, the EMT [18,19] has been used to calculate the thermoelectric properties. It approximates the complex network of different phases with an effective homogeneous material and calculates its resulting macroscopic properties based on the individual properties of the components and their volume fractions. Due to the separated phases in composites with large contrasts in their local transport properties, the charge transport occurs heterogeneously. The current flow passes through the material, following the paths with the highest conductivity, while regions with low conductivity contribute less to the current. In the extreme limit of infinite contrast between the local conductivities, this behavior can be described as a percolation problem. Here a phase with a volume fraction higher than a specific threshold v c is connected through the entire network and dominates the electric or thermoelectric charge transport. The value of v c depends on the dimensionality as well as the connectivity (lattice type) of the specific network. [20,21] One disadvantage of EMT is that it does not account for the rich physics that occurs in the vicinity of the percolation transition. To address some of these shortcomings, a GEMT was introduced to properly account for the value of the percolation threshold and the microscopic structure and morphology of the composite. [17] While GEMT offers an excellent description of sample averaged values of S, , and hence is a significant improvement over the earlier EMT, it ignores the effects of sample-to-sample fluctuations due to percolation effects. [22,23] In this work, we investigate the role of sample-to-sample fluctuations in transport properties in composites with large material property contrasts. As expected for the behavior near a critical point, fluctuations become increasingly important in the vicinity of the percolation threshold. We demonstrate that most individual samples possess transport properties that are dramatically different from the mean values. This increased variety in the allowed transport properties can then be exploited to obtain samples with much improved thermoelectric performance. To this end, we use a numerical model that is based on RRNs to simulate the thermoelectric properties of composites. RRNs can be used as an alternative to state-of-the-art EMTs [22,[24][25][26][27][28][29] and have proven to be powerful for describing effects from charge transport and heat transport in composite media to fluid flow in porous materials. By discretizing space into finite elements (sites), we model the material as a large linear network. We then randomly assign thermoelectric transport coefficients to the bonds between the sites of a phase with a probability equal to its volume fraction. This way, we conserve the information about the spatial distribution of the local potential and the local temperature. By using an expanded nodal analysis method, we calculate the heterogeneous charge and heat transport which we evaluate to compute the macroscopic thermoelectric properties. We validate our approach by averaging the results of numerous random networks and comparing them to the GEMT. Further studies of the networks reveal that for some compositions the average deviation of the materials' properties gets extremely large so that averaging and therefore the GEMT is no longer a good approximation of the macroscopic materials' properties. We investigate the influence of percolation on the macroscopic ZT value and point out that the percolating phase with the higher conductivity dominates the charge transport and dictates the effective ZT value.

Effective Medium Theories
In composite materials, the thermoelectric properties can be calculated based on the properties of the components as well as their volume fraction v i in the composite. The most common method is the well-established EMT stating: [18] Here, is a placeholder for an effective material property of the composite and i denotes the property of phase i. Landauer proposed Equation (2) in 1952 with = for the calculation of the effective electrical conductivity of a multi-phase system. [19] One year earlier, Odelevskii applied a similar theory to the effective thermal conductivity ( = ). [30] Finally, it was shown that the EMT can also be used to calculate the effective Seebeck coefficient by choosing = /S. [31,32] To address some of the shortcomings of EMT, McLachlan et al. proposed in 1990 a GEMT. [17] By applying principles of percolation theory to the already existing EMT a modified selfconsistency equation as proposed:  10×10 thermoelectric RRNs. Each site i is described by its electrical potential V i and its temperature T i . The transport coefficients G ij , ij , and Κ ij belong to the bond between sites i and j and are assigned at random between two distinct transport mechanisms TM I with G I , I , Κ I , or TM II with G II , II , and Κ II . The input and output nodes are marked in grey.
grain structure and morphology. For v c = 1/3 and t = 1 we obtain the original Equation (2) and hence Equation (3) extends Equation (2). In their work, McLachlan et al. were able to show that all principles of the EMT still apply and they showed its validity for the electrical and thermal conductivity. Recently, Sonntag showed that the GEMT was also suitable for calculating the Seebeck coefficient analogously to the EMT by = /S. [33] The GEMT offers, therefore, a full set of equations to calculate the ZT value for any composition of a multiphase composite. The percolation threshold v c and the t parameter in the GEMT account for the influence of the percolation of charge transport in a composite on the effective material properties, but the parameters need to be known beforehand or can only be obtained by fitting experimental data. [34]

Percolation Theory of Infinite Networks
Percolation theory is a vast research field utilizing graph theory and probability theory. It studies the formation of clusters in a connected network. The premise is the random assignments of connections (bonds) on a grid until two opposing sides are connected by creating a percolation path between them. In the natural and engineering sciences, percolation is a useful concept to describe a variety of phenomena from fluid flow in porous systems to charge transport in disordered solids. [20,21,23] In the practical implementation, the bonds are set randomly in a lattice with a binary probability distribution of P bond = v and P empty = 1 − v that two sites are connected. The percolation threshold P = v c is defined as being the probability where for the first time a single cluster of one phase spans the entire sample from one side to the other. A percolation cluster is a fractal object giving rise to the crucial role of rare connecting paths that dominate the transport behavior. Changing only a few sites in the network can strongly affect those connecting transport paths, giving rise to large sample-to-sample fluctuations for large but finite samples. Bond percolation occurs when the connections are along the edges of the grid as opposed to jumping from tile to tile (site percolation). The value of v c depends on the lattice type (squared, cubic, tetragonal, spherical, etc.) as well as the dimensionality of the network. The correct choice of the lattice type is therefore relevant to model experimental data correctly. [20] However, the existence of strong fluctuations at the percolation threshold is universally present, regardless of the details of the lattice.
Within percolation theory, one only captures the limit of infinite contrast in the local conductances. In real systems, the fact that the conductivities are finite for all bonds requires a numerical investigation. Still strong percolation-like behavior is present when the system has a high distinctiveness between the compounds resulting in II >> I or II << I . [17,21]

Thermoelectric Random Resistor Networks
To model the disordered structure of thermoelectric composite materials, we choose a squared lattice of size n × m as a discretization of a 2D geometry (Figure 1). Each site i in the lattice has attributed a local electrical potential V i and a local temperature T i .
Two neighboring sites i and j connect through a bond with a specific transport mechanism (TM) defined by a set of thermoelectric transport coefficients G ij = ij , ij = ij S ij T, and Κ ij = ij Τ which are derived from the local thermoelectric properties of the phase in which the bond ij lies. [8] We can express the injected current density J i and the injected heat flux q i into site i by applying Kirchhoff's first law [35] for all neighboring sites j: www.advancedsciencenews.com www.advtheorysimul.com T m is the mean temperature of the network. To allow a parallel injection and extraction of current densities into the network we add two additional sites, the input, connected to all sites at the x = 0 line, and the output, connected all sites at x = n line. The result is a large linear network that can be solved by an algorithm similar to those used for solving linear electrical networks. [36] We have expanded the method to accommodate the Seebeck effect, the Peltier effect, and the thermal conductivity. We can, therefore, describe the thermoelectric network as the system of equations below: where ⃗ J represents all the source current densities and ⃗ q all the source heat flows into the sites and ⃗ V and ⃖⃖⃖⃖⃖ ⃗ ΔT are the electrical potentials and temperature differences to the electrical and thermal ground.Ĝ,Â, andK are sparse matrices containing the information of all local transport parameters given byĜ ij = −G ij andĜ ii = ∑ i≠j G ij as well as the entries forÂ andK analogously.
Solving the system of equations for different boundary conditions enables us to calculate the macroscopic material properties of the system. First, we inject an electrical current density J in = J 0 into the network while keeping the global temperature difference ΔT = ΔT in − ΔT out = 0. From the resulting voltage distribution, we obtain the effective electrical conductivity by: Correspondingly, the effective thermal conductivity equals: with the initial heat current q in = q 0 and V = V in − V out = 0. Furthermore, we apply an initial heat current q in = q 0 and keep ΔT = ΔT in − ΔT out = 0. This results in Equation (9) as an expression for the effective Seebeck coefficient: [8] with Equations (7)-(9) we obtain an additional set of equations to calculate the ZT value for a thermoelectric network. For the simulation of composite material, we can simply assign different transport properties to the bonds, defined by their respective transport coefficients. The transport coefficients are distributed randomly in space with a probability equal to the volume fractions v I and v II = 1 − v I of their corresponding phases.

Results and Discussion
For reasons of simplicity, we focus in this work on a finite 2D squared lattice, since it is enough to demonstrate the principles of percolation. The threshold for bond percolation in an infinite 2D square grid is v c,2D = 0.5. [37] With no loss of generality, we study a two-phase composite using the RRN method with the local thermoelectric properties I , S I , and I in phase I and II , S II , and II in phase II, respectively ( Table 1). We chose those parameters as well as T m = 300 K inspired by a printable thermoelectric ink and its use in an energy-harvesting device at room temperature.
For statistical significance, we simulated 1000 random configurations for each composition with a network size of 500 × 500. We calculated the thermoelectric properties according to Equations (7)-(9) and varied the respective fractions of the phases. We took the arithmetic means < >, <S>, and < > of all randomly cast networks as the effective thermoelectric properties of the composites. For comparison, we calculated the effective thermoelectric properties using the GEMT. We fitted Equation (3) to the RRN results and obtained t = 1.2 and v c = 0.5 in agreement with the established values for a 2D square lattice. [37] Figure 2a shows the results of both methods in comparison as a function of v II . Both methods deliver essentially identical results. This result is in agreement with literature. [22] The methods only differ in a negligible numerical error stemming from the finite meshing of the network.
It can, therefore, be argued that the averaged transport properties < >, <S>, and < > derived from RRN calculation are in full agreement with the results of GEMT (Equation (3)) for all transport parameters despite the strongly different ratios of the transport parameters of the different phases II / I . Thus, GEMT offers an excellent description of the averaged values of the relevant transport quantities, including the region in the vicinity of the percolation transition.
However, EMTs fail to describe the large sample-to-sample fluctuations that emerge near the percolation transition. To this end we compare the effective properties of each of the N = 1000 individual network to < > we calculate the mean relative deviation: of the corresponding transport coefficient . Figure 2b shows Δ , Δ S , and Δ versus v II . For compositions far away from the percolation threshold Δ is negligible for all material properties. The effective material properties of the individual networks are therefore well approximated by < > and the GEMT. A completely different picture results for v II close to v c . Very large deviations from < > for the electrical conductivity and the Seebeck coefficient as shown in Figure 2b are obtained. The peak in Δ indicates, that the effective properties of the individual networks differ from < > on average by 175%. In this case < > is an unsuitable approximation and therefore the GEMT should not be used to describe the effective properties of the individual networks. The position of the maximum is slightly shifted from v c towards the side with a smaller < >. Furthermore, the thermal conductivity does not experience such a peak in Δ , due to the low ratio II / I = 2 in our example resulting in the absence of percolation effects.
To illustrate this fact, we introduce the criterion that GEMT adequately predicts for all practical purposes the macroscopic material's property of a network if the property lies within ±5% of the average < >. Based on this definition, we calculated the prediction rate of the GEMT as a function of the ratio II / I of the electrical conductivities of the two materials and the volume fraction v II with 200 randomly cast networks per point. Figure 3 shows the results, alongside the position of the maximum of the mean relative deviation v Δ ,max . Composite materials with ratios II / I < 10 3 can be predicted by the GEMT with a rate of 100% regardless of the composition. Above this threshold, the prediction rate drastically drops to 0% close to the percolation threshold (v II = 0.5). In this region, the effective properties of the network start to differ significantly from < > leading to the peak in Δ in Figure 2b.
We define the interval [v min , v max ] in which the prediction rate falls below 50%. Within these boundaries, the GEMT becomes inapplicable to yield meaningful information about the transport properties of an individual sample, since only the properties of less than half of the networks are predicted correctly. The interval boundaries v min ≈ 0.47 and v max ≈ 0.53 for a 2D square lattice stay constant for the ratio II / I once it exceeds the critical threshold the material system becomes binary and percolation effects dominate.
To understand the nature of the drop in the prediction rate and the peak in Δ , it is helpful to look at the distribution of the material properties of the individual networks close to the percolation threshold. Figure 4 shows the histograms of the effective electri- cal conductivities of 1000 networks for different representative volume fractions v II of the two-phase composite of Table 1. The networks with v II = 0.45 are well below v c , therefore it is guaranteed that only phase I percolates through the network and TM I will be responsible for the more dominant transport mechanism for effective electrical conductivity. The distribution of electrical conductivities of all networks is binominal with a very small standard deviation.
For v II = 0.48 TM I continues to be more dominant, however the standard deviation of the −distribution increases due to the increasing size of phase II clusters in the network. Larger clusters result in a system more sensitive to small changes in the local transport coefficients. This broadening causes the prediction rate to drop, however, the probability that a network has an effective electrical conductivity = < > remains non-zero. Additionally, < > steadily decreases with increasing v II due to the rising number of low conductive bonds in the network. The -distributions for v II = 0.55, v II = 0.52 can be explained analogously with phase I and phase II being swapped.
Compositions with volume fractions equal or very close to the percolation threshold (v II = 0.5 and v II = 0.507) experience a split in their distribution of into two sub-distributions. This split emerges for values of v II in which both phases have a significantly high probability of percolating through the network. Hence, one must carefully distinguish between several qualitatively different basic constellations that can occur in a given sample: in case A only phase I forms a percolation path and transport mechanism TM I will be more dominant. Analogously in case B only phase II percolates resulting in a dominance of TM II . Case C occurs when both phases experience percolation through the network forming a parallel connection of two clusters with TM I and TM II respectively. In this condition, the current will primarily flow through phase I due to its being by orders of magnitude higher conductivity. Therefore, here TM I will be the more dominant transport mechanism. Finally, case D occurs when no phase forms a percolative path resulting in a series connection of both phases in which TM II is more dominant since it is the phase with significantly lower conductivity. At the percolation threshold, all constellations occur equally likely with a probability of 25%. As a result, the effective electrical conductivity is predominantly influenced by either TM I or TM II and never by a mixture of both. This binarity leads to the split in the distribution of for v II = 0.5, v II = 0.504, v II = 0.506, and v II = 0.508 as seen in Figure 4. In this case, the average value < > is a bad approximation of the effective electrical conductivity because it naturally falls inside the gap of the -distribution where the probability for a network to have = < > is zero. The prediction rate of the GEMT also drops down to almost 0% and Δ increases drastically for split -distributions. Therefore, it is not possible to assign a single value for compositions with volume fractions close to the percolation threshold without considering the specific morphology and percolation paths of the composites or system with multiple transport mechanisms. If II / I < 10 3 , the two sub-distributions in thedistribution in the histogram are close enough and they overlap. In this case, the network does not have a binary characteristic and percolation effects are less pronounced.
Analogous to the electrical conductivity, the thermal conductivity and the Seebeck coefficient experience the same binary behavior, since they are all linear effects being governed by the same fundamental behavior.
Fluctuation effects near the percolation threshold have therefore strong influences on the effective ZT value of nearpercolation threshold composites. Exemplarily, in Figure 5a we selected four 100 × 100 networks with the material parameters listed in Table 2.
These parameters were explicitly chosen to enlarge percolation effects. All networks have the same composition of v II = v I = 0.5, however, the networks represent each one of the previously discussed cases. The bonds are colored blue (Phase I) and red (Phase II). The bonds with the darker shade belong to a cluster connected to the input node at x = 0 and the lighter bonds have no connection to the input.
Being of equal composition the EMT predicts for all four networks the same ZT value of 0.0028. However, since for all material properties II / I >> 10 3 the composite is a binary system in which the percolation effects dominate. The true ZTs of the networks, therefore, vary over several orders of magnitudes depending on which of the percolation cases is present. In the cases A and B, charge transport occurs through the only phase forming a percolation path resulting in ZT A ≈ ZT I and ZT B ≈ ZT II . These strong variations in the possible ZT values of a given sample can, therefore, be utilized to identify, in a given production process the subset of samples with the most desired transport properties.
While for a given realization, the predictability of the transport behavior is low, the yield of the subset of samples with the desired behavior is given by the distribution functions shown in Figure 4 and is sizable.
The influence of the non-percolating phase is marginal. In case C the two phases form a parallel connection in which phase I with  the higher electrical conductivity is dominating the charge transport resulting in ZT C ≈ ZT I . Analogously in case D the network has no percolation path and consists of a series connection of the phases. Here the lower electrically conductive phase II limits the charge transport, leading to ZT D ≈ ZT II . In no constellation, the ZT value calculated by the GEMT is thereby accurate. For the exemplary networks in Figure 5a the electrical conductivity is the key material property because it has the highest II / I ratio and is the most binary property.
To show that the ratio of the electrical conductivities is the significant material property that determines the macroscopic properties we can simply vary I over several orders of magnitudes while adjusting S I , keeping ZT I , I, and all properties of phase II constant. Figure 5b displays the effective ZTs of all cases, the ZT via GEMT as well as the ZTs of the phases versus II / I . These results show conclusively the arguments presented in this work. For II << I we have the previously discussed binary system of Figure 5a. When 10 −3 < II / I < 10 3 the ZT values of the four networks fall together at the ZT calculated by the GEMT. The system seizes to be binary and all percolation effects disappear. The effective material properties become the average of the properties of both phases. When II exceeds I by orders of magnitudes, the networks have again a binary character. In cases A and B, the ZT-dominating phases have not changed since they have only one phase percolating. However, in cases C and D the phases dominating the ZT have switched. Phase II having the higher electrical conductivity now dominates the charge transport in the parallel-connected phases and phase I exhibits the lower conductive charge transport in the series connection. This proves again the emergence of percolation effects in a binary system and the influence of those on the effective ZT value of a composite depending on the specific percolation constellation. Under these conditions, the GEMT fails to predict the correct effective material properties.

Conclusion
The need to use composite materials in printed thermoelectrics is usually accompanied by a significant reduction in the thermoelectric figure of merit if compared to the performance of the better of the two components. Here we demonstrated that one can exploit the large sample-to-sample fluctuations that occur in the vicinity of a percolation threshold to significantly improve the composite's thermoelectric figure of merit for a sizable subset of the produced samples. The effect is particularly relevant for finite-size systems. Calculating the effective thermoelectric properties of a composite material close to the percolation threshold by use of RRNs and expanded nodal analysis proved to be advantageous compared with the established GEMT. The conservation of information about phase morphology, current, and heat flow paths allows the study of morphological influences on the macroscopic properties from cluster formation and percolation.
For II / I > 10 3 percolation effects dominate in composite materials and influence the effective thermoelectric transport parameters for volume fractions close to the percolation threshold. Rather than treating the composite as an effective material, the RRN model allows an analysis of charge transport in the different phases. Depending on the percolation paths in the formed network the effective material properties and the figure of merit are dominated by either the only phase percolating, the phase with the highest conductivity within parallel-connected phases, or the phase with the lowest conductivity in the case of phases connected in series. We showed that the established GEMT cannot model percolation effects and therefore results in an incorrect prediction of the effective properties and its invalidity for volume fractions around the percolation threshold. The reduced predictability of the EMT has its roots in large fluctuations of individual samples that offer a significantly enhanced spectrum of allowed electric and thermoelectric transport coefficients.