Model of the Thermoelectric Properties of Anisotropic Organic Semiconductors

A model of charge hopping transport that accounts for anisotropy of localized states and Coulomb interaction between charges is proposed. For the anisotropic localized states, the degree of orientation relates exponentially to the ratio of conductivities in parallel and perpendicular directions, while the ratio of Seebeck coefficients stays nearly unaffected. However, the ratio of Seebeck coefficients increases if Coulomb interaction is screened stronger in a direction parallel to the predominant orientation of the localized states. This implies two different physical mechanisms responsible for the anisotropy of thermoelectric properties in the hopping regime: electronic state localization for conductivities and screening for Seebeck coefficients. This provides an explanation for the recent experimental findings on tensile drawn and rubbed polymer films.


■ INTRODUCTION
Recently, rubbing and tensile drawing have been proposed as methods to enhance thermoelectric efficiency of conjugated polymers. 1−6 This type of mechanical processing makes otherwise random orientation of long polymer chains become uniaxially aligned. A high degree of orientation has been directly observed by polarized optical microscopy, 1,2,5,6 transmission electron microscopy, 2,5 and wide-angle X-ray scattering. 3 For example, Untilova et al. 1 found that poly(3-hexylthiophene) (P3HT) doped with Mo(tfdCOCF 3 ) 3 reveals conductivity along the rubbing direction (σ ∥ ) that is 2.6 times larger than conductivity in the isotropic samples (σ) of the same compound and by the absolute value σ ∥ exceeds all of the previous experimentally measured conductivities on P3HT. In contrast, conductivity σ ⊥ was by an order of magnitude smaller in the perpendicular direction: for oxidation level 11%, σ ∥ = 681 S/cm, σ ⊥ = 50 S/cm, σ = 260 S/cm. 1 Similar ratios and large absolute values along the rubbing direction have been observed for the Seebeck coefficient: S ∥ /S ⊥ = 7.1 for oxidation level 11%. In the experiment carried out by Hynynen et al., 3 P3HT doped with Mo(tfdCOCF 3 ) 3 revealed strongly increased conductivity along the drawing direction, whereas the Seebeck coefficient was surprisingly unaffected. Apart from rubbing and tensile drawing, anisotropic P3HT films with similar properties have also been fabricated by a different technique, where fiber morphology was created via epitaxial growth and temperature-gradient crystallization when organic small-molecule 1,3,5-trichlorobenzene particles were added to the solution. 7 Experimentally 1−7 measured large conductivities make anisotropic polymer films attractive as an electrode material in printed electronics, 8 while their large power factors make them attractive for application in thermoelectric generators. 9 Theoretical treatment of enhanced thermoelectric properties of the anisotropic polymer is limited to kinetic Monte-Carlo 4 and resistor network studies, 3 in which several drawbacks are faced. For example, the latter approach is limited to 2D transport with positional dependence of the tunneling rates and interaction between charges both disregarded. Kinetic Monte-Carlo modeling in ref 4 included Coulomb interaction between charges. However, only "on-site" interaction was taken into account, which neglects long-range nature of the repulsive Coulomb force between charged particles. Neglecting the longranged part of Coulomb interaction results in an incorrect ground state and inability to capture the Coulomb gap, which is a fundamental property of the disordered system of the localized states 10−12 that has been confirmed experimentally. 13,14 The long-range part is particularly important for mediums with low dielectric permittivity to which organic semiconductors belong to. Typical permittivity of the organic semiconductor is about 3. 9,15,16 Furthermore, both theories 3,4 use a simplified model of spatial anisotropy of the localized states and also use exponential density of states (DOS), which contradicts a number of studies, which have pointed out that Gaussian DOS is more accurate for disordered organic semiconductors. 16−21 Both theories predict many-fold increase of σ ∥ /σ ⊥ , in agreement with experiments. 1−7 This increase was attributed to the increase of the anisotropy degree of the localized states. However, S ∥ /S ⊥ disagreed with experimental data in refs 1 2 4,−6 where this ratio was 3−7, while the theory 4 predicted only 1.1−1.4 for the parameter range where the agreement on σ ∥ /σ ⊥ was achieved. This implies different physical mechanisms responsible for S anisotropy. Thus, the origin of enhanced thermoelectric properties of anisotropic organic semiconductors observed in recent experiments 1−6 has remained an open question.
A typical organic semiconductor contains a blend of conducting polymer molecules, dopants, and insulating host molecules. 22,23 A polymer is a long molecule that is composed of many repeating subunits (C 4 H 2 S for P3HT) that join together by covalent bonds. Due to covalent bonding, the electron wave function is delocalized along the molecule. Therefore, the whole molecule can be represented by a single localized state that is highly anisotropic in space. For an isotropic material, a blend of randomly oriented states gives zero net anisotropy, similar to magnetic dipoles in a paramagnet in the absence of a magnetic field. Rubbing or tensile drawing 1−6 or the method proposed in ref 7in the first placechange orientation and not the degree of anisotropy of the localized states (molecular structure of the polymer). One of the aims of this study is to explore how orientation affects the thermoelectric properties of anisotropic materials.
Another aim of this study is to understand the effects due to long-ranged Coulomb interaction on the thermoelectric properties of anisotropic materials, when unperturbed (singleelectron) DOS is characterized by the Gaussian distribution.
To achieve these objectives, a model of charge hopping transport is formulated that accounts for the effects of the orientation of anisotropic localized states and long-range Coulomb interactions between those states. The results that are presented in this article show that σ ∥ /σ ⊥ increases exponentially with a degree of orientation, while S ∥ /S ⊥ remains nearly constant. Maximum power factors for making an efficient power generator can be attained at an intermediate orientation degree, where less geometrical constraints apply on the current flow. Visualization of the current flow reveals a predominant path along the host orientation. For typical material parameters and at room temperature, the charge transport in both parallel and perpendicular directions occurs as 3D hopping and at a crossover between variable range hopping (VRH) and nearestneighbour hopping (NNH) regimes. The Coulomb interaction does not affect the orientational dependencies of σ and S. However, in comparison to the noninteracting theory, the ratio S ∥ /S ⊥ can increase if the Coulomb interaction is screened more strongly in one direction than another, specifically in a direction parallel to the predominant orientation of the localized states. This anisotropic screening might be a result of the larger extent of the electron wave functions along polymer backbone chains when a system behaves like metallic. Screening of the electric field inside the anisotropic organic semiconductor, which depends on morphology as well as the presence of conducting layers nearby (like metal gate electrode), explains why in some experiments 3,7 S ∥ /S ⊥ ≈ 1 but in the others 1,2,4−6 S ∥ /S ⊥ > 1, while σ ∥ /σ ⊥ > 1 in all of them. These results provide a microscopic explanation for the thermoelectric properties of anisotropic polymer films in recent experiments. 1−7

■ MODEL
The hopping conduction between localized states in a disordered system is modeled by a resistor network. 24−26 The resistance between two states i and j is 10 where the average tunneling rate accounting for wave function anisotropy is i k j j j j j j j j j y { z z z z z z z z z with γ 0 being the electron-phonon coupling parameter, (x ij , y ij , z ij ) are coordinate components of the vector connecting i and j sites, E i is the energy of the i-th state, and μ is the chemical potential.
The localized electronic states are characterized by ellipsoids of revolution with semimajor and semiminor axes, ξ ∥ and ξ ⊥ , as shown in Figure 1a. Ratio ξ ∥ /ξ ⊥ describes the degree of anisotropy. Each ellipsoid is tilted with respect to the reference frame by a random pair of angles θ and φ. The unit vector ξ′ defines the direction of ξ ∥ , as shown in Figure 1a. The Cartesian coordinates of the localized states in 2 are the distances from the origin to the cross points of the ellipsoid surface with corresponding axes Extra subscripts i and j in 2 denote the lattice sites. Note that the minor principal axes are equal for an ellipsoid of revolution, ACS Physical Chemistry Au pubs.acs.org/physchemau Article but ξ x ≠ ξ y ≠ ξ z in general. For the isotropic state, ξ ∥ = ξ ⊥ = ξ x = ξ y = ξ z , and 2 reduces to a familiar expression for the tunneling rate. 10 An average orientation of the localized states is represented by the symmetric second-order tensor, which is calculated as the dyadic product of the unit vectors Ä The degree of orientation is thus conveniently written as to make O d ∈ [0, 1]; λ 1 is the largest eigenvalue. The orientation of individual states is obtained from random distribution of the tilt angles of ξ′; see Figure 1b for three representative orientations.
A network of the localized states is put on a simple cubic lattice (3D Cartesian grid) with the unit constant l. No positional disorder is assumed.
Here, two models for energies E i are considered: a noninteracting model and a model with Coulomb interaction between charged particles. In the noninteracting model, E i = E i 0 with E i 0 randomly generated from the Gaussian distribution having the standard deviation σ DOS . In the interacting model, the energies are additionally renormalized leading to the local meanfield equations 27 where r ij is the distance between sites i and j and Q b is the positive background charge equal to the relative charge concentration in the noninteracting model. The summation runs over all of the lattice indices and takes only the shortest distance between two sites in the repeated lattice; periodic boundary conditions are applied. Coulomb interaction results in the electrons moving in the average potential generated by all other electrons. This model is known 27 to correctly reproduce the Coulomb gap and Efros−Shklovskii VRH at low temperatures. 10,12 The method to find conductivity σ and current densities I is described in ref 28. The system is assumed to be in a linear Ohmic regime. The chemical potential determines the charge density where N 0 = l −3 is the concentration of sites and f is the Fermi-Dirac distribution function. The Seebeck coefficient, or thermopower, is given by If conductivity is determined by diffusion and drift of noninteracting particles, then the Einstein relation can be applied to write (11) where ρ is DOS at the chemical potential and D is the diffusion coefficient. The Seebeck coefficient 10 thus becomes Because D depends weakly on n, S is mostly determined by ρ and its slope over energy.

■ RESULTS AND DISCUSSION
The numerical calculations are performed for a parameter set that is typical for organic semiconductors. [15][16][17][18]29,30 In particular, ξ ∥ is chosen to be equal to the lattice constant, a value large enough not to bring the system into a strong localization (insulating) regime. The anisotropy of the localized states is ξ ∥ /ξ ⊥ = 4. The lattice constant l = 1 nm. The electronphonon coupling γ 0 = 10 13 s −1 . The strength of energetic disorder is σ DOS = 0.1 eV. T = 300 K. The disorder is assumed to be only energetic and orientational; the effect of positional disorder will be commented on later. The system size for the results presented below is 20 × 20 × 20 unless otherwise stated. Averaging is performed over 100 different disorder realizations. The calculations were also performed for different sizes, ξ ∥ and ξ ⊥ , and similar results were obtained.
To understand how the orientation of the localized states affects the thermoelectric properties, let us first consider the noninteracting theory.
Conductivity, Seebeck coefficient, and power factor S 2 σ (PF) for a system of anisotropic localized states follow concentration dependence similar to that of isotropic states, 29−31 as shown in Figure 2a−c. As charge concentration (or chemical potential) increases, the effects on σ and S go in opposite directions, so the maximum of PF occurs at some intermediate n that is referred to as an optimal doping level. 31 This is a desirable value for making an efficient thermoelectric generator. 9 n is directly proportional to the oxidation level measured in the experiments. 29,30 In the absence of a predominant orientation, the network of anisotropic localized states is a system with random spatial extents of the wave functions localized on the lattice sites and it acts as if it is made of the isotropic states but randomly distributed in position. The macroscopic quantities, such as σ, are thus direction-independent. When uniaxial orientation starts to reveal, it clearly manifests itself in σ ∥ and σ ⊥ , whose ratio exponentially increases with O d , as shown in Figure 2a The transport regime can be determined from the temperature dependence of the reduced activation energy, 32 as shown in Figure 3a, which suggests that for typical values of σ DOS , [15][16][17][18]29,30 the charge transport at room temperature occurs at a crossover between VRH and NNH. For a charge carrier, the phonon energy becomes insufficient to assist hopping to the nearest localized states and hopping to the distant states becomes energetically favorable. This can be directly observed in the current density visualization, both parallel and perpendicular to the alignment direction; see Figure 3b,c. The charge flow spans uniformly over device volume and reveals a substantial degree of anisotropy that reflects the underlying orientation of the localized states and their anisotropy. If the dominant orientation is perpendicular to the current flow direction, see Figure 3c, then the electrons adjust their path further in the perpendicular direction to find a site to hop in that is closer in energy. 28 The average hopping lengths might be estimated from the length of individual current segments, which are 2.12l and 2.07l, for Figure 3b,c, respectively.
The results of the noninteracting theory show that the uniaxial orientation of the anisotropic states results in an exponential increase of the ratio σ ∥ /σ ⊥ , while S ∥ /S ⊥ remains nearly constant and the transport regime is at the crossover of VRH and NNH. Similar ratios were observed in the experiments on P3HT films in refs 3 and 7. However, other experimental data on rubbed and tensile drawn polymer films, 1,2,4−6 including P3HT, 1,2,5 revealed a simultaneous manifold increase of σ ∥ /σ ⊥ and S ∥ /S ⊥ in the highly oriented state. One of the reasons that might cause the latter observations is the effect of electron−electron interaction, which was not taken into account in the abovementioned results. Figure 4 shows how Coulomb interaction affects the concentration dependence of σ and S. For the same parameter set as in the noninteracting modeling, the concentration dependence of σ occurs more weakly. This is a result of DOS renormalization due to Coulomb repulsion between charge carriers; see the inset in Figure 4b, where the Coulomb gap 10−12 at μ is clearly seen. Note that the existence of the Coulomb gap in the DOS spectrum of disordered systems was confirmed in the electron tunneling experiments. 13,14 At high concentrations, the number of states available for conduction is reduced and thus σ becomes smaller when compared to the noninteracting case. At n/N 0 ≈ 0.015, the effect due to electron interaction on σ diminishes and reverts at lower n. The DOS shape at the Coulomb gap is close to symmetric, even though the singleparticle ρ is the exponential function of energy, which results in smaller absolute values of the Seebeck coefficient, as can already be obtained from eq 12. Orientational dependence is not affected by Coulomb interaction because the renormalized energies 8 are scalars entering the tunneling rates (the second term in 2). For the sake of visual clarity, only O d = 0 and O d = 0.97 are presented in Figure 4.   33 In (b,c), the dots mark the hopping sites with size inversely proportional to the absolute energy difference to μ; the same energetic disorder in (b,c). The gray pads are the source and drain electrodes. T is chosen to be room temperature [RT in (a)], which for a typical organic semiconductor [15][16][17][18]29,30 relates to the degree of disorder as 4k B T = σ DOS . Lattice 10 × 10 × 5 is shown. n/N 0 = 0.1.
While the interacting theory alone cannot explain the findings on S ∥ /S ⊥ in refs 1 2 4,−6, an important result is that the Coulomb interaction causes a many-fold suppression of the Seebeck coef f icient in comparison to the noninteracting theory for the same system. Thus, experimental findings 1,2,4−6 might be explained to the dominance of the Coulomb interaction in only one (perpendicular) direction. This can in turn be rationalized by the fact that in the parallel direction of the highly ordered organic semiconductor, the wave functions overlap strongly and screening of the electric field is more effective. Note that the mean-field theory cannot capture this effect because it is formulated for point-like charges that interact via unscreened (isotropic) Coulomb potential. Screening might be introduced into the theory via the dielectric constant κ by replacing e 2 by κ e 2 in 8 and further assuming that it is direction-dependent. (Examples of materials with anisotropic κ include barium titanate, black phosphorus, and nematic liquid crystals.) In the case of perfect screening, κ ∥ → ∞, and S ∥ is given by the noninteracting theory. If κ ⊥ = 1, S ⊥ is given by the result from the interacting theory. Therefore, S ∥ /S ⊥ would be the ratio between the values in the noninteracting and interacting theories and, for n/N 0 = 0.1, S ∥ /S ⊥ ≈ 2.5, as shown in Figure 4b. This estimate should be lower for the organic semiconductors, which have κ ≈ it might be argued that electron−electron interactions are suppressed due to, for example, a nearby gate electrode or another conducting layer. It is then straightforward to implement a verification of this theory: place a metal near the sample and check whether S ∥ /S ⊥ decreases or not.
Anisotropic screening in P3HT or PBTTT might be deduced from consideration of the polymer structure where the dopant ions such as F4TCNQ-are intercalated in-between the layers of π-stacked polymer chains. Hence, Coulomb potential should have a strong periodic variation perpendicular to the chain direction, whereas the polaronic charge delocalization along the chain direction should help in the screening of the electric field in that direction.
Transport regimes obtained in the noninteracting model should still be valid for the experimental conditions, 1−7 even though electron interactions are included into theory because Efros−Shklovskii VRH occurs at much lower T. 34 The crossover between Mott's VRH and NNH depends primarily on wave function localization, which can be demonstrated using Mott's hopping theory. 33 In the VRH regime, the average hopping distance is given as follows 10 (13) where ρ is DOS at the chemical potential and  (14) is the energy at the percolation threshold. At the crossover, r ≈ l. Thus, the following relation holds which implies that as the wave functions become more localized on the trapping sites, a lower temperature is required to bring the system from NNH into the VRH regime. The relation 15 was obtained for ρ being constant over at least the scale of k B T. This is not fulfilled for Gaussian distribution, but the analysis similar to ref 35 can be applied to show that this relation should also be valid for Gaussian DOS. Also, note that the Mott's law with exponent 1 4 holds for Gaussian DOS in Figure 3a. The VRH-NNH crossover depends on wave function localization but not on the degree of orientation O d . This seemingly contradicts the argument in ref 5 that different transport mechanisms exist in parallel and perpendicular directions. As shown in Figure 3, the VRH regime corresponds to 3D transport because the parameters of the anisotropic system still allow for significant hopping rates in the transverse direction. For 1D VRH to take place, the degree of anisotropy should be much stronger.
The orientational dependence of PF shown in Figure 2f implies that a more effective power generator can be built by aligning the localized states to some moderate degree but not to the full extent. This might naturally be the case in the experimental samples due to misalignment at the boundaries between crystalline grains. 23 It is interesting that an ideal single crystal for organic thermoelectrics (in the hopping regime) is not the solution because it would reduce PF. An experiment of gradual uniaxial alignment, including much higher degrees of alignment and concentration measurements, would be of interest to verify the results presented here.
The effects due to misalignment at the boundaries between crystalline grains might be included into the model by assigning different orientations O d ′ to specific parts of the computational domain. Figure 5 shows an example of the misaligned interlayer that divides the highly oriented sample into eight spatial portions. The interlayer here might be considered as an amorphous spacer separating highly oriented crystals of the same orientation O d . As O d ′ decreases, the ratio σ ∥ /σ ⊥ quickly drops, which is a result of the interlayer itself acting as a lowresistance shunt specifically in the perpendicular direction. This morphology model might be further improved by choosing more accurate values of ξ and more relevant geometries; however, such modeling is outside of the scope of this study. For heterogeneous morphology, the model by Kaiser 22,36 provides accurate description of oriented polymer films.
The theoretical results presented above agree qualitatively with the experimental data 1−7 for the ratios σ ∥ /σ ⊥ and S ∥ /S ⊥ . However, absolute values of σ and S differ; in particular, σ is a few orders of magnitude smaller. This can be explained by the usage of typical parameters [15][16][17][18]29,30 for modeling without any additional tweaking. Quantitative agreement for σ and S is left for the future study. Note that experimental data 1−7 largely vary from one sample to another because of the strong sensibility of morphology to the preparation process. For P3HT, the resulting polymer film might generally be amorphous, crystalline, or a mix of the two. 23 Thus quantitative estimation of the parameters entering the theory should be done for every sample individually.
Apart from the effects introduced by Coulomb interaction in comparison to the noninteracting theory as shown above, it is instructive to consider the consequences of truncating the longranged part of the electron−electron interaction. This can be done using the classic argument for the existence of the Coulomb gap in the distribution of energy levels of localized electrons in strongly disordered semiconductors. 12 Let us consider a pair of states i and j, respectively, above and below μ. The stability criterion for the ground state requires that 10 where Δ ij is positive and signifies that the energy of the ground state cannot be lowered by promoting an electron from j to i. The states on opposite sides of μ that differ in energy by less than a small value η must be separated in space by a distance larger than e 2 /η. Hence, the spatial DOS vanished at least as fast as (η/ e 2 ) d , where d is the dimensionality of the system. If long distances are removed from consideration, the smallness of η is never achieved and DOS does not vanish at μ. In structurally anisotropic materials, there might be other scattering mechanisms responsible for anisotropy in the thermoelectric properties. For example, when band transport dominates, 22 anisotropy in electron effective mass leads to anisotropy in electric field screening. 37 This present study is limited to inelastic transport, when momentum and phase are not conserved between consecutive hopping events and phonons act in the transport assisting way. 10 It would be interesting to make comparison with other theories, both on inelastic and elastic transport in anisotropic organic semiconductors, and it is hoped that the present exploration will stimulate development of such theories and comparison efforts.
Several final comments follow. First, if positional disorder is added to the modeling, the orientational dependence does not change appreciably. Second, both noninteracting and interacting models predict S−σ dependence to have a fall-off shape, similar to other hopping theories, 1,4 but in contrast to experimentally observed trend S ∝ σ −1/4 . 38 That was shown to be due to a limitation of the VRH model itself. 39 Third, the hopping rates 2 assume electrons or holes as charge carriers. These rates are modified when polaron effects become important. Those effects, however, are expected to be small for the system parameters and the linear Ohmic regime studied here. 40

■ CONCLUSIONS
A charge hopping model is presented that accounts for the correlation between deformational and orientational degrees of freedom of the localized states, energetic disorder, electron− electron interactions, and charge concentration. For a parameter set similar to typical organic semiconductors, [15][16][17][18]29,30 it is shown that an increase of the degree of orientation causes an exponential increase of the ratio of conductivities in parallel and perpendicular directions, while the ratio of Seebeck coefficients remains nearly unaffected. However, the ratio of Seebeck coefficients can increase if Coulomb interaction is taken into consideration and charge screening in the direction parallel to the predominant orientation of the localized states is stronger than in the perpendicular direction. The regime of charge transport occurs at the crossover between VRH and NNH in both directions. These findings provide a microscopic explanation for the thermoelectric properties of anisotropic polymer films in recent experiments 1

■ ACKNOWLEDGMENTS
This work was supported by SNIC 2020/13-95. It is a pleasure to acknowledge discussion with X. Crispin.