Effects of Disorder on Thermoelectric Properties of Semiconducting Polymers

TL;DR

This study models how different types of disorder in semiconducting polymers influence their thermoelectric properties, revealing that positional disorder can enhance performance while energetic disorder hampers it.

Contribution

It introduces a Gaussian disorder model combined with Pauli's master equation to analyze disorder effects on thermoelectric properties of polymers, providing new insights into optimizing TE performance.

Findings

Positional disorder improves electrical conductivity and Seebeck coefficient.

Energetic disorder hinders charge transport and reduces TE efficiency.

Correlation in positional disorder affects conductivity, but energy correlation does not influence TE properties.

Abstract

Organic materials have attracted recent interest as thermoelectric (TE) converters due to their low cost and ease of fabrication. We examine the effects of disorder on the TE properties of semiconducting polymers based on the Gaussian disorder model (GDM) for site energies while employing Pauli's master equation approach to model hopping between localized sites. Our model is in good agreement with experimental results and a useful tool to study hopping transport. We show that stronger overlap between sites can improve the electrical conductivity without adversely affecting the Seebeck coefficient. We find that positional disorder aids the formation of new conduction paths with an increased probability of carriers in high energy sites leading to a simultaneous increase in electrical conductivity and Seebeck coefficient. On the other hand, energetic disorder leads to increased energy gap between sites, hindering transport, and adversely affects both conductivity and Seebeck. Furthermore, positional correlation primarily affects conductivity, while correlation in energy has no effect on TE properties of polymers. Our results also show that the Lorenz number increases with Seebeck coefficient, largely deviating from the Sommerfeld value, in agreement with experiments and in contrast to band conductors. We conclude that reducing energetic disorder while increasing positional disorder can lead to higher TE power factors.

Equations

$$g(E)=\frac{1}{\sqrt{2\pi {\Gamma }_E^2}}exp\left(-\frac{E^2}{2{\Gamma }_E^2}\right),$$

$$\frac{d{p}_i}{dt}=0=\sum _j[W_{ij}{p}_i(1-p_j)-W_{ji}{p}_j(1-p_i)],$$

$$W_{ij}=\frac{\pi }{\hslash }\sum _q\mid M_{ij,q}{\mid }^2\left[\left(N\left({\omega }_q\right)+1\right){\rho }_{FCWT}\left(\Delta E_{ij}+\hslash {\omega }_q\right)+N\left({\omega }_q\right){\rho }_{FCWT}\left(\Delta E_{ij}-\hslash {\omega }_q\right)\right],$$

$$W_{ij}=\frac{\pi }{\hslash }\sum _q\mid M_{ij,q}{\mid }^2\left[\left(N\left({\omega }_q\right)+1\right)\delta \left(\Delta E_{ij}+\hslash {\omega }_q\right)+N\left({\omega }_q\right)\delta \left(\Delta E_{ij}-\hslash {\omega }_q\right)\right].$$

$$W_{ij}={\beta }^2{\gamma }_{ij}^2\left[N\left(\Delta E_{ij}\right)+\frac{1}{2}\pm \frac{1}{2}\right]D_{ph}\left(\Delta E_{ij}\right)/\Delta E_{ij}$$

$$W_{ij}=v_{0}exp(-2{\gamma }_{ij}{R}_{ij})\left[N\left(\Delta E_{ij}\right)+\frac{1}{2}\pm \frac{1}{2}\right]$$

$$W_{ij}=v_0\sqrt{\frac{1}{4\pi E_0{k}_{B}T}}exp\left(-2{\gamma }_{ij}{R}_{ij}\right)exp\left[-\frac{{\left(-\Delta E_{ij}+E_0\right)}^2}{4E_0{k}_{b}T}\right],$$

$$J=\frac{e}{a^{3}N}\sum _{i,j}{W}_{ij}{p}_i(1-p_j)R_{ij,x}$$

$$S=\frac{(E_F-E_T)}{eT}$$

$$E_T=<E_i>=\frac{{\sum }_{i,j}{E}_i{W}_{ij}{p}_i(1-p_j)R_{ij,x}}{{\sum }_{i,j}{W}_{ij}{p}_i(1-p_j)R_{ij,x}},$$

$${\kappa }_o=\frac{{\sum }_{i,j}(E_i-E_F{)}^2{W}_{ij}{p}_i(1-p_j)R_{ij,x}}{e{T}^2}$$

$$L=⟨{\left(\frac{E_i-E_F}{eT}\right)}^2⟩-{⟨\left(\frac{E_i-E_F}{eT}\right)⟩}^2.$$

$$\text{E}={\mathcal{F}}^{-1}\left(|\mathcal{S}|e^{i\varphi }\right).$$

Full text

Thermoelectric, polymers, disorder

I Introduction

The increased interest in electricity generation from waste heat and natural sources BellSci08 , as well as large-scale Peltier coolers SharpIEEE06 , driven by the need for alternative energy sources, has led to extensive research on more cost-effective and energy-efficient TE materials. TE conversion efficiency of a given material is governed by its dimensionless figure of merit $ZT=S^{2}\sigma T\left/\kappa\right.$, where S is the Seebeck coefficient, $\sigma$ is the electrical conductivity, T is the absolute temperature, and $\kappa=\kappa_{l}+\kappa_{e}$ is the total thermal conductivity composed of electronic ($\kappa_{e}$) and lattice contributions ($\kappa_{l}$). The electronic component of thermal conductivity is related to the electrical conductivity through the Wiedemann-Franz law $\kappa_{e}=L\sigma T$ Franz53 , where L is the Lorenz number. Most doped semiconductors have a high $\sigma$ and a moderate Seebeck coefficient, but ZT is still limited by their high lattice thermal conductivity. There has been much research effort to reduce the thermal conductivity by nanostructuring GonzalezRev13 , embedding nanoparticles in the semiconductor KimPRL06 and surface roughness HochbaumNAT08 ; MartinPRL09 ; MaurerAPL15 , which increase scattering and limit phonon transport.

Organic materials offer several advantages as they have an inherently low thermal conductivity on account of their disordered structure and thus do not require further processing such as nanostructuring to reduce it. Thermal conductivity in $\pi$-conjugated polymers such as poly(3-hexylthiophene) (P3HT), polypyrrole (PPY), polycarbazoles, polythiophenes, polyaniline (PANI) and poly(3,4-ethylenedioxythiophene) (PEDOT)GonzalezRev13 ; RothJCP91 ; ZhangSM10 ; DubeyJPP11 is at least an order of magnitude lower than inorganic compounds, with lattice thermal conductivity typically $<$1 Wm*-1K-1* KimNMat13 ; SinghIEEE01 ; LiuACS15 . They can be chemically oxidized or reduced (doped) to induce high levels of carrier density and thus improve their electrical conductivity HeegerPRL78 ; KimNMat13 . They have many additional advantages, such as low cost due to inexpensive fabrication methods and large-area production, solubility in common solvents, and solution processing RogersAdvMat99 . Thin films can be easily manufactured using techniques such as solution casting RogersAdvMat99 , vacuum evaporation HorowitzSSC89 . and printing technologies PedeMSE98 . They are also potentially disposable and have lower potential for negative environmental impact SinghIEEE01 , which makes them a very attractive choice for commercially viable TE applications CulebrasRev14 ; MingEnerEnvSc13 ; DubeyJPP11 .

Doping polymers in order to improve the electrical conductivity has the undesirable effect of significantly reducing their Seebeck coefficient to a range in the order of tens of $\mu$VK*-1*.Sun10 Therefore, a long-standing problem in TEs has been to effectively decouple conductivity from the Seebeck coefficient and control them independently. Nonetheless, doping polymer blends with a minor additive component can result in a simultaneous increase of conductivity and Seebeck coefficient Sun10 . Also, polymer-CNT composites and organic-inorganic composite materials have been shown to possess higher ZT values due to increased electrical conductivity of these materials YuACS11 ; GonzalezRev13 ; TongaAMI17 . Further improvements by mixing with graphite/graphene, carbon nanotubes, or inorganic TE nanoparticles have also been observed.CulebrasRev14 Doping PEDOT with counter ions has shown to increase the TE efficiency up to two orders of magnitude depending on the counter-ion size CulebrasJChem14 . Recently, PEDOT:Tos with enhanced Seebeck coefficient of 210 $\mu$V, electrical conductivity of 70 Scm*-1* with a resulting ZT$\sim$0.25 has been reported BubnovaNmat11 , making them a viable alternative to inorganic materials.

Polymer systems do not possess the continuous order found in their inorganic counterparts; they are inherently disordered and charge transport can be described as a hopping process PautmeierSynMat90 ; BasslerJCP91 ; TesslerAdvMat09 . Simulation is a powerful tool to understand transport in these materials and investigate the effect of disorder on the TE performance of these materials. One of the most prominent phenomenological models is based on variable-range hopping (VRH) of electrons between the polymer chains and called the Gaussian Disorder Model (GDM) BasslerJCP91 . This model has been widely used to study charge transport in polymers and polymer based devices BasslerPSS93 ; PasveerPRL05 ; GrovesJCP08 ; FishchukPRB07 ; VandelholstCh2 . More recently, charge transport in organic systems have also been analyzed by combining ab initio calculations with classical molecular dynamics (MD) and kinetic Monte Carlo techniques KirkpatrickPRL07 ; NelsonJCP10 ; WangNL09 . In the work by Mendels and Tessler TesslerPCL14 , positional disorder was implemented as spatial variations without including the variation in orbital overlap. However, the TE performance of these materials and the combined effects of disorder and correlation on TE transport, especially within a hopping model, has not been fully explored. Also, there have been few studies to determine the Lorenz number and its relationship with Seebeck coefficient in semiconducting polymers.

We investigate the TE properties of disordered organic semiconductors employing a model based on the GDM and use Pauli’s master equation (PME) approach to calculate site occupational probabilities. Most studies based on the GDM and its variants compute the hopping rate between adjacent sites using the Miller-Abrahams rate equation. We present a comparison of the Miller-Abrahams hopping rate with the Marcus hopping rate, which considers the additional energy penalty to hopping due to polaronic binding. We explore the effect of various manifestations of disorder, including positional disorder, energetic disorder, as well as correlation in both energy and wave-function overlap distributions, on the electrical conductivity, Seebeck coefficient and Lorenz number. We find that the overlap has an enormous impact on electrical conductivity whereas spatial variations have negligible effect. We find that positional disorder aids the formation of new conduction paths with an increased probability of carriers in high energy sites leading to increase in electrical conductivity. In contrast, energetic disorder leads to increased energy gap between sites hindering transport and adversely affects conductivity, however, the increased energy gap also leads to a lower average site energy and a small increase in Seebeck coefficient. Consequently, positional correlation negatively affects conductivity, while correlation in energy has no effect on TE properties of polymers.

Experimental studies on iodine-doped polyacetylene MermilliodJPhys80 and PEDOT:PSS LiuACS15 have shown that the Wiedemann-Franz law holds and the Lorenz number is close to the Sommerfeld value. However, Weathers et al. WeathersAdvMat15 showed the electronic contribution to thermal conductivity was higher than previously reported, consistent with a large Lorenz number, while Lu et al. NianduanJAP16 reported a large deviation from the Wiedemann-Franz law under the effect of temperature, carrier concentration, energetic disorder, and electric field. We find that Lorenz number increases with Seebeck coefficient, largely deviating from the Sommerfeld value, and it increases further with increasing positional and energetic disorder. The presence of disorder leads to inherently different transport and the design of efficient polymer TEs requires consideration of both positional and energetic disorder, as well as the Lorenz number, in addition to optimizing the doping concentration.

II Polymer Theory

Conjugated polymers are positionally disordered systems in which the polymer chains typically interact through weak van der Waals (vdW) interactions. The charge transport within a chain occurs through the covalent framework, and between chains the interactions is through the pi-pi orbitals orthogonal to the chain axis. The displacement of the states about the lattice points causes disruption in the overlap of the pi orbital wave functions termed ‘positional disorder’ BasslerPSS93 . The interactions between orbitals of adjacent segments are very weak, and the strong electron-phonon coupling in these materials can destroy the coherence between neighboring sites, which causes electrons to become localized to that region and reduces the delocalization range. The decay parameter $\gamma^{-1}$ of the localized electron wave function, or the localization length, is typically 1 $\AA$ $<$ $a\gamma^{-1}$ $<$ 5 $\AA$ BaranovskiiPRB04 . The vdW and dipole-dipole interactions cause variation in the electrostatic environment BasslerJCP91 . Furthermore, the dopant molecules Coulombically interact with the carriers localized to a site, thus broadening the density of states (DOS) ArkhipovPRB05 ; this is called ‘energetic disorder’ and shown on the left of Fig. 1. The site energies are described by a Gaussian distribution of width $\Gamma_{E}$ BasslerJCP91 ; BasslerPSS93 and the DOS is given as

[TABLE]

where $\Gamma_{E}$ accounts for the degree of energetic disorder in the structure. Positional disorder is modeled as random variations in the wave function overlap parameter $\gamma$ between two sites, as depicted schematically in Fig. 1. Hence, $\gamma_{ij}=\gamma_{i}+\gamma_{j}$, where $\gamma_{i}$ and $\gamma_{j}$ are the site specific contributions obtained from a Gaussian distribution of width $\Sigma_{ij}$, where the width $\Sigma_{ij}$ accounts for the variation in the electronic wave function coupling due to variation of both the intersite distance and mutual orientation of molecules BasslerJCP91 .

II.1 Transport Model

In inorganic semiconductors, ordered crystal lattice with relatively weak electron-phonon coupling leads to a band transport where the interaction between electrons and lattice vibrations (phonons) can be described by perturbation theory. In contrast, polymers do not have long-range periodicity of the atomic structure; electrons are localized and transport is described as a hopping process. The hopping transport process is dependent on temperature LuPCCP16 . molecular structure, and inter-molecular packing of the material KiessSpringer92 . Carriers hop from one localized state to another through three possibilities: the electron hops to another state of equal energy by a tunneling process, it hops to a lower energy site while the difference in energy is compensated by the emission of a phonon, or it hops to a state of higher energy and the additional energy required is provided by absorbing a phonon, as illustrated by the green arrows in Fig. 1.

Our model describes the probability that a site is occupied by an electron in terms of Pauli’s master equation (PME), which is a differential equation that describes the time rate of change of each site occupation probability due to electrons hopping into and out of it. In the steady-state, the time rate of change of occupation probability will go to zero and PME is given as a sum over all possible transitions into and out of a site

[TABLE]

where $p_{i}$ is the occupation probability of a site $i$ and $W_{ij}$ is the hopping transition rate from site $i$ to $j$, summed over the neighboring sites $j$. The PME is solved for the site occupations using a non-linear iterative solver as described in Methods, after which relevant quantities like mobility and current can be calculated.VandelholstCh2 The initial site occupation is given by the equilibrium Fermi-Dirac distribution $p_{i}^{0}=\left[\exp\left(\frac{E_{i}-E_{F}}{k_{B}T}\right)+1\right]^{-1}$.

The general hopping rate from site $i$ to site $j$ is given as FornariJCP15

[TABLE]

where $\rho_{FCTW}\left(\Delta E\right)$ is a function that depends on the Franck-Condon factors, $\omega_{q}$is the energy of the phonon mode $q$, $N_{q}$ is the number of phonons in that mode, given by the Bose-Einstein distribution $N(\omega_{q})=\left[\exp\left(\frac{\hbar\omega_{q}}{k_{B}T}\right)-1\right]^{-1}$ where $T$ is the absolute temperature. The $M_{ij,q}$ is the phonon-electron coupling constant between sites $i$ and $j$ due to phonon mode $q$, $\Delta E_{ij}=E_{j}-E_{i}-eF\Delta R_{ij,x}$ where $E_{i}$ and $E_{j}$ are the energies of sites $i$ and $j$ and $F$ is the externally applied electric field. In the limit where there are no phonons with different equilibrium positions in sites $i$ and $j$, such that only transitions $q=q^{\prime}$ can take place, the function $\rho_{FCTW}$ becomes a Dirac delta function and we obtain WangNL09

[TABLE]

Calculation of these rates can be computationally challenging as it requires that we first calculate the electronic wave functions, phonon modes, and the electron-phonon coupling constants. We can simplify them further by approximating $M_{i}$ to be proportional to the overlap of the wave functions $\gamma_{ij}=\int d^{3}r\mid\psi_{i}(r)\mid\dot{\mid}\psi_{j}(r)\mid$, which then yieldsMladenovicAFM15

[TABLE]

where $\beta$ is the coupling constant between electron-phonon coupling constants and wave function overlap, and $D_{ph}(E)$ is the phonon DOS normalized such that $\int_{0}^{\infty}D_{ph}(E)dE=1$. For hops upwards in energy ($E_{j}>E_{i}$) by absorption of a phonon, it is $-\frac{1}{2}$ and for downward hops with the emission of a phonon, it is $+\frac{1}{2}$ in the rate equation. Further simplification can be made by assuming that the wave function overlap decays exponentially with distance, and if we ignore the energy dependance we get the Miller-Abrahams rate equation,MillerPRev60

[TABLE]

where $v_{0}$ is the attempt-to-escape frequency and $\Delta R_{ij}$ is the distance between the sites.

The Miller-Abrahams rate equation considers only bare charge transport. Since the phonon-electron coupling is strong in organic polymers, it is important to consider the affect of polaron transport, and analyze its effect on TE properties. As a polaron moves through different states, there is deformation of the molecule as the polaron arrives and leaves, and the energy associated with the relaxation of the molecule upon charge transfer is called the binding or reorganization energy. Using $\rho_{FCTW}\left(\Delta E\right)=\sqrt{\frac{1}{4\pi E_{0}k_{B}T}}exp\left[-\frac{\left(\Delta E+E_{0}\right)^{2}}{4E_{0}k_{B}T}\right]$,FornariJCP15 where $E_{0}$ is the reorganization energy in Eq. 3, and further simplifying we get

[TABLE]

which is the Marcus rate equation.MarcusRPC64 It is important to note that we can also obtain the Miller-Abrahams rate by taking the limit $E_{0}\to 0$ in the Marcus rate.

The non-linear PME is solved using these rates on a 35$\times$25$\times$25 lattice of sites with an average distance between adjacent lattice points $a$=1 nm. The number of ‘neighbors’depends on the hopping distance and lattice described; we have a cubic lattice and consider up to the third-nearest neighbor, which implies hopping to the nearest 26 sites and maximum hopping distance of $\sqrt{3a}$. The electronic wavefunctions are localized so their overlap has an exponential decay with distance; hence the probability of hopping to neighbors further than the third-nearest neighbor is very small and does not contribute as significantly to transport. The $(1-p_{i/j})$ factor accounts for the exclusion principle requiring that only one carrier is occupying a particular site, due to the high Coulomb penalty for the presence of two charges on one site. All simulations are run under low-field conditions with field $F$=106 Vm*-1*, attempt to jump frequency $v_{0}$=1012 s*-1*, energy distribution width $\Gamma_{E}$=3kBT, overlap $\gamma$=3 and $T$=300 K unless stated otherwise. We obtain the current density $J$ by a summation over all the carriers in the direction of the applied field PasveerPRL05 (here the x direction)

[TABLE]

where $e$ is the electron charge, and then for conductivity we take $\sigma=J/F$.

The Seebeck coefficient (or thermopower) is the voltage built up in response to an applied temperature gradient, given by $S=-\frac{\Delta V}{\Delta T}|_{I=0}$, as carriers that respond to an electric field can also be elicited by a temperature gradient. While each carrier carries a charge $e$, it also carries an ‘excess’ energy $E-E_{F}$ ZvyaginPSSB80 ; KangNMat17 , and the Seebeck coefficient can be calculated as the average carrier energy

[TABLE]

where $E_{T}$ is the average transport energy, calculated from TesslerPCL14

[TABLE]

where the brackets $<.>$ denote an average over the sites. The Lorenz number is related to the open-circuit electronic thermal conductivity WangJAP18

[TABLE]

through $L=\kappa_{o}/\sigma T-S^{2}$ and thus can be analogously obtained from the variance of the excess energy WangJAP18

[TABLE]

III Results and Discussion

First, we consider the impact of doping and overlap between neighboring sites on TE properties. Fig. 2a shows the dependence of electrical conductivity and Seebeck coefficient on the carrier density. We increase the carrier density by moving the Fermi level $E_{F}$ closer to the center of the Gaussian energy distribution, analogous to electrochemical doping of polymers. We can clearly see the inverse relation that exists between these parameters and the charge density, and the challenge it poses to obtaining high power factor values. At low concentration, the Seebeck is in the range of hundreds of $\mu$VK*-1*, but increasing the carrier concentration causes it to decrease to a few tens of $\mu$VK*-1*. However, the conductivity has an appreciable value only at high concentrations and the highest power factor is achieved at a carrier concentration of $2\times 10^{20}$ c$m^{-3}$, shown in Fig. 2b, which corresponds to 20% of sites being occupied by carriers.

We plot Seebeck coefficient vs. conductivity for different values of the overlap parameter in Fig. 2c, where each point on the curve represents the parameters computed at different carrier densities. The advantage of such a plot is that one can readily see the effect of both carrier density and overlap parameter on Seebeck and conductivity, wherein a curve bulging more towards the top right indicates higher power factor. We find that electrical conductivity is strongly dependent on the overlap parameter whereas it has negligible effect on the Seebeck coefficient, as seen in Fig. 2a. A smaller value of the overlap parameter, which implies stronger electronic orbital overlap between adjacent sites, is favorable for the hopping process, and thus increases the electrical conductivity. In Fig. 2d we compare the Seebeck and conductivity computed from Miller-Abrahams and Marcus rate equations. The curve shifts right for higher values of reorganization energy $E_{0}$, showing improved conductivity with increasing polaronic binding, while as we decrease $E_{0}$ and approach the limit $E_{0}\to 0$, it falls back to the curve obtained from Miller-Abrahams rate.

Next, we explore the effect of varying degrees of energetic and positional disorder in the system on its TE performance. A larger variation of the Gaussian energy distribution ($\Gamma_{E}$), leads to a decrease in conductivity and small increase in Seebeck (see Fig. 3a). Decreasing the width of the energy distribution, meaning a more sharply peaked DOS, lowers the spread in the site energies, leading to a smaller difference $\Delta E_{ij}$ between energies of neighboring sites. A favorable network of nearly equal energy sites thus forms, alleviating the required thermal assistance by absorption of phonons. It is widely thought, based on the work on Mahan and Sofo MahanPNAS96 that a narrower DOS leads to a higher thermoelectric figure-of-merit ZT, with a delta-function DOS being ideal.

Combining Eqns. 9 and 10 we see that the Seebeck coefficient can be viewed as the average “excess” (away from the Fermi level) entropy per carrier $S=<E_{F}-E_{i}>/eT$. Increasing energetic disorder broadens the DOS and makes it flatter but does not affect the shape of the Fermi-Dirac distribution function. Consequently, there are nearly as many states near the DOS peak as there are away from it, which pushes the average $E_{i}$ away from the center of the DOS. Reducing energetic disorder has the opposite effect: the DOS becomes more sharply peaked with many more states near the peak, resulting in an average $E_{i}$ closer to the middle of the DOS. If we fix the Fermi level $E_{F}$ and compare, then larger energetic disorder would imply smaller Seebeck and vice versa. However, if we compare while keeping the carrier concentration $n(E_{F},T)=\int{g(E)p_{i}^{0}(E,T)dE}$ constant, then the opposite trend emerges: increasing energetic disorder pushes the Fermi level away from the center of the DOS, countering the change in $E_{i}$ and thus slightly increasing the Seebeck, shown on the left axis in Fig. 3a. More quantitatively, based on the Mott formula for the Seebeck coefficient,HeremansSci08 $S=-(\pi^{2}/3)(k^{2}T/q)\frac{\partial}{\partial E}ln[\sigma(E)]|_{E=E_{F}}$ for carriers near the Fermi level, which can be expanded for hopping transport as $S\propto\frac{dln[g(E)]}{dE}+g\frac{d[\mu(E)]}{dn}$.IhnatsenkaPRB15 In the case of GDM and Gaussian DOS, the first term is dominant and becomes $-E_{F}/\Gamma_{E}^{2}$, while the second term is typically small. Overall, Fig. 3b shows that the effect of disorder on conductivity is more pronounced and reducing energetic disorder leads to higher TE power factors.

Next, we vary the amount of positional disorder in the system by varying the width of the Gaussian overlap distribution ($\Sigma_{ij}$). Positional disorder primarily affects the conductivity, while the Seebeck is mostly sensitive to the distribution of energies (DOS) and less on their relative positions (Fig. 3c). The conductivity increases with increasing positional disorder, and shifts the curves right in Fig. 3d, signifying higher power factor values. This is due to the increase in overlap of approximately half the near-neighboring sites in the system, which aids the formation of conduction paths, consequently increasing the probability of hopping into higher-energy sites. Thus, larger positional disorder but smaller energetic disorder are desirable for polymer TEs.

In the GDM, site energies are distributed independently with no correlations occurring over any length scale. Nonetheless, spatial fluctuations and corresponding correlation of energy arising from charge dipole interactions and molecular density fluctuations should affect transport NovikovPRL98 . It has been shown that energy correlation leads to Poole-Frenkel field dependence of mobility over a wide range of fields GartsteinCPL95 ; NovikovPRL98 . Spatial correlation in the energetic landscape, modeled by averaging energy over neighboring sites, was also shown to lower the transport energy and decrease the Seebeck coefficient TesslerPCL14 . However, the impact of correlation length has not been firmly established in the context of TE properties, nor has the effect of correlation on positional disorder through the overlap parameter been studied.

To further explore the effect of long-range correlation on TE parameters, we use an inverse fast Fourier transform (IFFT) method WuTB00 ; BuranIEEE09 ; RamayyaPRB12 to generate autocorrelated distributions of energy and overlap parameter with a specific correlation length. We start from the standard exponential form of the autocorrelation function of the site energies $\mathcal{C}(R_{ij})=<E(R_{i})E(R_{j})>=\Gamma_{E}^{2}\exp\left(-\sqrt{2}R_{ij}/L_{corr}\right)$ where $R_{i/j}$ are the positions of the $i/j$’th sites, $R_{ij}$ is the distance between the two sites, and $L_{corr}$ is the correlation length. The spectral density is Fourier transform of the autocorrelation function $|\mathcal{S}|^{2}=\mathcal{F}(\mathcal{C})$. The autocorrelated distribution is obtained by multiplying a random phase ($e^{i\phi}$), having angle $\phi$ uniformly distributed between (0, 2$\pi$), with the square root of spectral density and then taking the inverse Fourier transform MaurerAPL15

[TABLE]

The same procedure is applied to obtain a spatially-correlated distribution of overlap parameters $\gamma$. This method allows us to vary the correlation length independently of the variance. Fig. 4a shows the Seebeck coefficient vs. electrical conductivity for uncorrelated and correlated energy distribution with different correlation lengths and fixed $\Gamma_{E}$=3 and 6 kBT. The same is shown in Fig. 4b for correlation in overlap distribution. We find that energy correlation has virtually no effect on the TE parameters, while the curve shifts left with increasing correlation length for the overlap distribution. We conclude that smaller correlation between sites leads to better conductivity but the Seebeck remains largely unaffected, thus effectively decoupling these parameters.

We compare our model to experimental data from several measurements, including three PDPP4T samples (Fig. 5a) that we chemically doped with iodine and measured as described in Methods, PEDOT:Tos from Ref. BubnovaNmat11, (Fig. 5b), and P3HT sample from Ref. JungSciRep17, (Fig. 5c). Our results are in good agreement with the data and the fit is obtained by varying relevant parameters including overlap, average distance between adjacent lattice points, energetic, and positional disorder. The data for PEDOT:Tos reported to have a ZT$\sim$0.25 obtained from Ref. BubnovaNmat11, , is a close fit to the $\gamma$=0.1 and $\Gamma_{E}$=2.25 $k_{B}T$ curve, implying stronger electronic orbital overlap between adjacent sites and small energetic disorder, which explains the exceptional conductivity observed in these samples beyond what is obtained from $\gamma$=3, a value commonly used in calculations.

Lastly, we turn our attention to the Lorenz number. In most materials, the Lorenz number ranges from a value close to the Sommerfeld value found in metals and degenerate semiconductors $L_{0}=\pi^{2}/3(k_{B}/e)^{2}=2.44\times 10^{-8}$ W$\Omega$ K*-2* Sommerfeld27 decreasing to the value $L_{0}=2(k_{B}/e)^{2}=1.49\times 10^{-8}$ W$\Omega$ K*-2* found in single-parabolic-band materials when acoustic phonon scattering is dominant ThesbergPRB17 . It has been shown that a first-order correction to the degenerate limit $L=1.45+exp(-|S|/116)$ (where L is in 10*-8* W$\Omega$K*-2* and S in $\mu$VK*-1*) is a good approximation and holds across all known band semiconductors KimAPL15 . In contrast, we see an opposite trend in organic semiconductors, where $L$ increases with Seebeck coefficient as shown in Fig. 6. Increasing the overlap, positional disorder and polaronic binding energy increases the value of Lorenz number further, but the largest impact is seen with energetic disorder when $\Gamma_{E}$ is increased from 3 to 5 kBT. Experimental data also confirm L is significantly larger than $L_{0}$ in PEDOT:Tos  BubnovaNmat11 .In the limit when electronic thermal conductivity dominates ($\kappa_{e}>\kappa_{l}$), the ZT goes to $ZT=S^{2}/L$; therefore, a simultaneous increase in S and L could have a negative impact on ZT and hence design of effective TE materials with polymers requires consideration of the Lorentz number as well, carefully balancing the roles of disorder and correlation.

IV Conclusion

Polymers, with their inherently low thermal conductivity and low cost of manufacturing, are a promising choice for TE applications. However, for commercial success the figure of merit of these materials needs to be improved further. We study the effects of disorder and correlation on their TE properties using a hopping transport model. We find that positional disorder leads to a moderate increase in the electrical conductivity whereas the Seebeck remains unaffected. Energetic disorder has an adverse affect on conductivity but leads to a moderate increase in Seebeck coefficient at lower doping concentrations. Consequently, positional correlation primarily affects conductivity, while correlating the nearby site energies has no effect on the TE properties. We conclude that controlling energetic and positional disorder would allow us to decouple conductivity and the Seebeck coefficient. Minimizing energetic disorder and correlation while increasing positional disorder results in a higher TE power factor. Our results also show that the Lorenz number increases with the Seebeck coefficient, more so with increasing disorder, in contrast to the universal trend observed across all materials exhibiting band transport. We find that numerical transport models can play a key role in predicting the optimum structural characteristics and aid the design and development of novel materials for TE applications.

V Methods

V.1 Solving the non-linear PME

We solve the non-linear PME using a standard iterative non-linear solver. First, we cast the PME as zero-finding for a system of equations $F_{i}(\mathbf{p})=\sum_{j}[W_{ij}p_{i}(1-p_{j})-W_{ji}p_{j}(1-p_{i})]=0$, which can be written in terms of the in- and out-scattering as $F_{i}(\mathbf{p})=p_{i}S_{out}(\mathbf{p})-(1-p_{i})S_{in}(\mathbf{p})$ with $S_{out}(\mathbf{p})=\sum_{j}W_{ij}(1-p_{j})$ being the out-scattering term and $S_{in}(\mathbf{p})=\sum_{j}W_{ji}p_{j}$ being the in-scattering. The expression for $F_{i}(\mathbf{p})$ is nonlinear because both in- and out-scattering terms depend on the unknown $\mathbf{p}$. Previous studies typically rearranged this equation to obtain a fixed point iteration for the $p_{i}$ such as $p_{i}^{n+1}=S_{in}(\mathbf{p}^{n})\left/[S_{in}(\mathbf{p}^{n})+S_{out}(\mathbf{p}^{n})]\right.$ with the initial condition $p_{i}^{0}$ being the equilibrium Fermi-Dirac distribution given earlier. However, a fixed point iteration can stall, resulting in poor convergence for some cases. For this reason, we follow the same fixed-point iteration procedure here for the first few iterations and then use the resulting estimate of $p_{i}$ as an initial guess in the next step where we numerically solve for the zero of $F_{i}(\mathbf{p})$. Rather than solving for the site occupancies $p_{i}$, we solve for their deviation away from equilibrium $\Delta p_{i}=p_{i}-p_{i}^{0}$. Combining this with Eq. 2 and simplifying for $\Delta p_{i}$, we get $F_{i}(\mathbf{p})=\Delta p_{i}S_{out}(\mathbf{p})-(1-p_{i}^{0})S_{in}(\mathbf{p})=0$. We arrange the 35$\times$25$\times$25 array of $\Delta p_{i}$’s into a column vector $\mathbf{p}$ and compute the Jacobian matrix of derivatives of $F_{i}$ with respect to $p_{j}$ as $J_{ij}=dF_{i}/dp_{j}=-W_{ij}p_{i}-W_{ji}(1-p_{i})$. Then we apply the Levenberg-Marquardt algorithm,NumRecBook as implemented in MATLAB’s fsolve function, with the known Jacobian matrix, which requires a linear solve at each iteration but typically converges in a few iterations due to its high rate of convergence. The linear solver is a preconditioned Conjugate Gradients algorithm with a banded preconditioner based on an incomplete Cholesky factorization.

V.2 PDPP4T sample preparation and characterization

PDPP4T (molecular weight Mn=72.8 kDa and dispersity Đ=4.4) was synthesized according to previously reported procedures LiACS11 . Film preparation. A solution 8 mg/mL PDPP4T in chloroform was stirred for no less than 4 h prior to dropcasting 0.23 mL of the polymer solution onto a handcut, 1.1 cm $\times$ 2.2 cm glass coverslip that was preheated to 45 *∘*C on a hot plate. This was immediately covered with a watch glass to impede the escape of chloroform vapors and slow the rate of evaporation during dropcasting. After 10 min, the sample heating element was turned off and the sample was let stand under ambient conditions for no less than 24 h to further evaporate the chloroform. Doping with iodine vapor. 50$\pm$5 mg iodine crystals were loaded into a 1 mL glass vial, and this vial was loaded into a 20 glass mL vial and sealed with a cap closed for no less than 12 h to allow iodine vapor to diffuse within the vial. Iodine doping was carried out by placing the polymer film into this 20 mL vial, immediately resealing by capping the vial, and letting stand for 24 h.

V.3 Characterization of thermoelectric properties

Samples were transferred from their iodine doping chamber to a custom-built apparatus for thermoelectric characterization in a timely fashion, since they began dedoping immediately and rapidly in the absence of iodine vapor. The details of this apparatus are reported elsewhere TongaAMI17 . In brief, a thermal gradient is established between two copper blocks by heating one of them and leaving the other at ambient temperature. The polymer film was placed onto a glass slide located on these copper blocks, and Pt wire electrodes embedded in a PTFE block were firmly secured onto the sample. A LabView program was used to interface with a digital dual input thermometer with k-type thermocouples, a Keithley 6182 nanovoltmeter, and a Keithley 2440 source meter to repeat measurements of the temperature gradient $\Delta T$, voltage gradient $\Delta V$, and I-V characteristics across the sample sequentially and every 2 min. $\Delta V$ was taken to be the average of 1000 voltage measurements from the Keithley 6182 nanovoltmeter and the Seebeck coefficient was taken to be $\Delta V/\Delta T$. Only measurements for which the standard deviation of the 1000 voltage measurements is less than 1% of their mean are considered to calculate the Seebeck coefficient and the rest discarded to ensure only reliable estimates are used to compare to our model. The conductance was taken to be the slope of the $I-V$ curve and electrical conductivity was taken to be $G\times l/A$. $l$ is the length between two platinum electrodes and $A$ is the area of the device and is taken to be $w\times t$ where $w$ is the width of the film and $t$ is the thickness obtained from profilometry after thermoelectric measurements were finished. Since the iodine doped polymer film is initially heavily doped, and spontaneously dedopes over time, these repeated measurements collect Seebeck coefficients and conductivities across a broad range of doping. We will describe this methodology in more detail in an upcoming report.

VI Data Availability.

The datasets and codes generated and/or analysed during the current study are available from the corresponding author on reasonable request.

VII Acknowledgements

The authors thank Dr. Feng Liu and Prof. Thomas P. Russel of the University of Massachusetts Amherst for providing PDPP4T for these studies. The authors also thank Profs. Frank E. Karasz and Paul M. Lahti for the thermoelectric measurement equipment and Dr. Ljiljana Korugic for useful discussions.

VIII Author Contributions

Z.A. and D. V. conceived the idea. Z. A. implemented the simulation code. M. U. performed the simulations and created the graphs. D. V. and C. J. B. designed and implemented the experimental setup. C. J. B. performed the measurements. All authors co-wrote the manuscript.

IX Competing interests

The authors declare no competing interests.