We present a computational approach to model hole transport in an amorphous semiconducting fluorene-triphenylamine copolymer (TFB), which is based on the combination of molecular dynamics to predict the morphology of the oligomeric system and Kinetic Monte Carlo (KMC), parameterized with quantum chemistry calculations, to simulate hole transport. Carrying out a systematic comparison with available experimental results, we discuss the role that different transport parameters play in the KMC simulation and in particular the dynamic nature of positional and energetic disorder on the temperature and electric field dependence of charge mobility. It emerges that a semi-quantitative agreement with experiments is found only when the dynamic nature of the disorder is taken into account. This study establishes a clear link between microscopic quantities and macroscopic hole mobility for TFB and provides substantial evidence of the importance of incorporating fluctuations, at the molecular level, to obtain results that are in good agreement with temperature and electric field-dependent experimental mobilities. Our work makes a step forward towards the application of nanoscale theoretical schemes as a tool for predictive material screening.

## I. INTRODUCTION

Amorphous semiconducting polymers, due to their proven versatility, offer substantial advantages in terms of ease of processing and low production costs, with promising applications in the field of organic electronics, such as solar cells, light emitting diodes, and field effect transistors. Recent investigations suggest that charge transport in conjugated polymers is not simply governed by the degree of crystallinity, but in general by the presence of aggregates with short-range order, which promote inter-chain charge transport.^{1} Experimental evidence that aggregation of polymer chains with sufficient short-range order and interconnected domains, in a globally amorphous structure, is sufficient to promote charge transport^{2,3} points towards the refinement of models conventionally employed for describing charge transport in crystalline organic semiconducting systems.

As polymer chains are bound by weak interactions with multiple degrees of freedom, parameters influencing hole mobility that should be included in the modeling of hole transport include (i) the length of polymer chains, (ii) the conformational freedom along the chain subunits, (iii) the conformation space of the polymer, (iv) the number of possible neighboring units, and (v) the orientation of each chain with respect to the direction of the external electric field.^{4–8} The theoretical frameworks developed to address the role of morphology on electronic transport in semiconducting polymers demonstrate that intra-chain transport dominates at shorter length scales with high mobilities, up to 10-100 cm^{2} V^{−1}s^{−1},^{9} while inter-chain transport prevails at larger displacements, with mobilities, a few orders of magnitude lower. The slowest inter-chain transport defines the limiting charge carrier mobility attainable at mesoscopic device scale.^{5,10} Such theoretical approaches, in conjunction with Molecular Dynamics (MD) simulations, where reasonably accurate polymer morphologies can be obtained,^{11,12} provide substantial information about the effect of polymer morphology and local structure on charge transport properties of the polymeric system under study.^{13}

The main complexity originates from the need of incorporating the whole morphological variability of the macromolecular system into charge transport modeling. Indeed, a theoretical study at the atomistic level, which considers both intra- and inter-chain transport and explicitly takes into account all possible fluctuations over many time scales and morphologies, is a challenging task, both scientifically and in terms of computational resources.^{14} Such frameworks are particularly unpractical for preliminary material screening, whereas a simplified methodology, able to provide reliable results with respect to experimental observations at a low computational cost, would be more appropriate.

In what follows, we describe a step forward in the realization of this objective, in which, with a combination of MD and KMC simulations, we systematically compare computational predictions with experimental temperature-dependent mobilities. As a case study, we consider the [(9,9-dioctylfluorenyl-2,7-diyl) (4,$4\u2032$-(*N*-(4-*sec*-butylphenyl)))] diphenylamine alternating copolymer, also known as TFB (Fig. 1). TFB is an amorphous polymer^{15} with relatively high hole mobility^{16,17} which has been the subject of several fundamental studies,^{18–22} and it is employed in practical applications such as light emitting diodes^{23} and solar cells.^{24–26} Our approach is based on two main and interdependent assumptions (i) that the holes are localized along short segments of the polymer chain and (ii) that hole transport can be described as a hopping process between those oligomeric units. Evidence supporting these hypotheses is provided, and their consequences on calculated charge carrier mobilities are discussed in detail.

## II. COMPUTATIONAL DETAILS

### A. Molecular dynamics simulations

We focused on TFB oligomers composed of five 9,9-dioctylfluorene (FLU) units, alternated with four butyl triphenylamine (TPA) co-monomer units (see Fig. 1). This specific oligomer length was chosen as it corresponds to the upper limit of the conjugation length of TFB, as estimated on the basis of quantum chemical calculations reported by Sancho-García *et al.*^{27} Our choice of this oligomer approach^{11,28–30} in modeling charge transport in TFB was also motivated by practical reasons, such as the possibility of dealing with reasonable box sizes and short equilibration times.

We employed a united-atom^{31} force field complemented with quantum chemical calculations for atomic charges and torsional potentials, as detailed in the supplementary material. All molecular dynamics simulations were performed with the NAMD software.^{32} In order to simulate amorphous TFB samples, the starting geometries for the simulation were built by introducing a certain degree of randomness, according to the following procedure: (i) the starting geometry was generated by placing the oligomers with random orientations on the nodes of a cubic mesh of 5 $\xd7$ 5 $\xd7$ 5 sites, with a lattice constant of 80 Å; (ii) the low-density sample was compressed at 1000 K and 100 atm up to a rough volume stabilization, which occurs in about 1 ns; (iii) high-temperature annealing was performed for 10 ns of simulation at 1000 K and 1 atm; (iv) equilibration was carried out for 20 ns of simulation at *T*_{MD} = 300 K and 1 atm; (v) this procedure was applied to simulate four different and uncorrelated samples composed of 125 oligomers (33 380 united atom centers). A snapshot of one of the four samples after equilibration at 300 K is shown in the right panel of Fig. 1. The final density of the four samples is $0.987\xb10.006$ g/cm^{3}. It is worth noting that since the equilibration temperature is well below the experimental glass transition, i.e., 393-415 K depending on film thickness,^{15} significant changes on the morphology are not expected at lower temperatures. For this reason and for keeping the approach computationally efficient, all charge transport simulations were conducted on the 300 K morphologies, neglecting the density variation with temperature. It is worth noting that this preparation scheme does not take into account any solvent inclusion in the polymer and therefore leads to ideal, completely amorphous morphologies in which the number of contacts between the chains is maximized.

### B. Charge transport model

Charge transport in amorphous TFB was described with a hopping model. Specifically, we resort to the Marcus formula for hopping rates,^{33} describing non-adiabatic hole transfer in the weak coupling regime, which has been widely employed to compute the charge transfer rates between adjacent molecules in organic semiconductors.^{34} The charge transfer rate in the semi-classical Marcus formalism reads as

where *e* is the elementary charge, $\mathit{\epsilon}$_{i} and $\mathit{\epsilon}$_{j} are the energies of sites *i* and *j* involved in the charge transfer, ΔG_{ij} is the energy difference between the initial and final states, *J*_{ij} is the electronic coupling (or transfer integral), *λ* = *λ*_{i} + *λ*_{e} is the reorganization energy, *k*_{B} is the Boltzmann constant, and *T* is the temperature. In the presence of an applied electric field ($E\u2192$), ΔG_{ij} also includes the contribution $eE\u2192\u22c5R\u2192ij$, with $R\u2192ij$ being the site-site distance vector. As positions of the sites, we used the center of the Mulliken charges evaluated from single point AM1 calculations^{35} for positively charged oligomers at their instantaneous MD geometries. The reorganization energy consists of (i) an internal contribution (*λ*_{i}) associated with the change in geometry upon charge transfer and (ii) an external contribution (*λ*_{e}) that reflects the structural and electronic changes in the surrounding medium. The intramolecular part was set *λ*_{i} = 0.1 eV on the basis of DFT calculations carried out on oligomers of increasing length (see supplementary material), a value very similar to the one reported for the well-known poly(3-hexylthiophene) (P3HT)^{36} and for the TPA moiety.^{37} The external contribution, a challenging quantity to calculate exactly,^{38} was set to *λ*_{e} = 0.2 eV in the majority of the calculations, a value expected for a low dielectric constant medium according to classical continuum model calculations,^{39} and in line with theoretical estimates for oligoacene crystals.^{38}

For all molecules of the four MD samples, the following conditions were applied: (i) for quantum chemistry calculations, the hydrogen atoms, absent in the MD model, were added to simulated molecular structures using geometrical criteria; (ii) hole transfer integrals *J*_{ij} were calculated as one-electron HOMO-HOMO couplings using the dimer projection method^{40} at the ZINDO level of theory with the Mataga-Nishimoto potential;^{41,42} (iii) periodic boundary conditions were employed to select the pairs of molecules in contact, using 6 Å as the interatomic cutoff distance between atoms belonging to different molecules; (iv) site energies $\mathit{\epsilon}$_{i} and $\mathit{\epsilon}$_{j} were approximated by ZINDO HOMO energies, including the effect of distortion of the molecular geometry in realistic conformations.

Although the energetic disorder introduced by intermolecular electrostatic interactions, due to permanent and induced dipoles, can play a key role in improving charge separation at heterojunctions^{43–45} and charge-carrier localization,^{46} it was not included in this work because of the nonpolar nature of TFB, which implies a low energetic disorder as confirmed by Malliaras *et al.*^{18}

### C. Charge propagation

The propagation of charges is simulated with the Kinetic Monte Carlo (KMC) technique. The charge is initially located at a randomly chosen site *i*. The hopping rates *k*_{ij} from donor site *i* to any potential acceptor site *j*, computed using Eq. (1), are employed to generate instantaneous hopping times *τ*_{ij} in the framework of the “first reaction” method,^{47}

where *ξ* is the random number drawn from a uniform distribution within the interval 0 to 1. The destination site $i\u2032$ is selected from the set of available neighbor sites *j* as the one having the smallest reaction time $\tau ijmin$. The simulation progresses by following the subsequent hopping steps during charge propagation and the simulated elapsed time is increased by $\tau ii\u2032$ at each hopping step. All the KMC simulations were run assuming a linear voltage profile across the sample, by varying the magnitude of the electric field from $E=1\u22c5104$ V/cm to $30\u22c5104$ V/cm. In each simulation, a single charge is propagated, using 3D periodic boundary conditions, until it covers a fixed distance *d* = 4 *μ*m along the direction of the electric field vector. Since only one charge is present at a time, the simulation targets the limit of very low charge densities following experimental evidence of charge density-independent mobility for a polymer very similar to TFB.^{48} For each charge propagation run *k*, mobility is then estimated as follows:

where *τ*_{k} is the sum of all hopping times $\tau ii\u2032$, i.e., the time needed for a charge propagating in the field direction to cover the distance *d*. Unless explicitly stated, mobility values reported in this work are averages over 500 up to 750 KMC runs performed with the electric field vector parallel to *x*, *y*, and *z* Cartesian directions and over the four MD samples. Average values are reported either as arithmetic means, $\mu =(\u2211k\mu k)/N$, or as logarithmic means $\mu ln=exp(\u2211kln\mu k/N)$, as the latter is a better estimator of the mobility in small samples.^{49} We noticed that for our sample sizes and propagation distances, the two means are almost coincident, indicating a satisfactory convergence of the KMC simulations.

### D. Analysis of charge mobilities

According to Bässler,^{50} at intermediate fields and in conditions of thermally assisted hopping, the dependence of mobility on the electric field modulus can be described by the following empirical relation:

where $\beta \u2009>\u20090$ is the Poole-Frenkel factor. Using zero field, temperature-dependent mobilities *μ*_{0}(*T*) derived from Eq. (5), the “apparent” diagonal energetic disorder *σ*_{A} and the infinite temperature mobility ($\mu \u221e$) may in turn be obtained by extrapolating the temperature dependence of *μ*_{0} with

The Poole-Frenkel factor depends on both apparent diagonal (*σ*_{A}) and off-diagonal (Σ_{A}) disorders,

with *C* being an empirical proportionality factor. The combination of Eqs. (5)–(7) provides an approximate^{51} but universal equation for the dependence of mobility on the temperature, field, and energetic disorder. We use this theoretical framework to compare our data with experimental results obtained by Malliaras and co-workers.^{18,52}

In the original Gaussian disorder lattice model (GDM),^{50} *σ*_{A} and Σ_{A} correspond to the standard deviations of the distributions of site energies (the density of states, DOS) and transfer integrals. In the present work, we label these two parameters as “apparent” since they are effective quantities extracted by fitting the mobility obtained from our off-lattice KMC simulations with the equations above, which are strictly valid only for the Gaussian disorder model. The apparent tag serves to distinguish the phenomenological disorder parameters obtained by such a fit from the microscopic parameters employed in our KMC simulations. Although the GDM model prescribes a 1/*T*^{2} dependence of the logarithm of the mobility, it has been shown experimentally that often the temperature dependence of the mobility of organic semiconductors at low fields,^{53} and for exponential densities of states,^{54} can be well reproduced with an Arrhenius law, a behavior also well established for inorganic disordered semiconductors,^{51}

We used the latter equations to fit the temperature dependence of the zero field mobility, obtained by extrapolating simulation and experimental data with Eq. (6). The infinite temperature mobility in Eq. (8) is primed because it is different (higher) with respect to the one obtained by fitting with Eq. (6).

## III. RESULTS AND DISCUSSIONS

### A. Charge transport parameters

We initially analyzed the fluctuations of site energies in the samples, which generate the so-called diagonal (energetic) disorder, known to strongly influence charge transport.^{46,50,55–58} The distribution of site energies (HOMO levels) calculated at MD geometries has an approximately Gaussian shape with standard deviation $\sigma \mathit{\epsilon}$ = 50.2 meV (Fig. S5 of the supplementary material), consistent with the assumptions of the Gaussian disorder model. This value is also close to the effective disorder *σ*_{A} = 65.9 meV reported by Fong *et al.*^{18} and was hence used in our KMC simulations.

Other fundamental parameters in determining the transfer rate, besides the reorganization energy discussed previously, are the magnitude of the electronic couplings between the different sites (the nodes of the charge percolation network) and the average number of pathways linking one node to the others, i.e., the number of neighboring sites with non-zero electronic couplings. When examining the calculated transfer integrals, it is interesting to notice (Fig. S6 of the supplementary material) the very low values of the couplings, with an average value of 0.40 meV, when, for instance, values of the order of $\u223c$10 meV are found for crystalline and amorphous fullerenes at short intermolecular distances.^{46} In addition, as it happens also for crystalline systems along the directions with small couplings,^{59} the standard deviation *σ*_{J} = 1.43 meV is higher than the mean value, suggesting an important contribution of dynamic fluctuations on the mobility values.^{60} Despite the low coupling values, the number of neighbors with non zero coupling ($J>0.1$ meV) ranges between 10 and 20 per oligomer, ensuring a high number of percolation pathways for charge transport. Such couplings most often correspond to triphenylamine-triphenylamine contacts, while fluorene-fluorene contacts are less frequent because of the presence of the bulky dioctyl substituent on one side of the indeno group (cf. Fig. 1). The chemical design principle of attaching the solubilizing side chains on the fluorene unit is of course beneficial for charge transport since the hole is mostly localized on a triphenylamine unit (87% according to AM1 calculations), and then the HOMO-HOMO contacts are maximized if this region of the backbone is not hidden by the alkyl chains.

### B. Mobility calculated for static networks

In Secs. III C–III E, we address in detail the relation between charge transport simulation parameters and calculated mobility as a function of temperature and electric field, utilizing the comparison with available experimental data as an assessment of the simulation predictions. In doing this, we implicitly assume that an oligomeric representation of the polymer is sufficient to describe the charge transport for longer chains; this means that the conjugation length of the charge is assumed to be shorter than the oligomer length and that inter- and intra-chain charge transport mechanisms are equally well described by oligomer-oligomer hopping and have identical characteristic times.^{5} A low conjugation length is expected for an amorphous and relatively flexible polymer like TFB, which was verified at the AM1 level by the observation that for our structures the extra charge is indeed localized on average within 1.5 units only. A similarly short localization length was reported by Sancho-García *et al.* for polyfluorene and TFB.^{27} The molecular origin of the short conjugation length are the equilibrium values of the torsional angles along the chain, which are far away from the planar conformation, see also Figs. S3 and S5 of the supplementary material. Working at the low molecular weight limit implies that trapping at the chain ends is maximized,^{7} hence we would expect to somewhat underestimate experimental mobilities for long polymer chains, which for TFB are on the order of 10^{−2} cm^{2}/Vs at room temperature.^{18}

We started with running KMC simulations on the original MD samples (denoted as *R*^{1} in Table I) with site energies and transfer integrals calculated for each sample at the ZINDO level ($\mathit{\epsilon}$_{Q}, *J*_{Q}). In these simulations, charges are therefore propagating in a transport network whose properties are static in time, corresponding to single snapshots of the samples. As discussed in the literature,^{36,61,62} transfer integrals fluctuate over time scales often shorter than the typical hopping times and taking a static picture eventually leads to an underestimation of materials mobility. The KMC output at different electric fields is reported in Fig. 2. Apart the underestimation of the experimental mobilities of one or two orders of magnitude, what is striking in Fig. 2 is the absence, at least at the higher temperatures, of the Poole-Frenkel behavior ($ln\mu \u221dE1/2$) reported experimentally. To demonstrate that this offset is not originated by an unrealistically high value of the reorganization energy, we report in the same panels also mobilities calculated with *λ* = 100 meV (empty symbols), which can be considered as a lower limit for this parameter.^{63} The calculated mobilities for *λ* = 100 meV are of course larger than the ones obtained at *λ* = 300 meV (filled symbols), consistent with Eq. (1), but they are still one order of magnitude lower than the experimental values and show a decreasing trend of mobility with the field at high temperatures.

Symbol . | Definition . |
---|---|

R^{1} | Original four MD samples (125 sites each) |

R^{2} | 2 $\xd7$ 2 $\xd7$ 2 replica of MD samples (1000 sites) |

R^{3} | 3 $\xd7$ 3 $\xd7$ 3 replica of MD samples (3375 sites) |

$\mathit{\epsilon}$_{Q} | HOMO energies calculated at the ZINDO level at the MD |

geometry, static in time with $\sigma \mathit{\epsilon}$ = 50.2 meV | |

$\mathit{\epsilon}$_{G} | HOMO energies drawn at every KMC step from a Gaussian |

distribution with $\sigma \mathit{\epsilon}$ = 50.2 meV | |

$\mathit{\epsilon}GS$ | HOMO energies drawn from a Gaussian distribution at the |

beginning of each KMC simulation, with $\sigma \mathit{\epsilon}$ = 50.2 meV | |

J_{Q} | Electronic couplings fixed at their quantum mechanics value |

J_{G} | Electronic couplings fixed at their quantum mechanics value |

+ random fluctuations x extracted at every KMC step from a | |

Gaussian distribution with σ_{J} = 1.43 meV: J = J_{Q} + x | |

$JQT$ | Electronic couplings augmented according to a temperature- |

dependent disorder: $J2=JQ2+\sigma J2\u22c5(T/TMD)$ |

Symbol . | Definition . |
---|---|

R^{1} | Original four MD samples (125 sites each) |

R^{2} | 2 $\xd7$ 2 $\xd7$ 2 replica of MD samples (1000 sites) |

R^{3} | 3 $\xd7$ 3 $\xd7$ 3 replica of MD samples (3375 sites) |

$\mathit{\epsilon}$_{Q} | HOMO energies calculated at the ZINDO level at the MD |

geometry, static in time with $\sigma \mathit{\epsilon}$ = 50.2 meV | |

$\mathit{\epsilon}$_{G} | HOMO energies drawn at every KMC step from a Gaussian |

distribution with $\sigma \mathit{\epsilon}$ = 50.2 meV | |

$\mathit{\epsilon}GS$ | HOMO energies drawn from a Gaussian distribution at the |

beginning of each KMC simulation, with $\sigma \mathit{\epsilon}$ = 50.2 meV | |

J_{Q} | Electronic couplings fixed at their quantum mechanics value |

J_{G} | Electronic couplings fixed at their quantum mechanics value |

+ random fluctuations x extracted at every KMC step from a | |

Gaussian distribution with σ_{J} = 1.43 meV: J = J_{Q} + x | |

$JQT$ | Electronic couplings augmented according to a temperature- |

dependent disorder: $J2=JQ2+\sigma J2\u22c5(T/TMD)$ |

### C. Effect of sample size and of dynamic diagonal disorder

In order to understand the factors that contribute to the large offset between simulations and experiment, KMC simulations were carried out by considering different conditions: (i) sample size, varied by replicating in space the original MD box (R^{1}) two or three times along each Cartesian direction (*R*^{2} and *R*^{3}, respectively); (ii) diagonal disorder, with site energies assigned from quantum chemistry calculations ($\mathit{\epsilon}$_{Q}), or randomized by extracting site energy values from a Gaussian distribution having the same standard deviation of the calculated HOMO energies, either at the beginning of each KMC simulation (static, $\mathit{\epsilon}GS$) or at each KMC step (dynamic, $\mathit{\epsilon}$_{G}); (iii) fixed (*J*_{Q}) or dynamic off-diagonal disorder (*J*_{G}). While the fixed case for transfer integrals has an exact correspondence with the fixed site energies, in the dynamic case, disorder is added on top of the quantum chemistry values, using random numbers drawn from a Gaussian distribution with *σ*_{J} = 1.43 meV, corresponding to the standard deviation of ZINDO values. Finally, for (iv), we ran KMC simulations in which off-diagonal disorder, which we estimated at *T*_{MD} = 300 K, is assumed to be completely dynamic and is rescaled to the simulation temperature *T* using a factor (*T*/*T*_{MD})^{1/2}, according to the electron-phonon coupling model.^{42,64} Furthermore, the effect of increasing the mean values of the transfer integral is considered in an average way, by adding the variance to the square of the quantum mechanical (QM) transfer integrals^{65,66} ($JGT$ conditions). The meaning of the symbols corresponding to different simulation conditions is detailed in Table I.

The limitations of calculating mobilities from small-sized systems, yielding unrealistically high values and excessively broad distributions of single charge velocities, have been well discussed in the recent literature.^{28,49,56,67} Among others, Andrienko and co-workers^{56,68} underlined the relationship with the calculated mobility dispersion, energetic disorder, and sample size in amorphous systems and showed that for large disorder values ($\sigma \mathit{\epsilon}>3kBT$), only samples containing millions of sites can lead to correct results. The problem originates from the fact that the thermally accessible, low energy tail of the density of states, which becomes populated if the charge propagation experiments are adequately long, is not well represented by a discrete distribution of site energies (Fig. S5 of the supplementary material), in particular, if the site energies are fixed in time and the distribution is broad (i.e., the disorder is high).

In order to alleviate sample size effects and test their importance as a possible cause of the mismatch with experimental mobilities, we (i) introduced a static energetic disorder by randomizing the site energies at the beginning of each KMC run, (ii) we performed the simulations for periodic replicas of the original MD boxes, so as to have a larger number of sites and then a smoother density of states, and (iii) we calculated mean logarithmic mobilities as suggested by Bobbert and co-workers.^{49} By comparing the simulated mobilities for replicated samples *R*^{2} and *R*^{3}, we first notice that the differences between the two systems are minimal. In fact for TFB, the sample size is not expected to have a significant impact on mobility due to the small energetic disorder present in these systems ($\sigma \mathit{\epsilon}<3kBT$) at every temperature studied. The convergence with the system size is confirmed by the similarity between arithmetic and logarithmic averages (points and lines in Fig. 3, respectively), and consequently, in the remainder of the article, only logarithmic mobilities are reported. More interestingly, by comparing Figs. 2 and 3, it turns out that the static randomization of site energies scales down the mobility by a factor of about five. This reduction can be explained by the removal of (spurious) static correlations between the site energies which is instead present if $\mathit{\epsilon}$_{Q} values are used; in fact, as shown, for instance, by Kordt and Andrienko,^{69} uncorrelated disorder generally leads to lower mobilities.

To investigate the existence of spatial correlation in *R*^{1}|$\mathit{\epsilon}$_{Q}|*J*_{Q} simulations, we take a closer look into the mobility anisotropy for one of the four original MD samples (Fig. 4, filled symbols). Quite strikingly, there are marked differences among the three directions of application of the electric field, which in principle is not expected for an amorphous material. This behavior can be attributed to the relatively small size of the simulation samples (with respect to the oligomer size) since a few hopping events are sufficient for a charge to travel across the whole box. The randomization of site energies (gray-shaded and empty symbols) seems to fix the issue and improve the validity of transport landscape by yielding more isotropic mobility values, albeit two orders of magnitudes lower than the experimental ones and decreasing with increasing electric field. The mobilities obtained with frequent randomization of site energies ($\mathit{\epsilon}$_{G}) are as expected slightly larger (and more isotropic) than the ones with static random site energies ($\mathit{\epsilon}GS$),

### D. Mobility for dynamic networks

We showed so far that despite adopting reasonable approximations and parameters, calculated mobilities for static networks substantially differ from the experimental time-of-flight values. One possible source of error could be the adoption of the semiempirical ZINDO Hamiltonian in the calculation of the electronic couplings. ZINDO is in fact known to underestimate the couplings with respect to DFT calculations (which in turn are functional-dependent^{70}); however, the scaling factor ranges approximately from one to two depending on the material,^{46,71} which is not sufficient to recover two orders of magnitude in predicted mobilities. The other possible source of discrepancy with respect to experimental observations is the neglect of time fluctuations of the couplings, i.e., the static description of the percolation network. Cornil and co-workers elegantly showed how including fluctuations in KMC simulations increases the calculated mobility when the average value of the couplings is smaller than their standard deviation ($\u27e8J\u27e9/\sigma J<1$).^{60} The existence of these fluctuations, essentially generated by low frequency vibrations and with typical correlation times in the subpicosecond range, were observed by many authors.^{36,61,62,72–74}

It then appears that a sensible approach to introduce dynamic disorder on transfer integrals, as an effective way of considering small fluctuations of atomic positions around their equilibrium value, the so-called positional disorder.^{50} The electronic couplings calculated from the four MD samples indeed provide an instantaneous photograph of the network which magnifies the hindering effects on mobility due to trap sites. In order to improve the realism of our model, we accordingly reshuffled the transfer integrals at each KMC step by adding to each ZINDO-calculated coupling a random number extracted from a Gaussian distribution with *σ*_{J} = 1.43 meV. It is worth noting that this scheme (*J*_{G}) (i) assumes that the correlation time of transfer integral fluctuations is faster than the average hopping time; (ii) preserves the topology of the network, but allows for a time-dependent probability of the different pathways; (iii) accounts for the distance dependence of the transfer integrals; (iv) presumes that spatial fluctuations are a measure of the fluctuations in time, as is expected for a crystal, and (v) neglects the temperature dependence of the disorder. A more rigorous way of accounting for dynamic off-diagonal disorder, at least in the framework of the Marcus theory, is to augment all the squared transfer integrals of the dynamic variance^{65} and to consider that this variance is proportional to temperature^{42} ($JQT$ conditions).

In Fig. 5, left panel, we plot the results of KMC simulations performed on replicated sample *R*^{3} and compare them with the experimental data on the right panel. Clearly, introducing dynamic fluctuations on the electronic coupling increases significantly the simulated mobility and is the key to improve the agreement with experimental data, which confirms the importance of including dynamic fluctuations in charge transport studies of amorphous systems characterized by small values of electronic couplings.

The last step in this systematic review of the KMC simulation parameters concerns the temperature dependence of the positional (off-diagonal) disorder. So far we implicitly made the assumption of a disorder independent on temperature, but it is worth testing the opposite case and how much this would affect our conclusions ($\mathit{\epsilon}GS|JQT$ conditions). We see in Fig. 5 that the temperature-dependent disorder on transfer integrals does not produce qualitative changes since it marginally increases the mobility and the separation between the curves at different temperatures, improving the agreement with the experimental data from Malliaras and co-workers.^{18} The key factor for recovering the experimental order of magnitude of mobility and the Poole-Frenkel behavior is indeed the introduction of the fluctuations in the electronic couplings.

### E. Analysis of field and temperature dependence of mobilities

Analytical equations providing relatively simple formulae to fit the variation of mobility with field and temperature,^{50} as well as charge density,^{51,75} are often used in the experimental literature, with the double purpose of validating theoretical models and extracting intrinsic material parameters under the tacit assumption of the validity of the models. Here we perform the same exercise using simulated mobilities, in order to compare directly with the parameters extracted from experiment and to verify the internal consistency of the various theories developed to describe charge transport since some of the microscopic parameters, and notably the magnitude of the energetic disorder, are known beforehand.

In Fig. 6, we analyze the temperature dependence of zero field mobilities utilizing either Eq. (6), which is appropriate for hopping in a Gaussian density of states (left panel), or Eq. (8), typical of an exponential DOS, or of a Gaussian DOS in the presence of extrinsic charges (right panel).^{53} Even if the fitting with the Gaussian disorder model is more accurate as expected, also the Arrhenius equation gives satisfactory results although not in the whole temperature range. The two fits produce very different infinite temperature mobilities ($\mu \u221e$ $\u226a$ $\mu \u221e\u2032$ in Table II) as a consequence of the different temperature dependence; however, the apparent energetic disorders *σ*_{A} and activation energies Δ are rather consistent, with $\Delta \u22432\sigma A$ for all KMC setups and experiments. Changing the reorganization energy from *λ* = 200 meV to 300 meV does not alter significantly the activation energy, consistent with what was reported in Ref. 76 for oriented polymers and fields perpendicular to the direction of chain alignment, i.e., inter-chain transport.

Quantity . | Expt.^{18}
. | Expt.^{16}
. | R^{1}|$\mathit{\epsilon}$_{Q}|J_{Q}
. | R^{3}|$\mathit{\epsilon}$_{G}|J_{Q}
. | $R3|\mathit{\epsilon}GS|JG$ . | R^{3}|$\mathit{\epsilon}$_{G}|J_{G}
. | $R3|\mathit{\epsilon}G|JG$^{a}
. | $R3|\mathit{\epsilon}GS|JQT$ . |
---|---|---|---|---|---|---|---|---|

$\mu 0\u22c5103$^{b} | 9.6 | 0.31 | 0.33 | 0.14 | 3.1 | 4.5 | 14.7 | 5.9 |

$\beta \u22c5103$^{b} | 0.61 | 3.8 | -0.18 | -0.62 | 0.79 | 0.75 | 0.77 | 0.59 |

$\mu \u221e\u2032\u22c5103$ | 2200 | 17 | 46 | 920 | 870 | 1130 | 1990 | |

Δ | 140 | … | 100 | 155 | 149 | 135 | 111 | 154 |

$\mu \u221e\u22c5103$ | 190 | … | 2.4 | 2.5 | 77 | 71 | 130 | 160 |

σ_{A} | 66 | … | 54 | 67 | 71 | 64 | 62 | 69 |

$C\u22c5104$ | 2.8 | … | … | … | 3.1 | 1.3 | 1.4 | 3.3 |

Σ_{A} | 2.2 | … | … | … | 2.1 | 1.9 | 1.9 | 2.1 |

Quantity . | Expt.^{18}
. | Expt.^{16}
. | R^{1}|$\mathit{\epsilon}$_{Q}|J_{Q}
. | R^{3}|$\mathit{\epsilon}$_{G}|J_{Q}
. | $R3|\mathit{\epsilon}GS|JG$ . | R^{3}|$\mathit{\epsilon}$_{G}|J_{G}
. | $R3|\mathit{\epsilon}G|JG$^{a}
. | $R3|\mathit{\epsilon}GS|JQT$ . |
---|---|---|---|---|---|---|---|---|

$\mu 0\u22c5103$^{b} | 9.6 | 0.31 | 0.33 | 0.14 | 3.1 | 4.5 | 14.7 | 5.9 |

$\beta \u22c5103$^{b} | 0.61 | 3.8 | -0.18 | -0.62 | 0.79 | 0.75 | 0.77 | 0.59 |

$\mu \u221e\u2032\u22c5103$ | 2200 | 17 | 46 | 920 | 870 | 1130 | 1990 | |

Δ | 140 | … | 100 | 155 | 149 | 135 | 111 | 154 |

$\mu \u221e\u22c5103$ | 190 | … | 2.4 | 2.5 | 77 | 71 | 130 | 160 |

σ_{A} | 66 | … | 54 | 67 | 71 | 64 | 62 | 69 |

$C\u22c5104$ | 2.8 | … | … | … | 3.1 | 1.3 | 1.4 | 3.3 |

Σ_{A} | 2.2 | … | … | … | 2.1 | 1.9 | 1.9 | 2.1 |

^{a}

Simulation results for *λ* = 0.2 eV, *λ* = 0.3 eV otherwise.

^{b}

Quantity determined at T = 295 K.

The variations of a few meV of *σ*_{A} upon changing the simulation parameters are too small for drawing conclusions about their origin, but it is interesting to notice that all simulations produce similar apparent energetic disorder *σ*_{A}, in line with the experimental value of 66 meV, and only slightly larger than the intrinsic energetic disorder $\sigma \mathit{\epsilon}$ = 50.2 meV. The most significant change with KMC conditions, as shown already in Figs. 2, 4, and 5, occurs to the infinite temperature mobilities $\mu \u221e$; in particular, the best match with the experiments of Malliaras *et al.*^{18} is achieved with static diagonal disorder and temperature dependent off-diagonal disorder ($R3|\mathit{\epsilon}GS|JQT$). In Fig. 6, we also plot the zero field mobilities as obtained with *R*^{3}|$\mathit{\epsilon}$_{G}|*J*_{G} conditions but with a lower reorganization energy (*λ* = 0.2 eV, green pentagons), with the dual objective of demonstrating the robustness of the results versus the partly arbitrary choice of this parameter, and of showing how the experimental data are well bracketed by the KMC results with *λ* = 0.2 and *λ* = 0.3 eV (red squares).

To conclude our discussion, we focused on the material parameters that are extracted from fitting the temperature dependence of the Poole-Frenkel factors *β* with Eq. (7) (Fig. 7). The analysis could be performed only for some simulation setups because not all of them showed a Poole-Frenkel behavior [Eq. (5)]; the approach consists in utilizing previously derived *σ*_{A} values in Eq. (7), obtaining the parameter *C* from the slope of $log\beta $ vs (*σ*_{A}/*k*_{B}*T*)^{2}, and finally obtaining the square of the positional disorder *Σ*_{A} as the intercept of the fitting line with the x-axis. Figure 7 shows that experiment and simulations follow closely the behavior predicted by Bässler;^{50} however, the angular coefficient obtained from simulations is about half the experimental one when site energies are dynamically changed during KMC (see also Table II). If the diagonal disorder is instead assumed as completely static, we see indeed in Fig. 7 and Table II that larger values of *C* are obtained, very close to the experimental one. The positional disorder parameter *Σ*_{A}, which is a measure of the spread of the transfer integrals, shows instead a better overall agreement between experiment and simulations, independently from the details of the KMC algorithm, indicating that our model adequately captures the physical nature of the actual system.

## IV. CONCLUSIONS

We presented a theoretical study of the charge transport properties of a TFB fluorene-triphenylamine copolymer, a system with remarkably high hole mobility among amorphous organic semiconductors. KMC simulations based on atomistic input, i.e., atomistic molecular dynamics to build the morphologies and quantum chemistry to extract charge hopping rates, are employed to calculate the temperature and electric-field dependence of the hole mobility. Our analysis allows rationalizing why amorphous semi-conducting systems, characterized by poor intermolecular packing, relatively low energetic disorder, and small charge transfer integrals, can have reasonably high hole mobility, highlighting the crucial role of dynamic fluctuations on charge transport. The introduction in the model of dynamic off-diagonal disorder in fact boosts mobility by two orders of magnitude, allowing the achievement of an almost quantitative agreement with experimental mobility values, and the recovery of the so-called Poole-Frenkel behavior of the mobility against the electric field. These results provide a strong evidence of the importance of dynamic energetic disorder for the effective simulation of charge transport in amorphous semiconductors and introduce a relatively simple and cheap protocol for the computational screening of this class of materials.

## SUPPLEMENTARY MATERIAL

See supplementary material for a detailed description of the molecular mechanics force field parametrization; radial and dihedral distributions from MD simulations; distributions of transfer integrals and site energies; additional mobility vs electric field plots.

## ACKNOWLEDGMENTS

The work at University of Bordeaux was financed by the French national Grant No. ANR-10-LABX-0042-AMADEus managed by the National Research Agency under the initiative of excellence IdEx Bordeaux programme (Reference No. ANR-10-IDEX-0003-02). Computer time for this study was provided by “Mésocentre de Calcul Intensif Aquitain” computing facilities at Université de Bordeaux and Université de Pau et des Pays de l’Adour. The work at Georgia Tech was supported by the Office of Naval Research (Award No. N00014-17-1-2208). Y.Y. acknowledges financial support from the National Natural Science Foundation of China (21373229). This project has also received funding from the European Union Horizon 2020 research and innovation programme under Grant Agreement No. 646176 (EXTMOS).