Effects of Disorder on Thermoelectric Properties of Semiconducting Polymers

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 an increase in electrical conductivity while leaving the Seebeck unchanged. On the other hand, energetic disorder leads to increased energy gaps between sites, hindering transport. This adversely affects conductivity while only slightly increasing Seebeck and results in lower TE power factors. Furthermore, positional correlation primarily affects conductivity, while correlation in site energies 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 and positional correlation, while increasing positional disorder can lead to higher TE power factors.


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 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' 29 . 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 γ −1 of the localized electron wave function, or the localization length, is typically 1 < aγ < 5 41 . The vdW and dipole-dipole interactions cause variation in the electrostatic environment 27 ; furthermore, the dopant molecules Coulombically interact with the carriers localized to a site, thus broadening the density of states (DOS) 42 ; this is called 'energetic disorder' and shown on the left of Fig. 1. The site energies are described by a Gaussian distribution of width Γ E 27,29 , and the DOS is given as where Γ E accounts for the degree of energetic disorder in the structure. Positional disorder is modeled as random variations in the wave function overlap parameter γ between two sites, as depicted schematically in Fig. 1. Hence, γ ij = γ i + γ j , where γ i and γ j are the site specific contributions obtained from a Gaussian distribution of width Σ ij , where the width Σ ij accounts for the variation in the electronic wave function coupling due to variation of both the intersite distance and mutual orientation of molecules 27 . www.nature.com/scientificreports www.nature.com/scientificreports/ 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 43 , molecular structure, and inter-molecular packing of the material 44 . Carriers hop from one localized state to another through three possibilities: 1) the electron hops to another state of equal energy by a tunneling process, 2) it hops to a lower energy site while the difference in energy is compensated by the emission of a phonon, or 3) 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 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 33 . The initial site occupation is given by the equilibrium Fermi-Dirac distribution The general hopping rate from site i to site j is given as 45 where ρ FCTW (ΔE) is a function that depends on the Franck-Condon factors, ω 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 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, ΔE ij = E j − E i − eFΔ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′ can take place, the function ρ FCTW becomes a Dirac delta function and we obtain 36 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 , which then yields 46 www.nature.com/scientificreports www.nature.com/scientificreports/ where β is the coupling constant between electron-phonon coupling constants and wave function overlap, and . For hops upwards in energy (E j > E i ) by absorption of a phonon, it is − 1 2 and for downward hops with the emission of a phonon, it is + 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 47 , where v 0 is the attempt-to-escape frequency and Δ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 effect of polaron transport, and analyze its impact 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.
, where E 0 is the reorganization energy in Eq. 3, and further simplifying we get which is the Marcus rate equation 48 . It is important to note that we can also obtain the Miller-Abrahams rate by taking the limit E 0 → 0 in the Marcus rate. The non-linear PME is solved using these rates on a 35 × 25 × 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 a 3 . 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 = 10 6 Vm −1 , attempt to jump frequency v 0 = 10 12 s −1 , energy distribution width Γ E = 3k B T, overlap γ = 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 30 (here the x direction) where e is the electron charge, and then for conductivity we take σ = J/F. The Seebeck coefficient (or thermopower) is the voltage built up in response to an applied temperature gradient, given by = − | Δ Δ = S V 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 49,50 , and the Seebeck coefficient can be calculated as the average carrier energy where E T is the average transport energy, calculated from 37 where the brackets 〈 . 〉 denote an average over the sites. The Lorenz number is related to the open-circuit electronic thermal conductivity 51 through L = κ o /σT−S 2 and thus can be analogously obtained from the variance of the excess energy 51 www.nature.com/scientificreports www.nature.com/scientificreports/

Results and Discussion
First, we consider the impact of doping and overlap between neighboring sites on TE properties. Figure 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 μVK −1 , but increasing the carrier concentration causes it to decrease to a few tens of μ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 × 10 20 cm −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 → 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 (Γ 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 Δ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 52 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 www.nature.com/scientificreports www.nature.com/scientificreports/ 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 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 53 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 (Σ 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 55 . It has been shown that energy correlation leads to Poole-Frenkel field dependence of mobility over a wide range of fields 55,56 . 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 37 . 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 [57][58][59] 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 ener- 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 www.nature.com/scientificreports www.nature.com/scientificreports/ autocorrelation function | | = ( ) 2 S F C . The autocorrelated distribution is obtained by multiplying a random phase (e iφ ), having angle φ uniformly distributed between (0, 2π), with the square root of spectral density and then taking the inverse Fourier transform 8 i 1

F S
The same procedure is applied to obtain a spatially-correlated distribution of overlap parameters γ. This method allows us to vary the correlation length independently of the variance. Figure 4a shows the Seebeck coefficient vs. electrical conductivity for uncorrelated and correlated energy distribution with different correlation lengths and fixed Γ E = 3 and 6 k B T. 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. 25 (Fig. 5b), and P3HT sample from Ref. 60 (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 ~ 0.25 obtained from Ref. 25 , is a close fit to the γ = 0.1 and Γ 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 γ = 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 = π 2 /3(k B /e) 2 = 2.44 × 10 −8 WΩ K −2 61 decreasing to the value L 0 = 2(k B /e) 2 = 1.49 × 10 −8 WΩ K −2 found in single-parabolic-band materials when acoustic phonon scattering is dominant 62 . 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ΩK −2 and S in μVK −1 ) is a good approximation and holds across all known band semiconductors 63 . 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 Γ E is increased from 3 to 5 k B T. Experimental data also confirm L is significantly larger than L 0 in PEDOT:Tos 25 . In the limit when electronic thermal conductivity dominates (κ e > κ 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.

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 www.nature.com/scientificreports www.nature.com/scientificreports/ 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.

Methods
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 , which can be written in terms of the in-and out-scattering as . We arrange the 35 × 25 × 25 array of Δp i 's into a column vector 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 64 , 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.
PDPP4T sample preparation and characterization. PDPP4T (molecular weight M n = 72.8 kDa and dispersity = 4.4) was synthesized according to previously reported procedures 65 . 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 × 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 ± 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 Figure 5. Comparison of our model to experimental data from (a) PDPP4T samples from our measurements, (b) PEDOT:Tos from Ref. 25 , and (c) P3HT sample from Ref. 60 , showing good agreement across multiple data sets. We plot two curves (solid and dashed lines) on the top and botton of the experimental data to show that the values would fall between these two extremes and account for possible error bars.