Crossover from Hopping to Band-Like Charge Transport in an Organic Semiconductor Model: Atomistic Nonadiabatic Molecular Dynamics Simulation.

The mechanism of charge transport (CT) in a 1D atomistic model of an organic semiconductor is investigated using surface hopping nonadiabatic molecular dynamics. The simulations benefit from a newly implemented state tracking algorithm that accounts for trivial surface crossings and from a projection algorithm that removes decoherence correction-induced artificial long-range charge transfer. The CT mechanism changes from slow hopping of a fully localized charge to fast diffusion of a polaron delocalized over several molecules as electronic coupling between the molecules exceeds the critical threshold V ≥ λ/2 (λ is the reorganization energy). With increasing temperature, the polaron becomes more localized and the mobility exhibits a "band-like" power law decay due to increased site energy and electronic coupling fluctuations (local and nonlocal electron-phonon coupling). Thus, reducing both types of electron-phonon coupling while retaining high mean electronic couplings should be part of the strategy toward discovery of new organics with high room-temperature mobilities.

Forming the active layers of organic field effect transistor, light-emitting diodes, and organic photovoltaics, organic semiconductors (OSs) have enabled disruptive technologies in the plastics electronics industry and in the renewable energy sector. OSs combine many advantages: they are flexible and lightweight and some of their properties can be easily tuned using synthetic organic chemistry. Yet, a disadvantage is their limited charge carrier mobility compared to, e.g., inorganic semiconductors such as Si. This limits the range of applications of OSs and/or their efficiency. For instance, it is thought that higher carrier mobilities could decrease charge recombination in organic photovoltaics and improve efficiency. Therefore, the search for high mobility but cheap OSs remains an active field of research involving experimental, theoretical, and computational approaches.

Progress toward this objective is hampered by our incomplete understanding of charge transport (CT) in OSs including the simplest single-crystalline materials. CT in OSs is often in a difficult regime where standard transport theories, such as band theory or charge carrier hopping, do not apply.[1−5] The reason is that none of the parameters that determine the coupled electron–nuclear dynamics in these materials are “small” and can be neglected:[1,6−9] electronic coupling between the molecules (V = 10–200 meV), reorganization free energy or local electron–phonon coupling (λ = 100–200 meV), vibrational energies (ℏω = 10–200 meV), and thermal energy (kBT = 10–50 meV) often differ by not more than an order of magnitude, placing CT “between” the ranges of validity of band and hopping theory. For an in-depth discussion of this issue, we refer the interested reader to review articles.[1,3−5] There are successful approaches that reconcile the effects of electronic and nuclear motions in extended polaronic band theories through inclusion of both local and nonlocal electron–phonon couplings at the level of model Hamiltonians (Holstein–Peierl Hamiltonian).[2,10,11] However, at high temperature, the polaronic band description, which assumes delocalized states, no longer applies as the mean-free path of the charge carrier becomes comparable to the lattice spacing.[1]

Arguably, direct propagation of the electron–nuclear dynamics is one of the most appealing approaches to CT in OSs as it is, in principle, free of model assumptions and permits an “ab initio” view into the CT mechanism. In practice, approximations have to be made, but they can often be improved, sometimes systematically. A common approximation is the classical treatment of nuclei leading to so-called mixed quantum–classical nonadiabatic molecular dynamics methods (MQC-NAMD), where the charge carrier is propagated quantum mechanically in the time-dependent potential of the classical nuclei.[12,13] Previous MQC-NAMD simulations of CT in OS models differed in the level of theory used for the electronic Hamiltonian (Holstein, Peierl, or Holstein–Peierl type,[14−16] approximate DFT,[17] SCC-DFTB[18,19]) and the way nuclear dynamics, especially the feedback from electronic to nuclear degrees, was treated (no feedback,[14,17] Ehrenfest,[20] fewest switches surface hopping (FSSH)[15,16]). Some of the most notable findings in this respect were reported by Troisi et al.[1,14] and later by Wang and co-workers.[15,16] Propagating the carrier wave function along a chain of sites described by a model Hamiltonian, both predicted a decrease in mobility with increasing temperature in the large electronic coupling regime (defined below), similar to polaronic band-like transport but without assuming delocalized carriers. While Troisi’s work emphasized the crucial role of thermal fluctuations of electronic coupling (nonlocal electron–phonon coupling, Peierl term),[14,21] Wang concluded that thermal site energy fluctuations (local electron–phonon couplings, Holstein term) are the major source of polaron localization and mobility decay in this regime.[15]

In our previous works, we developed an efficient FSSH method termed fragment orbital-based surface hopping (FOB-SH)[22−24] that includes atomistic and electronic structure details in the dynamics, thereby going beyond previous model Hamiltonian studies.[14−16] A surface hopping method, it also includes feedback from the electronic to the nuclear degrees of freedom, which was often missing in previous propagation schemes.[14,17] In our quest for a better understanding of CT in OSs, we apply here an improved version of our FOB-SH method (details below) to investigate the temperature dependence of electron–hole mobility in a chain of ELMs in parameter regimes that are typical for application-relevant OSs.

We find that in the regime of small electronic coupling strengths V ≪ λ/2 (λ is the reorganization energy), the hole wave function is strongly localized on a single molecule at all temperatures and moves slowly along the chain via infrequent hops. In this regime, the temperature dependence of hole mobility goes through a maximum at around room temperature, which can be explained by a simple hopping model with rates obtained from electron transfer (ET) theory. Importantly, as electronic coupling is increased beyond the critical threshold V ≥ λ/2,[1,5,25] the wave function forms a polaron that is delocalized over several sites (as reported experimentally for some OSs[26−28]) and that moves rapidly along the chain. With increasing temperature, the polaron becomes more localized concomitant with a decrease in mobility according to a power law, μ ∝ T–1.2. We find that polaron localization and mobility decay in this regime are primarily due to site energy fluctuations. Electronic coupling fluctuations lead to a further increase in polaron localization and decrease in mobility, but this additional effect is not very large. Our atomistic FOB-SH simulations confirm the qualitative picture obtained from previous model Hamiltonian studies of Troisi et al.[14] and Wang and co-workers.[15,16] They reinforce the view that the “band-like” mobility decay with temperature, observed experimentally,[27] is due to increasing dynamic localization of the polaronic charge carrier with temperature.

In the following, we briefly introduce the previously developed FOB-SH methodology[22−24] and describe two new algorithmic developments for detection of trivial crossings and for removal of unphysical long-range charge transfer that is introduced by common correction schemes for electronic decoherence. In FOB-SH, it is assumed that the complicated many-body dynamics of an excess electron or electron–hole can be effectively described by a time-dependent one-particle wave function, Ψ(t). The latter is expanded in an orthogonalized basis of frontier molecular orbitals that mediate the CT, {ϕ}. For hole transport, the orbitals ϕ are the highest occupied molecular orbitals (HOMOs) of the molecules forming the material, and for electron transport, they are the lowest unoccupied molecular orbitals (LUMOs),where R(t) denotes the time-dependent nuclear positions. We refer to the set of orbitals {ϕ} as diabatic basis even though the nonadiabatic coupling elements (NACEs) d are in general nonzero, albeit small along the entire dynamics (typically ∼0.3 meV/ℏ). Insertion of eq into the time-dependent Schrödinger equation gives the familiar equation of motion for the expansion coefficients u in the diabatic basiswhere are the elements of the diabatic electronic Hamiltonian matrix and are the elements of the matrix of NACEs, . A key feature of FOB-SH is that explicit electronic structure calculations of the elements H are avoided during time propagation, which allows us to investigate large systems and long time scales. The diagonal matrix elements of the electronic Hamiltonian, H = ⟨ϕ|H|ϕ⟩, are the total electronic energy of the system when the charge is localized on molecule k (site energy). It is approximated here with a classical force field (see the Supporting Information). The off-diagonal matrix elements, H = ⟨ϕ|H|ϕ⟩, and the NACEs, d, are calculated using our recently developed analytic overlap method (AOM).[29] AOM couplings have been benchmarked before against wave function theory-validated fragment orbitals-based DFT (FODFT)[29] and against approximate coupled cluster (SCS-CC2)/Generalized Mulliken Hush calculations[30] on a diverse set of π-conjugated organic donor–acceptor compounds. It was found that errors are less than a factor of 2 with respect to reference coupling values that spanned 5 orders of magnitude. We refer to recent papers for a detailed discussion of this method in the context of FOB-SH.[22,24] The thermal average of the off-diagonal elements H, V = ⟨H2⟩1/2, is simply referred to as electronic coupling. We have dropped the indices from V because it is, to a good approximation, the same for each neighboring pair k,l in a chain of identical molecules.

As shown before,[31] in FSSH, the nuclear degrees of freedom should be propagated on one of the adiabatic potential energy surfaces E, , , with the unitary matrix transforming from the diabatic {ϕ} to the adiabatic electronic state basis {ψ}, . For an explicit expression of the nuclear forces on the adiabatic states and for the nonadiabatic coupling vectors (NACVs) between them, we refer to our recent work.[24] The probability to hop from the current (active) adiabatic state i to another state j needs to be evaluated at every MD time step Δt and is given according to Tully[32] bywhere a = cc* is the electronic density matrix, c the expansion coefficients of the charge carrier wave function eq in the adiabatic basis, , and dad the NACEs between the adiabatic states, dad = ⟨ψ|ψ̇⟩. Whenever a surface hop occurs, the velocities are adjusted in the direction of the NACVs between the adiabatic electronic states.[33] This adjustment massively improves energy conservation and detailed balance compared to isotropic rescaling of velocities, as shown previously by Coker and Xiao[34] and by our group in the context of FOB-SH.[24]

Straightforward application of the FOB-SH algorithm outlined, i.e., solving eq for the electron dynamics and propagating the nuclei on one of the adiabats E as determined stochastically via eq , is hampered by a number of problems, especially when simulating CT:

(1) Trivial or unavoided crossings occur when two (or more) adiabatic potential energy surfaces with vanishing derivative couplings cross between subsequent time steps, e.g., due to thermal fluctuations. This causes a change in the state ordering that, if undetected, gives rise to serious problems: continuation of the nuclear dynamics on an incorrect active state, discontinuity in the nuclear forces deteriorating energy conservation, erroneous calculation of time derivatives, especially dad, deteriorating excited-state population and detailed balance, and, most seriously in the context of charge transport, spurious and unphysical long-range charge transfer events between spatially distant states of similar energy. The problem is magnified in the case of large systems with tens, hundreds, or more electronic states as, e.g., in organic semiconductors, where adiabatic states may be localized and spatially distant (hence have a small derivative coupling) but very close in energy.

The simplest solution to the trivial crossing problem is to use a sufficiently small MD time step, but this is impractical in most applications. Here we track the identity of the states along the dynamics by calculating the overlap O between adiabatic electronic states l at time t and i at time t – Δt, ∀l,i = 1,MTo find the state at t – Δt that corresponds to state l at time t, we determine the maximum overlap, Omax = max|O|, assumed to occur for state i = k. If 1 – Omax ≤ ϵ (ϵ is a constant set to 0.1), state k at time t – Δt is mapped on state l at time t. This mapping procedure is carried out for each state l at time t. Finally, the subset of states l at time t for which 1 – Omax > ϵ and the subset of unmapped states i at time t – Δt are ordered by state index and consecutively mapped onto one another. Similar state tracking strategies have been proposed before by Tretiak and co-workers[35] and Thiel and co-workers.[33] In addition to the state tracking algorithm, we adopt the self-consistent FSSH (SC-FSSH) method suggested by Prezhdo and co-workers[36] to improve estimation of the hopping probability eq between the active state and the state next higher in energy, which is prone to numerical inaccuracies at trivial crossings. We note that the SC-FSSH algorithm alone, i.e., without application of the state tracking algorithm, did not sufficiently cure the trivial crossing problem and its manifestations in our simulations.

(2) Electronic Decoherence. A well-known problem of FSSH is the overcoherence of the electronic wave function.[37−39] Adopting a decoherence correction is paramount to reach good internal consistency between the surface and wave function population and to avoid unphysical delocalization of the wave function when nuclei approach avoided crossings. In the present work, we use a decoherence algorithm based on exponential damping of all except the active (subscript “a”) adiabatic electronic states (i ≠ a)[39]The coefficient for state a, c, is scaled appropriately to ensure norm conservation.[40] Different energy-based[39,40] or force-based[41] decoherence times containing various adjustable parameters have been proposed based on different physical arguments. In the present work, we adopted the simplest possible, parameter-free, Heisenberg principle-related decoherence time, τ = ℏ/|E – E|, and note that other common choices give very similar results for charge mobility (see the SI).

(3) Decoherence Correction-Induced Spurious Long-Range Charge Transfer. In small systems with only a few electronic states, surface hops between localized but spatially distant electronic states are unlikely because of the small derivative coupling in this case. In large systems with a high density of electronic states, the probability for a single transition is still small, but the large number of such transitions (due to the high density of localized states) makes them more likely to occur. In this case, the active adiabatic electronic state is reassigned within one MD time step from an adiabatic wave function localized in one region of space, say ψ′(t – Δt), to another adiabatic wave function localized in a different region of space, ψ(t). While such transitions are not an artifact of the FSSH algorithm per se, the problem is that the decoherence correction eq tends to quickly collapse the charge carrier wave function eq from Ψ′ ≈ ψ′ to Ψ ≈ ψ. This results in unphysical long-range charge transfer and in divergent charge mobilities with respect to system size. The same problem has also been reported by Wang and co-workers in a very recent study.[16] The authors refer to it as the decoherence-enhanced trivial crossing problem. We prefer to call it decoherence correction-induced spurious long-range charge transfer to emphasize that its origin lies in the approximate treatment of decoherence.

Here we implement a simple correction scheme to this problem. We define a moving and flexible active region containing virtually all of the charge carrier density |Ψ|2, eq . In practice, we choose the region to contain a fraction of 0.999 of the density and update this region at every MD step. After applying the decoherence correction eq , the changes of amplitudes of Ψ on all sites i outside of the active region, Δu, are set to zero and the amplitudes inside of the active region are scaled appropriately to preserve the norm. Hence, the decoherence correction is only applied locally within the active region by removing its unphysical long-range charge transfer components. Propagation of the wave function according to eq remains unaffected by the presence of the active region.

Our aim is to investigate the charge transport mechanism in the wide parameter regime that is characteristic of OSs. In model Hamiltonian studies, this is usually achieved by changing the numerical values for electronic and electron–phonon coupling in the Hamiltonian.[15,16,36] We apply a similar strategy in the current atomistic simulations. Here, electron–hole transport is simulated along linear 1D chains of ethylene-like molecules (ELMs), one of the simplest π-bonded extended systems. We refer to the molecules as “ethylene-like” because only their nuclear geometries correspond to real ethylene molecules. Properties determining the CT dynamics, such as reorganization free energy and electronic coupling strength of an ELM pair, are chosen as explained below to simulate hole transfer in different CT regimes.

In our model, hole transfer is mediated by a set of (orthogonalized) HOMO orbitals of the M ELMs comprising the chain, ϕ, i = 1,M. They are obtained from DFT calculations and used for calculation of H, k ≠ l = 1,M, with the AOM method.[29] The orbitals are propagated along the FOB-SH trajectories as explained in ref (22). The diagonal elements or site energies, H, are calculated with a force field, and λ was set to 200 meV in all simulations via (small) displacement of the C=C double bond equilibrium distance in the charged state.

FOB-SH hole transport simulations are carried out for three different electronic coupling strengths between the ELMs, giving average values of V = ⟨H2⟩1/2 ≈ 8, 30, and 120 meV along the trajectories. The different coupling strengths were obtained by scaling of the orbital overlap, and they are denoted in the following as “small”, “medium”, and “large”. We note that the CT parameters chosen here for the ELM chain are comparable to those of OS materials, where λ is typically between 100 and 200 meV and V between 1 and 100 meV.[9] Further details on the system setup and ELM parametrization can be found in the SI.

The steric effect of the extended environment is replaced by position restraints acting on the centers of mass of the ELMs. They prevent the chain from bending and keep the center of mass (COM) separations between the ELMs at approximately 4 Å, as they would in a condensed phase system. The ELMs are free to vibrate/rotate around their lattice site otherwise. While at low temperatures (100 K) the ELMs are preferentially oriented parallel, at room temperature (300 K) their orientation becomes disordered, and at the highest temperature investigated (1000 K) any orientation is about equally likely. Such behavior is typical of molecular OS materials and a consequence of the weak van der Waals interactions that hold the molecules together.

Initial configurations for FOB-SH simulations are prepared by thermally equilibrating the chain on diabatic potential energy surface H11. The hole carrier wave function is initially placed on ELM1, Ψ(0) = ϕ1(0), and the active adiabatic state ψ(0) is chosen randomly from the manifold of adiabatic states ψ(0), i = 1,M, with probability . The charge mobility is obtained by means of the Einstein relationwhere e is the elementary charge, kB is the Boltzmann constant, and D is the diffusion coefficient of the electron–holewhere MSD is the mean-squared displacement of the hole obtained from FOB-SH simulations. For each FOB-SH trajectory n, the center of the propagated hole wave function Ψ(t) (eq ) is calculated as a function of time and the resultant MSDs are averaged over all Ntraj trajectorieswhere uν, and xν, are, respectively, the expansion coefficients of wave function eq and the position of the COM of ELM ν in trajectory n and x1,(t = 0) = 0. Error bars for mobility were determined by block averaging of 1000 FOB-SH trajectories with a block size of 200 trajectories.

We first investigate the convergence of mobility with chain length and MD time step for the aforementioned different electronic coupling strengths; see Figure . We find that in all cases the mobility converges at a time step of 0.1 fs. The minimum chain length required for convergence increases with increasing coupling or mobility values from 5 ELMs for low coupling to 50 ELMs for high coupling strengths. The excellent convergence behavior is due to the state tracking algorithm applied for detection of trivial crossings and the projection algorithm applied to remove decoherence correction-induced artificial long-range charge transfers. When either one of the two is switched off, the mobility steadily increases and diverges with chain length (dotted and dashed–dotted lines in Figure ). This is because with increasing chain length the density of electronic states increases and, with it, the probability for trivial crossings and decoherence correction-induced artificial long-range charge transfer. The two algorithms that we employ here can deal with this problem very effectively, in particular, when considering that the results shown in Figure are for very high temperatures (1000 K), which is the most challenging regime as nuclear motion is fast and the probability to miss a trivial crossing is high. Moreover, we note that removal of the decoherence correction-induced artificial long-range charge transfer does not adversely affect the surface hopping dynamics as characterized by the internal consistency; see SI Figure 3. The average quantum amplitudes are, to a very good approximation, equal to the average surface populations at any point in time, just as in decoherence-corrected surface hopping without removal of artificial charge transfer.

Figure 1

Convergence of hole mobility with respect to size of the chain (in number of ELM molecules) and the MD time step. FOB-SH simulations were carried out for (A) large, (B) medium, and (C) small electronic coupling strengths V, with average values 120, 30, and 8 meV, respectively, for three MD time steps (0.5, 0.1, 0.05 fs). λ = 200 meV and T = 1000 K in all simulations. In panel (B), we show results if reordering of states via the state tracking algorithm (NO REORD) or decoherence-induced spurious charge transfer correction is switched off (NO SPTC).

We now come to the main result of this Letter: the hole mobility as a function of temperature, shown in Figure A for small, medium, and large electronic coupling strength. For small and medium coupling values, we observe thermally activated transport at low temperatures and band-like decay at higher temperatures. As the coupling strength is increased, the activated regime gradually disappears. For the largest coupling strength investigated, the mobility exhibits a band-like decrease for all temperatures according to a power law μ ∝ T–1.2. These results, obtained here for a fully atomistic model with electronic structure detail, are similar to the ones reported in the model Hamiltonian studies by Wang and co-workers,[15,16] implying that the latter provides a realistic coarse grain description for CT.

Figure 2

Temperature dependence of the (A) hole mobility and (B) inverse participation ratio (IPR, eq ) for small (red), medium (blue), and large (green) electronic coupling strengths V. The MD time step is 0.1 fs and λ = 200 meV in each case. FOB-SH mobilities eqs –8 are indicated by circles (○ symbols), and FOB-SH mobilities with off-diagonal elements H frozen to V are indicated by squares (“Frz. V”, □ symbols). Error bars are obtained from 5 blocks of 200 trajectories each. Hopping mobilities, obtained by solving a master equation with semiclassical ET rate constants, are indicated with solid lines for medium and low coupling strengths; see the SI for details. Best fits to an inverse power law are indicated by dashed lines for V = 120 meV. The resonance probability eq divided by T is shown with dotted lines. Snapshots of the polaron formed at 50 and 1000 K (indicated in (B) by blue and red arrows, respectively) are shown in Figure .

Figure 3

Snapshots of a FOB-SH simulation for a chain of 50 ELMs at large electronic coupling strength (V = 120 meV, λ = 200 meV) at (A) 50 and (B) 1000 K. The real part of the charge carrier wave function Ψ(t) eq is shown with red and blue isosurfaces and the imaginary part is shown with yellow and green isosurfaces. Note the increased localization of the polaron at higher temperature (panel B).

To obtain some insight into the molecular mechanism of charge transport, we characterize the localization length of the hole wave function by defining the inverse participation ratio IPR(t)The IPR measures the average number of molecules over which the wave function is delocalized. It converges rather quickly in our simulations, after 400 fs at high coupling and after 2 ps at medium and low coupling strengths (see SI Figure 1). The converged values, in the following simply denoted as IPR, are shown in Figure B.

For small coupling, IPR = 1 for all temperatures, indicating that hopping of a fully localized hole is the CT mechanism in this regime. Hole hopping from one molecule to the next is driven by thermal fluctuations in the site energies. When the energy of the site carrying the hole (say, site 1, H11) becomes resonant, i.e., within about H12 of the energy of the neighboring site 2, H22, the charge starts to oscillate between the two molecules and surface hops may take place. At the end of the lifetime of the resonance, the hole may settle on site 2, followed by thermal relaxation to a lower site energy as determined by reorganization energy λ. During this process, the active adiabatic electronic state that drives the nuclear dynamics remains very localized and changes from a state close to diabatic state ϕ1 to a state close to diabatic state ϕ2. The decoherence correction is paramount to achieve good internal consistency and localization of the hole wave function Ψ(t) ≈ ϕ2 after the hopping event has taken place. If the DC is switched off, Ψ(t) spreads unphysically, causing the IPR to diverge with time (see SI Figure 2).

The above observations suggest that CT in this regime may be well described by a simple hole hopping model, ELM1+–ELM2–ELM3– ELM1–ELM2+–ELM3– ···, with ET rates k obtained from ET theory and evaluated for the same values of electronic coupling and reorganization energy that were used in the FOB-SH simulations. To this end, we have solved a chemical master equation for hole hopping to obtain the time-dependent hole population of each site and the corresponding hole hopping mobility (see the SI for details). The results are shown in Figure A for the two coupling strengths where hole hopping is observed (low and medium coupling, solid lines). Agreement with the FOB-SH mobilities is very good, with deviations that are typically less than a factor of 3. The mobility saturation at around room temperature is because the ET activation free energy ΔA‡[23] becomes comparable to thermal energy at this point, ΔA‡ = (λ/4) – (V – V2/λ) ≈ 1–2kBT at T = 300 K and λ = 0.2 eV. For higher temperatures, the ET rate and diffusion coefficient become nearly temperature-independent (T–1/2 for nonadiabatic ET; T0 for adiabatic ET) and the mobility decay is dominated by the T–1 dependence of mobility on the diffusion coefficient in the Einstein equation (eq ).

The situation is strikingly different for large coupling strength where a finite activation barrier for hole transfer no longer exists (V > λ/2).[1,5,25] At low temperatures, the hole is delocalized over four ELMs on average, and as the temperature is increased, the IPR steadily decreases to 2. This decrease in IPR correlates remarkably well with the decrease in hole mobility. In Figure , we show a few snapshots of the hole wave function Ψ(t) along a randomly selected FOB-SH trajectory (red/blue and yellow/green isosurfaces indicate real and imaginary parts, respectively). At low temperature, the hole can quickly spread over several molecules to form a polaron that moves along the chain on the 100 fs time scale. By contrast, at high temperature, the hole wave function remains trapped over a smaller number of sites and cannot move as easily along the chain. Previously, localization of the polaron has been attributed to thermal fluctuations of electronic coupling (nonlocal electron–phonon coupling)[1,14,21] and site energy (local electron–phonon coupling).[15,16] To better understand the importance of these two contributions, we recomputed the mobility with the off-diagonal Hamiltonian matrix elements H frozen to the thermal average V. The results are shown in Figure (square symbols in green). We find that the qualitative trends in mobility and IPR are already well reproduced for frozen electronic couplings. Inclusion of coupling fluctuations enhances polaron localization and reduces mobility but does not change the qualitative picture.

In the following, we show that the trends observed for an M-site chain, the change in slope of mobility versus temperature at low couplings and the power law decrease at high couplings, correlate well with the temperature dependence of resonance in a simple two-site model. We consider the standard ET model comprised of two parabolic free energy curves for donor and acceptor diabatic states, A1(ΔE) and A2(ΔE), respectively, as a function of the site energy difference ΔE = H22 – H11 and assume constant electronic coupling V = H12. The free energy profile for the corresponding adiabatic electronic ground state can be written as[25]From the equation above, we can determine the corresponding probability of the site energy difference P(ΔE)Defining the ET resonance region (region of the transition state and avoided crossing) as −2V < ΔE < 2V, we obtain the probability to be in resonance, Pres, by solving the integral

The probability distributions eq are shown in Figure for small and large coupling strengths and for three different temperatures. The resonance probability eq is indicated by the shaded areas in each case, and Pres/T is superimposed on the hole mobility in Figure A (dotted lines). For small coupling, the peaks of the probability distribution are at ±λ for all temperatures, well outside of the narrow resonance region, explaining the low mobilities in the hopping regime (Figure A). Increasing the temperature widens the energy gap distributions and increases the resonance probability and therefore the hopping mobility. However, this effect saturates at around room temperature, causing the curve Pres/T to flatten and to cross over from a positive to a negative slope (Figure A (red and blue dotted lines)).

Figure 4

Probability density distributions eq for the site energy difference ΔE = H22 – H11 of a two-state model, λ = 200 meV. The resonance region for the charge transfer, ±2V, is indicated by black dashed vertical lines. (A) For small electronic coupling strength V ≪ λ/2, temperature broadening of ΔE increases the resonance probability Pres (eq , shaded area under curves). (B) For large coupling strength V ≥ λ/2, the temperature broadening reduces resonance probability. Pres/T correlates well with the computed mobilities for the ELM chain (dotted lines in Figure ).

The situation is markedly different at high coupling; see Figure B. In this case, the barrier for ET no longer exists and the transition state at ΔE = 0 becomes a minimum on the adiabatic free energy surface eq . Therefore, the peak of P(ΔE) is centered at ΔE = 0 for all temperatures, well within the resonance region. As the temperature is increased, the distributions again broaden, but this leads now to a reduction in resonance probability. The downward slope that the model predicts is in excellent agreement with the results from FOB-SH simulations; see Figure A (green dotted lines). It is remarkable that this very simple two-state model can reproduce the trends in mobility for an M-site chain in all relevant mobility regimes, bearing in mind that no dynamical effects are included, only equilibrium free energies. A natural extension of this model would be to include the possibility for M-site resonances (with M > 2)[25,42] and to derive an expression for the average drift velocity as required for prediction of mobilities.

Our results have implications for the screening of high-mobility crystalline OS materials:[21,43] (1) Electronic coupling should exceed half of the reorganization free energy (V > λ/2) along a connected path through the crystal lattice to ensure that CT is not in the temperature-activated hopping regime but mediated by mobile, delocalized polarons; (2) site energy fluctuations should be small (σ < V); and (3) electronic coupling fluctuations should be small (σ < V) to ensure that the molecules remain “resonant” for CT. In harmonic (linear response) theories, site energy fluctuations are related to reorganization free energy by σ = (kBTλ)1/2, implying that small λ’s favor both (1) and (2). Electronic coupling fluctuations are more difficult to predict as they depend on how much the HOMO orbital overlap between neighboring molecules (or the LUMO overlap for electron transport) changes in response to intermolecular vibrations, which in turn depends on the nodal shape of the orbitals and the amplitude of the vibrations. Ideally, one would lock neighboring molecules in a conformation with high electronic coupling (to fulfill (1) above) and restrict the amplitude of intermolecular vibrations, e.g., by introduction of bulky groups, to reduce electronic coupling fluctuations.

In conclusion, we have shown that FOB-SH is a powerful simulation methodology that can be used for prediction of charge carrier mobilities and CT mechanisms in strongly fluctuating media such as OSs. Based on an atomistic description, it naturally incorporates important effects determining electron–nuclear dynamics, e.g., local and nonlocal electron–phonon coupling. Several algorithmic advances on top of standard FSSH were necessary to arrive at this point: detection of trivial crossings via state tracking, electronic decoherence correction, and removal of spurious charge transfer caused by the latter. The method is now ready for benchmarking against experimental mobilities, which is currently being carried out in our laboratory.