Flickering Polarons Extending over Ten Nanometres Mediate Charge Transport in High‐Mobility Organic Crystals

Progress in the design of high‐mobility organic semiconductors has been hampered by an incomplete fundamental understanding of the elusive charge carrier dynamics mediating electrical current in these materials. To address this problem, a novel fully atomistic non‐adiabatic molecular dynamics approach termed fragment orbital‐based surface hopping (FOB‐SH) that propagates the electron‐nuclear motion has been further improved and, for the first time, used to calculate the full 2D charge mobility tensor for the conductive planes of six structurally well characterized organic single crystals, in good agreement with available experimental data. The nature of the charge carrier in these materials is best described as a flickering polaron constantly changing shape and extensions under the influence of thermal disorder. Thermal intra‐band excitations from modestly delocalized band edge states (up to 5 nm or 10–20 molecules) to highly delocalized tail states (up to 10 nm or 40–60 molecules in the most conductive materials) give rise to short, ≈ 10 fs‐long bursts of the charge carrier wavefunction that drives the spatial displacement of the polaron, resulting in carrier diffusion and mobility. This study implies that key to the design of high‐mobility materials is a high density of strongly delocalized and thermally accessible tail states.


DOI: 10.1002/adts.202000093
OS. Nowadays the charge mobility of organic single crystals routinely exceeds the one in amorphous silicon films (1 cm 2 V −1 s −1 for electrons). [1] This progress can be first and foremost ascribed to ever more sophisticated crystallization, purification, and device fabrication techniques. [2][3][4] Yet, designing materials with mobilities exceeding 10 cm 2 V −1 s −1 to enable a lower power consumption and a faster operation of organic circuits remains a central challenge. [5] One reason for this is that rational design approaches are still hampered by our incomplete understanding of the nature and transport mechanism of excess charge carriers in organic materials.
It is now well established that the weak molecular interactions between the molecules of OSs give rise to slow intermolecular vibrations [6][7][8] that strongly impact charge carrier mobility and device performance. This (dynamical) disorder also places charge transport in a difficult regime where the relevant transport parameters (electronic coupling, reorganization energy, site energy and electronic coupling fluctuations, and thermal energy) are all on the same energy scale, typically 10-200 meV. [6,7,[9][10][11][12][13] In this regime the standard transport theories, that is band theory on the one hand, and small polaron hopping on the other, cease to be valid and can no longer be physically justified, see, for example, refs. [6,13] for a detailed discussion. In addition, a plethora of experimental studies involving electron spin paramagnetic resonance, [14] Hall-effect measurements, [15,16] and charge-modulation spectroscopy [17] point to an intermediate transport scenario that was interpreted as "between" the well established band and hopping limits or as a coexistence regime of localized and delocalized charge carriers.
The impact of dynamic disorder in OSs [18] as well as the role of nuclear tunneling effects [19,20] have been studied before using golden-rule rates and related perturbative approaches, as they can be straightforwardly combined with efficient kinetic Monte-Carlo schemes for multiscale modeling. [18,21] However, perturbative methods (with a distinctive rate prefactor proportional to the square of the electronic coupling) assume that electronic coupling is much smaller than reorganization energy and they are therefore not strictly applicable to high mobility OSs including rubrene and pentacene. In addition, polaronic band theory www.advancedsciencenews.com www.advtheorysimul.com has been extended for OSs through inclusion of both local and non-local electron-phonon couplings at the level of model Hamiltonians. [22,23] The theory predicts a band narrowing with increasing temperature as a consequence of the larger polaron mass thereby successfully reproducing a band-like temperature dependence as often seen in experiments on ultrapure crystals, e.g. naphthalene. [24] However, like band theory, polaronic band theory also becomes problematic at higher temperatures. [6] More recently, Fratini and co-workers have developed the transient localization theory (TLT), [12,25] which is based on the observation that charge carriers are localized, to a certain degree, by thermal fluctuations of electronic couplings (non-local electron-phonon coupling) between the molecules. Asserting the relaxation time approximation, a simple relation was derived predicting that carrier mobility increases linearly with carrier delocalization. The theory was successfully used to show that charge transport in the high mobility planes of typical singlecrystalline OSs is enhanced if electronic couplings between the molecules within the plane are isotropic and are of specific sign combinations. [7,26,27] TLT has also recently been connected to the standard band transport description in the case of small disorder. [28] Yet, we note that TLT does not give information on how the charge carrier moves across the material in real time due to the assumptions of this theory with regards to the dynamics. Moreover, it is currently unclear how this theory is extended to the small polaron hopping regime where the mobility is no longer proportional to localization length. Hence, although it is extremely useful and important to have analytic theories, it is paramount to develop also numerical schemes that give insight into the actual dynamics and that seamlessly bridge the gap between different mechanistic regimes. With this purpose in mind, we have developed an efficient fully atomistic non-adiabatic molecular dynamics approach, denoted fragment orbital-based surface hopping (FOB-SH), [29][30][31][32] which allows us to propagate the coupled charge-nuclear motion in realtime for condensed phase systems. The methodology is based on a DFT-parametrized tight-binding representation of the electronic Hamiltonian (updated on-the-fly) to naturally incorporate local and non-local electron-phonon couplings, thus encompassing in a non-perturbative manner a broad range of possible transport mechanisms. In two previous applications, this method has proven very successful in predicting charge mobility in organic single crystals using 1D models for specific crystallographic directions [33] as well as for columnar crystal phases of tetracene derivatives. [34] In this work we further improve the efficiency of our FOB-SH methodology by using a multiple time step algorithm to reduce the computational cost without sacrificing accuracy. This allows us to effectively double the system size at the same computational effort and to converge the full 2D mobility tensor for the high conductivity plane of organic crystals. Another important aspect of this work is that we analyze in detail the mechanism that drives the polaron across the material resulting in charge mobility and electronic conductivity. We also investigate the nature of the intermolecular modes that give rise to off-diagonal disorder limiting charge transport and the effect on carrier delocalization when these modes are suppressed. In fact, it was recently pointed out in ref. [5] that a single or few so-called "killer" modes are responsible for a large suppression of the mobility for some OSs. We explore whether any such modes exist in the common organic crystals studied in this work.
Here we investigate six structurally well-characterized organic crystals including 1,4-bis(4-methylstyryl)benzene (pMSB), naphthalene (NAP), anthracene (ANT), perylene (PER), pentacene (PEN), and rubrene (RUB) covering three orders of magnitude in charge mobility (see Figure S1, Supporting Information for their crystal structure). We obtain near-quantitative agreement with experimental measurements for both magnitude and anisotropy of charge mobility, in particular for systems where time-of-flight (TOF) data are available. The charge carrier wavefunction is best described as a flickering polaron constantly changing its shape and extensions due to thermal intra-band excitations (though not to be confused with the flickering resonance mechanism of ref. [35] as explained further below).
In the following we briefly describe the FOB-SH methodology and refer to the Experimental Section for a more detailed explanation and simulation details. The electronic structure of typical organic crystals exhibits narrow bands across the entire Brillouin zone and a minimal band dispersion. Hence, the electronic Hamiltonian constructed in the space of the relevant frontier orbitals of the molecules forming the crystal (highest occupied molecular orbital (HOMO) for hole transport and lowest unoccupied orbital (LUMO) for electron transport), i , is expected to give a good approximation to the true band structure of the highest valence or lowest conduction band, see, for example, also refs. [6,9,11,12,36]. These considerations motivate the use of simplified one-particle electronic Hamiltonians for hole or excess electron transport, Equation (1) in Experimental Section, with diagonal and off-diagonal elements referred to as site energies and electronic couplings. Indeed, we find that the density of electronic states of this Hamiltonian is in good agreement with the results from standard band structure calculations at the Kohn-Sham density functional theory (DFT) level with regards to both peak positions and bandwidth (see Figure S3, Supporting Information).
In FOB-SH the charge carrier wavefunction, Ψ(t), is expanded in terms of the HOMO or LUMO orbitals, Equation (2), and propagated in time within the valence or conduction band by solving the time-dependent Schrödinger equation, Equation (3), with the electronic Hamiltonian Equation (1). Importantly, diagonal and off-diagonal electronic disorder are included due to the dependence of the electronic Hamiltonian on the nuclear degrees of freedom. Along with the charge carrier wavefunction, the nuclei are also propagated in time on the potential energy surface (PES) of one of the valence or conduction band states (i.e., eigenstates) of the electronic Hamiltonian Equation (1) (referred to as active eigenstate a), and intra-band transitions of the nuclear dynamics ("hops") to another band state j are included using Tully's surface hopping probability. [37] As described in detail in our previous work, [31,32] we use FOB-SH in combination with three important extensions of the original surface hopping method: decoherence correction, removal of decoherence correction induced artificial long-range charge transfer and tracking of trivial surface crossings. These algorithmic advances improve a number of desirable properties of FOB-SH including Boltzmann occupation of the band states in the long time limit, internal consistency between charge carrier wavefunction and surface populations of  The experimental polaron size for pentacene, as obtained from electron spin paramagnetic resonance, is depicted in green dashed lines. [14] the band states, and convergence of charge mobility with system size and nuclear dynamics time step. [32] We refer to the Experimental Section for a specific discussion of the multiple time step algorithm as well as a more efficient propagation of the electronic Schrödinger equation introduced in this work to deal with large systems.
The charge carrier wavefunction is initially chosen to be fully localized on a single molecule i, Ψ(t = 0) = i , and propagated in time according to the FOB-SH algorithm. The mean-square displacement (MSD) for the centre of charge of Ψ(t) in the a-b high-mobility plane of PER-e − and RUB-h + is shown in Figure 1a,b, where the appendix "e − " and "h + " indicates electron and hole transport, respectively (see Figure S4, Supporting Information, for the other OSs.) At short times the MSD tends to increase rapidly until after a few 100 fs the increase is linear in time within the error bars of our simulations. The initial dynamics is due to the relaxation of the fully localized charge carrier, which is a mixture of excited band states, to states close to the valence or conduction band edge. During this quantum relaxation process the charge carrier "equilibrates" to form a polaron within the a-b plane accompanied by a non-linear increase in the MSD. A similar relaxation dynamics was found by Schnedermann et al. [38] by using optical transient absorption spectroscopy in the context of exciton and charge transport in a thin pentacene film. As we describe further below the polaron is a highly dynamical species constantly changing shape and extension. A representative snapshot is depicted in Figure 2 where the carrier delocalization is larger than ten molecules in the most conductive materials. We verified that the average polaron size and shape is insensitive to the choice of the initial carrier wavefunction while the time for relaxation to this state may, of course, vary. This is because FOB-SH maintains detailed balance in the long-time limit to a good approximation [30,32] ensuring that after initial relaxation, the equilibrium (Boltzmann) populations of electronic band states are reached, independently on the initial starting point. It should be noted that the finite size of the polaron is due to the thermal diagonal and off-diagonal disorder in the electronic Hamiltonian Equation (1). Without the disorder the carrier wavefunction would be a delocalized band state.
Returning to the MSD, the linear regime that follows at longer times is indicative of Einstein diffusion of the polaron and this is the part that is used for the determination of the diffusion and charge mobility tensors according to Equations (4)- (6). As shown in Figure 1c, system sizes of almost 1000 electronically active molecules within the a-b plane are necessary to converge charge mobility for the most conductive materials, while the polaron size defined by the inverse participation ratio (IPR), Equation (7), converges significantly faster, Figure 1d. On a more technical note, we point out that the good convergence of mobility and IPR with system size (Figure 1c,d) as well as with time step ( Figure S5, Supporting Information) implies that www.advancedsciencenews.com www.advtheorysimul.com spurious long-range charge transfer due, for example, to trivial crossings have been largely eliminated.
The computed charge mobilities for the a-b planes are shown in polar representation in Figure 2. As one might expect and in agreement with transient localization theory, the mobility in a given direction correlates well with the delocalization of the polaron in that direction as both depend primarily on the strength of electronic coupling between the molecules. Overall, we observe good to near quantitative agreement between FOB-SH (data in blue) and experimental mobilities (data in black) for both magnitude and anisotropy within the plane, especially when comparing to Karl et al.'s time-of-flight (TOF) data for NAP, ANT, and PER (panels b, c, and d, relative error averaged over a and b direction is 24%, 31%, and 30%). This technique probes the bulk mobility at dilute carrier concentration and corresponds most closely to the conditions for which FOB-SH simulations are carried out.
Unfortunately, to the best of our knowledge, TOF data are not available for pMSB, RUB, and PEN. Here comparison must be made to field-effect-transistor (FET) data, which are more characteristic of surface mobilities and may depend on the nature of the gate insulator. With these caveats in mind, we find good agreement between the computed highest in-plane mobility for RUB, 16 cm 2 V −1 s −1 along the a-direction, and the most reproducible experimental FET measurements, in the range 10-20 cm 2 V −1 s −1 . [39][40][41] The anisotropy is also reasonably well captured albeit somewhat underestimated. The computed mobility along the bdirection is 10 cm 2 V −1 s −1 compared to 4-8 cm 2 V −1 s −1 from experiment. [39][40][41] This shows that RUB should be considered as a 2D conductor with preferential conduction along the a-direction, not 1D as sometimes described in the literature. At first sight, this may be somewhat surprising because rubrene forms a slipped pistacked structure with electronic couplings along the b-direction (T 1 and T 2 ) about a factor of 4 smaller than along the a-direction (P). However, the thermal fluctuations of the electronic couplings T 1 and T 2 are also a factor of 3 smaller than for P (see Table S1 and Figure S2, Supporting Information). We believe that this partly offsets the smaller mean couplings leading to comparable mobilities in a and b directions. A discussion of our results for PEN can be found in the SI.
It is worth noting that all mobilities reported here from FOB-SH simulation for the full 2D plane are larger than the mobilities obtained previously [33] for 1D models along the same crystallographic direction, typically by a factor of 2-3 (see Table S3, Supporting Information, for a summary). This is in line with the finding of Fratini et al. who concluded that 2D systems with isotropic couplings exhibit higher mobilities than anisotropic systems and noting that the 1D models are perfectly anisotropic. [26] Intuitively, the charge carrier wavefunction can delocalize in a 2D plane more effectively than in a 1D chain because of the larger number of nearest neighbours. Therefore, when possible bottlenecks are present (e.g., neighboring molecules in a pose with unfavorable couplings), the charge carrier may bypass them in other directions. [21] Nevertheless, rapid evaluations of mobilities from 1D models can already provide good estimates of the intrinsic mobilities of the system at hand. [33,42] We discuss briefly the performance of alternative approaches including hopping, band theory, and TLT. In pMSB, NAP and ANT the electronic couplings are smaller than half the reorganization energy (see Table S1 and Figure S2, Supporting Information), hence a finite barrier for site-to-site hopping of a fully localized polaron exists and electron transfer (ET) rates may be defined. The mobility obtained by solving such a hopping model is shown in Figure 2 (green lines), see SI for details. It gives relatively accurate results for the OSs with the lowest mobilities, pMSB and NAP, but significantly overestimates the mobility for ANT. However, it is clear from the FOB-SH simulation that even in the system with the lowest mobility (pMSB) the charge carrier does not simply transfer via site-to-site hopping as in donor-acceptor electron transfer reactions. Moreover, as pointed out before, [26] since the hopping rates depend on the magnitude of the couplings, the effect of the sign relationship between the different couplings on the symmetry of DOS is not accounted for in this approach. At the other extreme, standard band theory calculations tend to strongly overestimate mobility in all systems, as discussed previously. [33] Mobilities from transient localization theory (TLT) are reported in Figure S6, Supporting Information, and compared to mobilities from FOB-SH and experiment (see Supporting Information for details). We note that the TLT mobilities are in excellent agreement with FOB-SH and experiment but only if the diagonal disorder is excluded in the calculation of the localization length (green symbols). When diagonal disorder is included to ensure a likefor-like comparison with FOB-SH, TLT mobilities tend to underestimate FOB-SH and experimental mobilities (red symbols). The deviation tend to become larger, in some cases up to an order of magnitude, for the OSs with the lowest mobilities. This is expected because TLT, at least in its present form, does not extend to the low mobility/strong localization regime. Another interesting point is that the average IPR (with which we quantify the extension of the polaron in a statistical manner) generally correlates well with the localization length. This was pointed out by Fratini et al. [26] and further investigated by us for real systems in our previous work. [33] However, it is important to mention that the IPR per se does not give dynamical information while the localization length is directly related to charge mobility in TLT.
Having quantified the computed charge mobilities we now turn our attention to the mechanism of the diffusive charge transport: how do the polaronic charge carriers shown in Figure 2 move within the organic semiconductor? Watching the real-time trajectories generated by FOB-SH, the charge Charge mobilities from FOB-SH (blue) and from a site-to-site hopping model using electron transfer rates (green) are compared to experimental mobilities (black). For NAP, [60,61] ANT [62] and PER, [63,64] time-of-flight (TOF) data along a and b directions are indicated by black circles (•). For pMSB, the field-effect transistor (FET) mobility is taken from Ref. [65]. For RUB, FET mobilities are taken from Ref. [39] (▶), Ref. [66] (◀) and Ref. [40] (▴). For PEN, FET mobilities are taken from Ref. [67] (▶) and Ref. [68] (▴) and SCLC mobilities from Ref. [  carrier wavefunction Ψ(t) may be best described as a flickering polaron continuously changing its shape and extensions (see animations provided in the Supporting Information) though with preference to delocalize in the direction with the highest coupling. Remarkably, we frequently observe events where the polaron expands to about 2-3 times its average size, followed by relocalization at a position about a polaron diameter apart from the original position (Figure 3c-e). These "diffusive jumps," as we termed them in ref. [33] in analogy to molecular diffusion in heterogeneous media, [43] are at the origin of charge mobility in OSs. They occur on the 100 fs time scale for all OSs investigated (similar to the parameter of TLT [26] ) and the short wavefunction expansions preceding the jumps, denoted "resonances" in ref. [33], last for about 10 fs. [33] Here we analyze the origin of the resonances and diffusive jumps in RUB and NAP in more detail, from the perspective of the charge-nuclear dynamics in the valence or conduction band, as simulated by FOB-SH. In essence, we will show that short-lived thermal intra-band excitations to delocalized band states drive the dynamics of the charge.
At first we consider the potential energy, E a , of the active valence band state, a , on which the nuclear dynamics is run in FOB-SH, (lines in red in Figure 3a,b, top of valence band at 0 meV). We observe frequent surface hopping events (dashed black lines) that ensure approximate Boltzmann population of the active valence band state on the fast ps time scale of present simulations. The average electronic energy is ⟨E a ⟩ ≈ −1.5 k B T (T = 300 K) and low-lying valence band states up to ≈ −200 meV below the top of the valence band are occupied for very short durations of time. Notably, there is a good correlation between energy and delocalization of the band states: the lower E a the higher the IPR of the active state, IPR a (blue lines). This is in accord with our analysis of the DOS ( Figure S3, Supporting Information) and the well-known fact that the states in the middle of the valence band are more delocalized than at the top of the valence band.
The actual hole carrier wavefunction of RUB and NAP, Ψ(t), closely tracks the active valence band state of the system, Ψ(t) ≈ states interact. Consequently, the IPR of Ψ(t) (blue lines in the middle panels of Figure 3a,b) also closely follows IPR a , though we observe that the former is smoother and peaks with a certain delay compared to the latter due to the finite decoherence time. Crucially, the IPR of Ψ(t) correlates very well with the instantaneous rate of charge carrier displacement or classical drift velocity v d a (t)(v d b (t)) contributing to the MSD (green lines, subscript denoting displacement along the a-direction (b-direction) of RUB (NAP)): the peaks in the IPR give rise to peaks in the instantaneous drift velocity in about 75% of cases -these are the productive resonances resulting in charge carrier displacement and mobility, as exemplified in Figure 3c-e at about 200 fs. There are also unproductive resonances where a peak in the IPR does not lead to a corresponding increase in drift velocity, for example, at 160 and 260 fs in Figure 3b for NAP). This happens, for example, when the polaron undergoing an expansion returns to its original position or when the expansion is near symmetric around the original position so that the center of charge does not move.
We note that the motion of the flickering polaron described here is the result of explicitly solving the time-dependent electronic Schrödinger equation coupled to nuclear motion and should not be confused with the flickering resonance model proposed by Skourtis and Beratan in ref. [35]. In this theoretical model mechanism the authors consider electron transfer (ET) between a pair of molecular (i.e., fully localized) donor (D) and acceptor (A) levels bridged by intermediate "bridge" levels (B) and assert that donor-acceptor ET occurs when all levels A, D, and B come into energetic resonance with one another. This model still assumes that reorganization energy is significantly larger than electronic coupling so that the ET is in the non-adiabatic or adiabatic regime. Indeed, for low mobility OSs with sufficiently small coupling the transport mechanism obtained from FOB-SH resemblances closely the flickering resonance model. [33] However, for medium to high mobility OSs the couplings are too large for the flickering resonance rate model to be valid, the polaron is permanently delocalized over several molecules and a clear dis-tinction between donor, bridge, and acceptor sites is no longer possible.
As pointed out before, the delocalization of the polaron and hence its mobility is limited by the detrimental effect of diagonal and off-diagonal electronic disorder. [10,44] While the diagonal disorder is due to intramolecular vibrations, the off-diagonal disorder is due to intermolecular vibrations and several attempts have been made to suppress them, for example, by introduction of bulky side chains or by inducing strain to limit the amplitude of these fluctuations. [45][46][47] Moreover, Schweicher et al. recently reported that in one OS, C8-DNNT-C8, 75% of the off-diagonal disorder is due to a single "killer mode" associated with a longaxis sliding motion. [5] It is therefore pertinent to explore whether such "killer modes" also exist in the OSs presented in this work. To this end we go one step further and calculate the Boltzmann and time averaged IPR of the band states, ⟨IPR⟩ (see Equation (10) in the Supporting Information), when all off-diagonal thermal fluctuations larger than a cutoff-frequency max are removed. The change in the ⟨IPR⟩ is then a proxy for the expected change in mobility.
The ⟨IPR⟩ as a function of max is shown in Figure 4, see Supporting Information for details of the calculation. For comparison, the average IPR of Ψ(t) obtained from FOB-SH runs with all frequencies included is shown in dashed grey lines. For OSs with medium to high mobilities, such as PER, PEN, and RUB, we find that ⟨IPR⟩ indeed strongly increases as coupling fluctuations slower than 50 cm −1 are removed, for example, Figure 4a (see Figure S7a, Supporting Information, for state-resolved IPR). Hence, removing coupling fluctuations for these materials will increase the transport efficiency. Yet, the increase in ⟨IPR⟩ is smooth and appears to be a collective effect of many modes -it is not possible from this analysis to single out a small set of discrete phonon frequencies that would be particularly important in increasing the delocalization. We also note that even when all coupling fluctuations are filtered out ( max = 0) the charge remains localized over a finite number of molecules due to diagonal disorder.

www.advancedsciencenews.com www.advtheorysimul.com
The situation is qualitatively different for the low mobility OSs pMSB and NAP: ⟨IPR⟩ now slightly decreases as the coupling fluctuations are removed, e.g. Figure 4b. In these materials the magnitude of electronic coupling is small and the low frequency fluctuations increase, but only very slightly, the delocalization of these states (see also Figure S7b, Supporting Information) thereby favoring formation of resonances and facilitating charge transport. [11] If one applies Marcus theory for these OSs (with the above caveats), this trend becomes immediately obvious: the hopping rate is proportional to ⟨|H kl | 2 ⟩ = ⟨H kl ⟩ 2 + 2 and therefore increases with increasing off-diagonal thermal fluctuations 2 . This analysis shows that the mobility of high-mobility OSs can be further increased by removing off-diagonal thermal disorder, but that low mobility OSs cannot be turned into high-mobility OSs via this strategy.
Our work shows that key to high-mobility OSs is a high density of thermally accessible delocalized states at the top of the valence band (hole) or bottom of the conduction band (electron). For OSs this is only the case when several conditions are met 1) electronic couplings larger than half the reorganization energy, V > ∕2, where V = √ ⟨|H kl | 2 ⟩ is the mean coupling, to avoid trapping and formation of small polarons, 2) small thermal fluctuations of electronic coupling (i.e. low off-diagonal disorder), < 0.3V (where 0.3 corresponds to the value for rubrene which is, somewhat arbitrarily, chosen as a reference here). 3) For 2D and 3D conduction, isotropic couplings with specific sign combinations favoring lowenergy delocalized electronic states. Importantly, we note that the same consensus on design rules has been reached before on the basis of transient localization theory. [26] The problem with most OSs is that one or more of these requirements are not fulfilled. While many OSs exist where (1) is fulfilled, it is still an open question how to best reduce offdiagonal disorder (2) without simultaneously diminishing the mean electronic coupling V. As mentioned before, reduction of thermal fluctuations has been attempted with various degree of success via core functionalization of organic molecules [5,45] and by application of external crystal strain and pressure. [46][47][48] Even more challenging is the design of specific sign combinations of couplings due to the complicated nodal shape of the relevant frontier orbitals. Consequently, the thermally accessible states in most materials tend to be rather localized resulting in modest instantaneous drift velocities and mobilities, as illustrated for NAP in Figure 3b. Arguably, among all known molecular OSs, RUB fulfils these criteria best, though its mobility would be even greater if the couplings were more isotropic and the off-diagonal electronic disorder was smaller. Nevertheless, the quantities appearing in rules (2)-(3) can be relatively straightforwardly calculated from DFT and used to refine high throughput screening studies that have previously focused on rule (1) alone. [49,50] In conclusion, we have reported the full 2D charge mobility tensors for six organic crystals and uncovered the real-time dynamics of the charge carriers using a powerful non-adiabatic molecular dynamics simulation methodology. We find that the charge carrier wavefunction forms a flickering, highly dynamic polaron that is delocalized over about 5 nm on average in the most conductive crystals and of finite size due to thermal energetic disorder. Thermal intra-band excitations lead to short, ≈ 10 fs-long bursts of the polaron during which it expands to 2-3 times of its average size, followed by deexcitation and relocaliza-tion. It is these short bursts that drive charge carrier diffusion in these materials. Hence, from this dynamical perspective, it is more suitable to describe charge carrier transport in OSs as a transient delocalization (rather than transient localization) process. The challenge for the future is to design stable materials that exhibit the three characteristics referred to above.

Experimental Section
FOB-SH: The valence (or conduction) band of the OSs is described by the following Hamiltonian, where, k = k (r, R(t)) is the HOMO (LUMO) of molecule k for hole (electron) transport, r is the position of the hole or excess electron, R(t) are the time-dependent nuclear coordinates, k = k (R(t)) is the site energy, that is, the potential energy of the state with the hole (excess electron) located at site k and H kl = H kl (R(t)) is the electronic coupling between k and l . All Hamiltonian matrix elements, that is, site energies and couplings, depend on the nuclear coordinates which, in turn, depend on time, R = R(t) as determined by the nuclear dynamics (see below). In the FOB-SH approach the hole (excess electron) is described by a time-dependent 1-particle wavefunction, Ψ(t), expanded in the same basis that is used to represent the Hamiltonian Equation (1), where u l are the expansion coefficients. Insertion of Equation (2) in the time-dependent Schrödinger equation gives the time-evolution of the charge carrier wavefunction in the valence (conduction) band, where d kl = ⟨ k |̇l⟩ are the non-adiabatic coupling elements. The nuclear degrees of freedom are propagated on one of the potential energy surfaces (PES) obtained by diagonalizing the Hamiltonian Equation (1) and denoted as E a ("a" for "active surface"). While nuclear motion couples to the motion of the excess charge via the dependences on R(t) in Equation (3), feedback from the excess charge to the nuclear motion is accounted for by transitions of the nuclear dynamics ("hops") from the PES of the active eigenstate a to the PES of another eigenstate j using Tully's surface hopping probability. [37] As pointed out in the introduction, the FOB-SH algorithm features several important improvements over the original surface hopping methodology [37] as necessary for robust and meaningful mobility calculations. Here we only describe in some detail the decoherence correction used and refer to ref. [32] for a detailed description of trivial crossing detection and elimination of spurious long-range charge transfer algorithms. The decoherence correction is based on exponential damping of all except the active band states (j ≠ a): c j → c j exp(−Δt∕ ja ). [51] The coefficients c j are the expansion coefficients of the charge charge carrier wavefunction Ψ(t) in terms of the eigenstates, j , of the electronic Hamiltonian, Ψ(t) = ∑ j c j (t) j . The coefficient for state a, c a is scaled appropriately to ensure norm conservation. We adopt, here, the low-cost and parameterfree Heisenberg decoherence time: ia = ℏ∕|E i − E a |, and note that other common choices, subjected to different physical arguments such as the force-based decoherence time, [52] give very similar results for charge mobility, as long as the decoherence time is fast enough to maintain internal consistency. [31] www.advancedsciencenews.com www.advtheorysimul.com Mobility Calculation and IPR: Solving Equation (3), one obtains the charge carrier wavefunction as a function of time, Ψ(t). This gives access to key dynamical properties, for example, the mobility tensor (Equation (4)), the extent of localization or delocalization of the charge carrier as a function of time and the mechanism by which the charge carrier moves within the material. The charge mobility can be expressed as a second rank tensor using the Einstein relation, where ( ) represent Cartesian coordinates, x, y, z. e is the elementary charge, k B the Boltzmann constant and T the temperature (T = 300 K in this work). The diffusion tensor components, D , can be obtained as the time derivative of the mean squared displacement along the nine Cartesian components (MSD ), where, In Equation (6) Figure 1; Figure S4, Supporting Information).
A common measure was used to describe the delocalization of the charge carrier wavefunction Ψ(t), the inverse participation ratio (IPR), The numerical value of the IPR is about equal to the number of molecules the wavefunction is delocalized over. A similar definition is used to describe the delocalization of the adiabatic states or eigenstates of the Hamiltonian Equation (1), i (t), where U ki,n are the components of the eigenvector i of the Hamiltonian Equation (1) in trajectory n (see Figure 1; Figure S5, Supporting Information, for convergence).

Site Energies and Electronic Coupling:
The force field parameters for calculation of the site energies of the Hamiltonian Equation (1), k , were taken from previous work. [33] In particular, the force field parameters for the charged state was parametrized to reproduce the intramolecular reorganization energy from DFT calculations: 255, 187, 142, 177, 152, 98 meV for pMSB-h + , NAP-h + , ANT-h + , PER-e − , RUB-h + , and PEN-h + , respectively. Such reorganization energies are typical for organic semiconductors and an order of magnitude smaller than those, for example, for redox processes in aqueous solution [53,54] or oxide materials. [55,56] Coupling (or off-diagonal) matrix elements, H kl , were obtained for molecular dimers forming the crystal using our well-established analytic overlap method (AOM). [57] This method is based on the observation that for conjugated molecules the coupling depends linearly on the orbital overlap, to a good approximation, H kl = CS kl where C is a constant determined from DFT calculations (sFODFT, giving C < 0 for all systems) andS kl is the overlap between HOMO (LUMO) orbitals on molecules k and l projected on a minimum Slater basis. For excess electron transport (PER-e − ), the site energies obtained from the force field correspond to electronic energy levels, hence the sign of C is the same sign as in DFT calculations, C < 0. For hole transport (all systems except PER-e − ), the site energies obtained from the force field correspond to hole energy levels and therefore the sign of C obtained from DFT calculations is inverted, C → −C. Consequently, in FOB-SH the energy levels for both excess and hole transport are approximately populated according to exp[−E i ∕(k B T)] where E i are the excess electron or hole energy levels (see also ref. [42]). Further discussion of the importance of the correct sign of H kl is given in the Supporting information. Reference sFODFT electronic couplings, AOM electronic couplings and C values are summarized in Table S1, Supporting Information, for all systems investigated. One should note that the good agreement in shape and bandwidth of the DOS of the FOB-SH Hamiltonians compared with KS-DFT ( Figure S3, Supporting Information) attests to the reliability of the Hamiltonians used in this study.
Force Calculation: The calculation of the nuclear forces on nucleus I in a given adiabatic state i and nucleus I, F I,i , was obtained from the gradient of the Hamiltonian matrix elements in the diabatic representation using the Hellmann-Feynman theorem, (9) where [∇ I ℍ] kl ≡ ∇ I H kl = ∇ I ⟨ k |H| l ⟩. The last identity in this equation has been shown explicitly in ref. [29]. In practice, the gradients of the diagonal elements (∇ I H kl with k = l) are obtained by using the gradient of classical force-field potentials, while the off-diagonal elements are found by finite differences of the orbital overlap that comes from the analytic overlap method (AOM), namely ∇ I H kl = C∇ ISkl . [57] As the number of atoms in the system increases and the off-diagonal elements to evaluate becomes larger, the calculation of the off-diagonal gradients of the Hamiltonian becomes the time-limiting step. Here, a multiple time step algorithm (MTS) was introduced to reduce the computational cost. In particular, all the gradients ∇ I H kl , with k ≠ l, are updated only every N MD time steps and kept unchanged between two updates. N must be chosen small enough to reproduce the time oscillations of the off-diagonal gradients well. Since the electronic couplings in OSs generally fluctuate with an oscillation period of ≈ 1 ps, [58] one can expect the gradients of the couplings to oscillate on the same time scale. It is worth mentioning that similar MTS approaches are often used in MD codes to efficiently speed-up different parts of the computation. [59] The quality of this algorithm is assessed in Figure S9, Supporting Information. It was noted that the same approach cannot be applied to the diagonal gradients without biasing the whole dynamics, as the site energies fluctuate in the order of the aromatic carbon stretching frequencies (≈ 20 − 30 fs).

Non-Adiabatic Coupling Element in the Localized Basis:
The nonadiabatic coupling elements (d kl ) in the localized orthogonal orbital basis that appear in Equation (3) would need to be evaluated together with the Hamiltonian in order to propagate the electronic equation of motion. In principle, d kl can be evaluated using the AOM approach as done in our previous works, [29][30][31][32] and detailed in the Supporting Information. It was found that calculating this term at each nuclear time step (and then linearly interpolating it at each electronic step while integrating Equation (3)) gives essentially the same dynamics as neglecting it completely in the propagation equation (see mobility and IPR in Figure S9, Supporting Information). This is somewhat expected as, in practice, d kl is always small (typically below 0.04 meV ℏ −1 for the investigated OSs) and smooth along the entire dynamics since the localized orthogonal basis { k (t)} is, in fact, quasidiabatic. Moreover the second term on the RHS of Equation (3) including d kl is, on average, only 0.5% of the first term including electronic couplings. Hence, for most practical purposes the term including d kl can be safely neglected. Importantly, avoiding the calculation of d kl means bypassing a matrix-matrix-matrix multiplication (see Supporting Information) at each nuclear time step and it permits a speed-up of almost a factor of 1.5 compared to the usual interpolation scheme when the system size reaches more than a thousand molecules. This efficient optimization combined with the MTS algorithm allowed one to almost double the system size at the same computational cost without sacrificing accuracy of mobility calculation (see details in Supporting Information and related benchmark in Figure S9, Supporting Information).
Simulation Details: For each OS, a series of supercells of increasing size were built from the experimental crystallographic unit cell. The dimensions of the largest supercells constructed are summarized in Table S2, Supporting Information. These supercells were equilibrated in periodic boundary conditions for the neutral state for at least 250 ps in the NVT ensemble to a target temperature of 300 K using a Nosé-Hoover thermostat, followed by at least 250 ps equilibration in the NVE ensemble. [32,33] From the NVE trajectory an uncorrelated set of positions and velocities were chosen as starting configurations for FOB-SH simulations. Molecules within a rectangular region of the a − b high mobility plane were treated as electronically active, that is, as molecular sites or fragments for construction of the electronic Hamiltonian (Equation (1)), with their frontier orbital (HOMO or LUMO) contributing to the expansion of the carrier wavefunction (Equation (2)). All other molecules of the supercell were treated electronically inactive and interacted with the active region only via non-bonded interactions. The initial carrier wavefunction is chosen to be localized on a single active molecule m, Ψ(0) = m and propagated in time according to the FOB-SH algorithm in the NVE ensemble. All FOB-SH simulations applied a decoherence correction, state-tracking for detection of trivial crossings, a projection algorithm for removal of decoherence correction-induced artificial long-range charge transfer and adjustment of the velocities in the direction of the non-adiabatic coupling vector in case of a successful surface hop. The nuclear time step, Δt, ranged from 0.01 to 0.1 fs depending on the size of the systems (see Figure S5, Supporting Information). The electronic time step for integration of Equation (3) using the Runge-Kutta algorithm to 4th order was t = Δt∕5. For each system at least 300 FOB-SH trajectories of length 1 ps were run. The components of the diffusion tensor Equation (6) were block averaged over three blocks (at least 100 trajectories each) for calculation of error bars. Here it was noted that for all crystals except pentacene (PEN) the Cartesian coordinates (x, y) of the supercell were chosen parallel to the crystallographic directions (a, b) that define the high mobility plane (this was possible since they form either orthorhombic (pMSB, RUB) or monoclinic crystals (NAP, PER, ANT)). In this representation the off-diagonal components of the diffusion tensor are zero due to symmetry. For pentacene (triclinic) the diffusion tensor was diagonalized. The number of active molecules required for convergence of charge mobility and the largest number of active molecules considered for each OS are summarized in Table  S2, Supporting Information. All simulations were carried out with our in-house implementation of FOB-SH in the CP2K simulation package. [59]

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.