On dynamical localization corrections to band transport

Bloch-Boltzmann transport theory fails to describe the carrier diffusion in current crystalline organic semiconductors, where the presence of large-amplitude thermal molecular motions causes substantial dynamical disorder. The charge transport mechanism in this original situation is now understood in terms of a transient localization of the carriers' wavefunctions, whose applicability is however limited to the strong disorder regime. In order to deal with the ever-improving performances of new materials, we develop here a unified theoretical framework that includes transient localization theory as a limiting case, and smoothly connects with the standard band description when molecular disorder is weak. The theory, which specifically adresses the emergence of dynamical localization corrections to semiclassical transport, is used to determine a"transport phase diagram"of high-mobility organic semiconductors.

Bloch-Boltzmann transport theory fails to describe the carrier diffusion in current crystalline organic semiconductors, where the presence of large-amplitude thermal molecular motions causes substantial dynamical disorder. The charge transport mechanism in this original situation is now understood in terms of a transient localization of the carriers' wavefunctions, whose applicability is however limited to the strong disorder regime. In order to deal with the ever-improving performances of new materials, we develop here a unified theoretical framework that includes transient localization theory as a limiting case, and smoothly connects with the standard band description when molecular disorder is weak. The theory, which specifically adresses the emergence of dynamical localization corrections to semiclassical transport, is used to determine a "transport phase diagram" of high-mobility organic semiconductors.

I. INTRODUCTION
The last decade has witnessed considerable progress in the understanding of charge transport in high-mobility organic semiconductors, with important milestones achieved in both experimental and theoretical research. On the experimental side, the widespread access and improved control on field-effect devices has provided a common ground for the systematic and reproducible measurement of carrier mobilities 1 . Initially restricted to crystalline rubrene 2,3 , which served as a prototypical material due to its outstanding performances and stability, results indicative of intrinsic charge transport are now obtained in a growing number of organic semiconductors [4][5][6][7][8][9][10][11] . On the theoretical side, it is now understood that the dominant intrinsic factor limiting the mobility is the presence of large thermal vibrations of the constituent molecules, which cause strong dynamic disorder 12 . The latter hinders the carrier motion in ways that differ substantially from what predicted by semiclassical scattering theories.
The dynamical nature of molecular disorder in organic solids makes the quantum theory of Anderson localization, that was developed for systems with static randomness 13,14 , of limited use per se. To deal with this original situation, the physical idea of charge carriers being coherently localized, but only over a limited timeframe, has emerged over the years 12,15,16 . By highlighting explicitly a connection with the physics of localized systems, the concept of transient localization could reconcile a number of puzzling features of organic semiconductors, most notably the observation of "bandlike" mobilities decreasing with temperature to values below the so-called Mott-Ioffe-Regel limit. The applicability of the transient localization approach, however, is by construction restricted to the regime of strong dynamical disorder. Materials with reduced localization effects, either featuring more isotropic two-dimensional band structures 17 or lower degrees of disorder 18 , are actively investigated. It can be expected that future organic compounds will progressively move away from the strong disorder regime, entering a crossover region for which there is yet no available theoretical description.
The aim of this work is to establish a unified theoretical framework that encompasses the whole range from Bloch-Boltzmann band theory, that applies in the limit of weak electron-phonon scattering [19][20][21][22] , to the transient localization (TL) regime relevant when dynamic disorder is strong 16 . Our approach, which is based on the evaluation of dynamical localization corrections to semiclassical transport, is valid regardless of the disorder strength, as confirmed by the comparison with available exact numerical data. We illustrate our findings on organic compounds of current interest, determining a general "transport phase diagram" for high-mobility organic semiconductors.

A. Dynamical localization corrections
Due to its semiclassical nature, Bloch-Boltzmann theory neglects localization processes altogether. While this is a viable approximation to treat electron-lattice interactions in the weak scattering limit, it is inappropriate to study the transport properties of electrons in strongly fluctuating environments, as is the case in organic semiconductors owing to the presence of large-amplitude molecular motions. Here we want to overcome this limitation by restoring those quantum processes that are missing in the semiclassical description. Our derivation is based on Kubo response theory, and builds on the formulation developed in Refs. 15,23 . The key quantity of interest is the time-dependent velocity-velocity anticommutator correlation function, C(t) = {V(t),V(0)} , withV the velocity operator for charge carriers in a given direction. The dynamical observables describing charge transport, i.e. the charge diffusivity, the charge mobility and the optical conductivity, can all be derived from the knowledge of this time-dependent quantity 15,23 .
Let us imagine that for a given system of electrons interacting with lattice vibrations, we are able to calculate the velocity correlation function using semiclassical Bloch-Boltzmann theory (see Appendix A), that we denote as C S C (t). By definition this quantity misses all those quantum processes that are instead present in the exact correlator C(t). It is then natural to define as dynamical localization corrections (DLC) the difference δC(t) = C(t) − C S C (t) between the exact correlator and the semiclassical result. δC(t) obviously entails all those velocity correlations that are not included in the semiclassical description. While this quantity is generally unknown, as its calculation requires the full solution of the problem, we now provide an approximation scheme that has a very broad validity and that turns out to be quantitatively accurate for the problem at hand, where the dynamical disorder is slow compared to the free-carrier dynamics.
The key observation is that the quantum corrections δC induced by dynamical lattice fluctuations are continuously connected to those that would be realized in a perfectly frozen lattice environment, δC 0 : the latter should be adiabatically recovered when the timescale of lattice fluctuations is sufficiently long as compared to all other timescales in the problem. The similarity between δC and δC 0 for slowly fluctuating disorder is advantageous, because solving a problem with frozen disorder (which can be done exactly via the diagonalization of a one-body disordered Hamiltonian) is a much easier task than solving the full dynamical problem, which is instead prohibitively difficult. The quantity δC 0 = C 0 − C S C evaluated in the frozen disorder limit contains most of the information that we need on localization corrections. The small difference between δC 0 and δC, originating solely from the dynamical nature of the lattice, can then be treated in an approximate way. The power of the method relies on the fact that no assumptions are made on the smallness of δC itself, so that the whole approach is not restricted to the weak disorder limit.

B. Decorrelation time
Disorder dynamics (here, the dynamics of atomic and molecular positions) are known to destroy the quantum processes at the very origin of wavefunction localization 13,14 . In the scaling theory of localization, the motion of the scattering centers on the dynamical timescale τ d gives rise to a finite cutoff length 24 which corresponds to the length traveled by semiclassical particles over this time 13 ; this cutoff length/time restores a finite diffusivity for the carriers, which would otherwise be vanishing. When translated to our problem, this implies that the quantum corrections δC 0 characterizing the localized system can only be sustained at times that are short compared to the timescale of dynamic disorder 25,26 , while they decay and vanish at longer times (see Appendix C). This is embodied in the following form: with the decorrelation time being set by the frequency of the relevant modes, i.e. τ d ∼ 1/ω 0 within a numerical factor (the optimal value of such prefactor will be discussed below). The corresponding velocity correlator is then Eq.
(2) constitutes the basis of our theoretical approach, that will now be benchmarked and applied to the study of charge transport in organic materials in the next sections.
Before proceeding further, we argue that the proposed approximation scheme should be valid even beyond the slow disorder limit initially assumed in the derivation. The reason is that Eq. (2) is able to interpolate from full localization all the way to the semiclassical limit depending on the value of τ d . Indeed, while letting τ d → ∞ obviously restores the correlator of the localized system, C(t) = C 0 (t), taking the opposite limit of fast disorder, τ d → 0, suppresses all quantum corrections, hence leaving C(t) = C S C (t). More generally, in cases where localization corrections are irrelevant to start with (because lattice fluctuations are weak, i.e. for weak electron-phonon interactions and at low temperatures), then δC 0 can be set to zero. Correspondingly, the correct semiclassical result will be trivially recovered by Eq. (2) regardless of the value of τ d .

C. Carrier mobility
The carrier mobility can now be obtained from Eq. (2) following the lines of Ref. 15. We first observe that the diffusion constant D is the long-time limit of the instantaneous diffusivity, defined as Applying the Laplace transform to Eq. (2) yields with p = 1/τ d , and hence D = D S C + δC 0 (p)/2. The mobility is then obtained from Einstein's relation µ = eD/k B T as with µ band = eD S C k B T and δµ = e 2k B T δC 0 (p). The above expression has the desired form of a semiclassical band mobility corrected by quantum processes.
To make the connection with previous works, we observe that when the term δµ dominates in Eq. (5), i.e. when localization effects are sufficiently strong, one is allowed to neglect the semiclassical terms altogether, which corresponds to setting C = δC in the previous derivation. Eq. (2) then becomes C(t) = C 0 (t)e −t/τ d , which is the form that was originally assumed in Ref. 15, setting the basis of transient localization theory. Repeating the same steps as above, Eq. (5) reduces to the usual TL formula for the charge mobility, µ = epL 2 (p)/2k B T , with the transient localization length defined by L 2 (p) =C 0 (p)/p 15,16 .
For the practical implementation of Eq. (5), we actually rearrange Eq. (4) as which allows us to take full advantage of the numerical tools already developed in the context of transient localization theory. The first term in Eq. (6) is then easily recognized as the TL result, which is readily evaluated with the methods described in Refs. 27 and 17. The remaining terms between brackets now only involve semiclassical quantities, which are evaluated following the procedure described in Appendix A. Finally, the mobility Eq. (5), is obtained making use of the explicit formula Eq. (B12), which follows directly from Eq. (6). Full details are presented in Appendix A, B and D.

A. Models
The theory developed in the preceding sections is totally general and can be applied to a variety of electron-phonon interaction models. For illustrative purposes we shall mainly consider the following class of tight-binding Hamiltonians 17 .
which is broadly representative of the physics of organic semiconductors. Eq. (7) describes charge carriers moving on a molecular lattice, with nearest neighbor transfer integrals J δ in the different bond directions δ that are linearly modulated by the coupling to intermolecular modes x i,δ . Unless otherwise specified, the modes are assumed to be uncorrelated between different bonds, as described by the Hamiltonian is the typical frequency of intermolecular vibration. Correlated bond fluctuations, as well as Holstein-type local interactions with slow intramolecular modes, of the form H I = i α H c + i c i x i , can also be studied within the present general theoretical framework (see below and Appendix A 3). The interaction with high-frequency intra-molecular modes, which is not treated explicitly here, can be included via a rescaling of the transfer integrals J δ , corresponding to the usual polaronic band narrowing 28 .
Following Ref. 17 we consider a two-dimensional, hexagonal molecular lattice of unit spacing a, with nearest neighbors δ =a,b,c (see the sketch in Fig. 2). We take J = J 2 a + J 2 b + J 2 c as the energy unit, which fixes the scale of the band dispersion. For clarity of presentation, we shall focus on the common situation encountered in high-mobility molecular semiconductors, where two bond directions are equal by symmetry, so that the set of transfer integrals can be characterized by a single parameter θ (J a = J cos θ, 2) 17 . Moreover, we shall be mostly concerned with the high-temperature regime where k B T ω 0 . In this case the mean square thermal fluctuation of the transfer integrals is readily evaluated in terms of the parameters of Eq. (7) to ∆J δ = α δ x 2 i,δ 1/2 = α δ (k B T/K) 1/2 . We can then introduce ∆J = ∆J 2 a + ∆J 2 b + ∆J 2 c as a measure of the overall energetic disorder induced by the intermolecular displacements. Unless otherwise specified, we fix the microscopic parameters to J = 0.1eV and T/J = 0.25, representative of high-mobility organic semiconductors at room temperature 17 , and set = 1.

B. Validation of the theory
Before exploring our general findings we validate the theory by benchmarking it against exact results available on simplified one-dimensional models. We consider the quantum Monte Carlo (QMC) data of Ref. 29, and the Finite Temperature Lanczos Method (FTLM) data of Ref. 30. Both approaches are in principle exact, and they fully include the quantum dynamics of the molecular vibrations without any assumptions. The case studied in Ref. 29 corresponds to the one-dimensional limit of the model Eq.  Fig. 1(a) shows the QMC data for the mobility versus temperature (squares) together with the result obtained from the present theory (DLC, red full and dotted lines) and band theory, as described in Appendix A [Eq. (A32)] (gray, thin). An excellent quantitative agreement with the QMC data is obtained if one takes a value p = 1/τ d = 2.2ω 0 in the DLC theory (red, full line). To illustrate the impact of the decorrelation time on the mobility, we also show the result obtained for p = ω 0 (red dotted). Reducing the value of p tends to overestimate quantum localization effects, and hence to underestimate the mobility. We observe that with the present choice of model parameters, which corresponds to a moderate amount of molecular disorder (∆J/J = 0.41 at T/J = 0.25), the QMC result is itself qualitatively similar to the band theory result in the whole temperature range explored, and the reduction of the mobility by localization corrections is less than 15%. rubrene We also note that the QMC data do not recover exactly the calculated band value in the low temperature limit, where molecular fluctuations are suppressed. This could signal either the presence of scattering processes not included in the perturbative band result, or a numerical artefact brought by the analytical continuation in the QMC calculation. Finally, the TL result (orange dashed) is quantitatively accurate around room temperature, but it becomes inappropriate in the weak disorder regime attained at lower temperatures (see next Sections). Fig. 1(b) shows the exact FTLM data for the optical conductivity per particle (black, thin), together with the result obtained from the present theory (DLC, red full and dotted lines) and from band theory (gray, thin) (details on the calculation of the frequency-dependent response are provided in Appendix B). Because in this example the level of disorder is quite large (the thermal fluctuation of the on-site molecular energy is evaluated to ∆/J = 1), the improvement brought by the inclusion of dynamical localization corrections is more striking: not only the theory corrects the gross overestimate of the mobility (i.e. the value at ω → 0) implied by band theory, but it very accurately captures the whole frequency response, including the emergence of a localization peak at ω 2J and the precise shape of the absorption tail at higher frequencies. The discrepancy observed at the low frequency absorption edge could instead be related to the small size of the cluster studied in Ref. 30, which is limited to 6 sites. As in the previous case, also here the choice p = 2.2ω 0 provides the overall best agreement with the exact result. We note that at this large level of disorder, the DLC result is essentially indistinguishable from the TL result (not shown).

C. Breakdown of band transport
Having validated the theoretical approach, we now study the emergence of dynamical localization corrections in the broad class of organic semiconductors, by analyzing the ensemble of models described by Eq. (7). Fig. 2(a) illustrates the evolution of the quantum correction term δµ of Eq. (5) relative to the band value, as a function of the electronic structure   The quantum correction to the band mobility is predominantly negative and expectedly increases in magnitude upon increasing the amount of disorder. The importance of processes beyond Bloch-Boltzmann theory is very much dependent on the band structure, which is in agreement with the general observation that different crystal structures are differently affected by disorder, with isotropic two-dimensional bands being the most resilient to localization processes 17 . The correction term is indeed maximum for one-dimensional structures (θ = 0, π), while it is minimized at the isotropic point where it becomes practically negligible at ∆J/J = 0.3. The structures around π − θ 0 are exceptions, exhibiting a positive correction term: there, due to a negative combination of the signs of the transfer integrals, J a J b J c < 0, a van Hove singularity of states resilient to localization arises close to the hole band edge. Thermal population of these states causes the mobility to rise above the band value.
The detailed dependence of quantum corrections on the amount of disorder is illustrated in Fig. 2(b), where we show the ratio of the total calculated mobility µ to the band prediction µ band , for selected electronic structures (solid lines and points). For each structure, one can identify a value of ∆J above which the mobility significantly deviates from the band value, signaling the emergence of quantum processes. For example, band transport holds up to ∆J/J 0.3 in the isotropic structure (red), while it breaks down already at ∆J/J 0.1 in the one-dimensional case (dark gray). The neglect of quantum processes beyond this point can lead to gross quantitative errors in the estimated mobility: at the largest values of ∆J/J studied here, for example, µ can deviate from the band value by up to a factor of three.
In Fig. 2(b) we also report the transient localization result, shown as dashed lines. The latter is applicable in the strong disorder regime, where it closely matches the result of Eq. (5). Upon reducing the disorder strength, however, the TL result does not tend to the band limit as required, but instead it incorrectly bends downwards. Band theory and TL theory therefore appear to be complementary, each addressing a different regime of parameters. Eq. (5) allows to bridge continuously between these two limiting regimes.

D. Localization corrections in the time domain
The buildup of quantum localization processes can be directly visualized by tracking the time-dependent diffusivity D(t), as shown in Fig. 2(c) for ∆J/J = 0.3 and θ = 0.21 (these are the values of the microscopic parameters calculated for rubrene in Ref. 17). The semiclassical diffusivity (gray) exhibits a monotonic evolution from ballistic at short times, showing no hint of localization. Localization processes are instead fully developed when considering the exact evolution in a frozen molecular environment (black thin line), as derived from the reference correlation function C 0 introduced in Sec. II A. Their onset can be identified with the locus of the maximum of D(t), that we denote as t loc , and they are responsible for the subsequent steady decrease of the diffusivity, which vanishes at long times. When the disorder is dynamic, such localization corrections are initially retained, leading to a partial suppression of the diffusivity (red curve and hatched region) w.r.t. the band value. The suppression of D(t) however stops at times t τ d , which follows from the fact that its derivative δC(t) vanishes (cf. Eq. (1)).
By tracking the time evolution of the diffusivity, we can now better visualize what controls the emergence of dynamical quantum corrections in Eq. (5). When the disorder is sufficiently strong (Fig. 2(c)) the condition t loc < τ d is fulfilled, and localization corrections can develop before they are suppressed by decorrelation due to the molecular dynamics. Reducing the disorder strength makes localization processes less efficient, resulting in an increase of the localization time t loc . When the latter reaches τ d , the quantum corrections cannot develop anymore, because they are cut off at their very onset by decorrelation, cf. Eq. (1). Quantum processes then become irrelevant and band theory applies. Therefore, the condition t loc τ d is what marks the breakdown of semiclassical behavior. According to this argument, the emergence of dynamical localization corrections is subject to the existence of either sufficiently strong (low t loc ) or sufficiently slow (large τ d ) dis-order fluctuations. Both conditions are naturally realized in organic crystals, where large amplitude, slow intermolecular fluctuations arise owing to the large masses of the molecular constituents and to the weak intermolecular binding forces.
The disappearance of localization corrections at low disorder is shown in Fig. 2(d), where we report the diffusivity calculated for ∆J/J = 0.1, well within the band regime. In this case t loc > τ d and the full diffusivity (red) is essentially indistinguishable from the band result, implying δµ 0.
Finally, Fig. 2(d) also shows that applying TL theory when t loc > τ d incorrectly underestimates the diffusion constant (dashed line), which is at the origin of the non-monotonic behavior exhibited in Fig. 2(b).

E. Localization corrections in the frequency domain
The conduction properties under a constant applied field and the response of charge carriers in the frequency domain are deeply intertwined. In the band transport regime, the optical absorption exhibits a simple Lorenztian (Drude) shape, which is a monotonically decreasing function of frequency (cf. Appendix A). Upon increasing the disorder strength, the suppression of the mobility (and hence the d.c. conductivity) by localization processes induces a dip in the absorption at ω = 0, as already shown in Fig. 1(b). This shifts the absorption maximum to finite frequencies 15,29,[32][33][34] , providing a direct and measurable signature of the breakdown of band transport.
To assess how the emergence of dynamical localization corrections is reflected in the optical conductivity, we take advantage of the following exact formula 15 σ(ω) = 2 tanh βω 2 ∞ 0 dt sin(ωt)D(t), whose second derivative reads It is clear from the above relationship that whenever the diffusivity is a monotonically increasing function of time, so that D(t) is always lower than its long time limit D, the curvature is necessarily negative and the optical conductivity remains peaked at ω = 0. This is the case in particular for the band diffusivity shown in Figs. 2(c) and (d), in agreement with the resulting Drude behavior. A necessary condition for the emergence of a finite-frequency peak is instead the existence of a region where D(t) > D, as indeed happens when quantum corrections are relevant (red curve in Fig. 2(c)). This, together with the analysis of the time-dependent diffusivity given in the preceding paragraphs, shows that the emergence of a dip in the optical absorption essentially coincides with the crossover condition τ d ∼ t loc . The change of sign of Eq. 8 can therefore be used to identify the breakdown of the band picture. The equivalence between these two conditions is further explored next.

F. Transport phase diagram
The locus of the breakdown of band transport as determined from the condition d 2 σ/dω 2 | ω=0 = 0 is reported in Fig. 3 as a function of the electronic structure parameter θ (shaded area). The symbols locate the electronic structures and disorder levels calculated for a number of high-mobility organic semiconductors of current interest 17,35 . For completeness we have included materials that either exactly fulfill the condition J b = J c , or are sufficiently close to it. The shaded area delimits the crossover lines calculated for values of the molecular fluctuation time τ d comprised in the interval 10 − 20 /J, which is the typical range encountered in materials. Remarkably, all the reported compounds are characterized by sizable dynamical localization processes, making band transport theory inappropriate. Several of them, however, are located very close to the crossover region. In the case of rubrene, in particular, this is in agreement with the fact that a finite-frequency peak is observed in optical absorption experiments at room temperature 36,37 , yet a normal Drude shape is recovered upon reducing the amount of thermal disorder 38 . We stress that the reported crossover line corresponds to the ideal situation where the only source of randomness is from transfer integral fluctuations. In real conditions, local site-energy fluctuations originating from the coupling to intra-molecular vibrations as well as extrinsic sources of disorder will likely shift the crossover to lower values of ∆J/J.
At any rate, we observe that while low levels of disorder and isotropic band structures have been independently achieved in current materials, no compound exists yet that is able to combine such optimal features together. If such a compound could be synthesized, we argue that it would fully enter the band transport regime, possibly opening new perspectives for organic-based applications.

IV. CONCLUDING REMARKS
Poor conduction properties are commonly observed in broad and diverse classes of solids, including disordered and amorphous metals, polymers, materials with strong electronphonon interactions and strongly correlated electron systems. Regardless of the microscopic origin, all these "bad" conductors and semiconductors have in common a fundamental breakdown of the weak-scattering hypothesis underlying band transport theory. In this respect, we speculate that the mechanisms described in this work could generally apply to other classes of materials, where slowly fluctuating degrees of freedom of electronic, magnetic or vibrational origin are coupled to the carrier motion. While the microscopic analysis of specific cases is clearly beyond the scope of this work, we can take advantage of the physical insight gained here to derive a simple phenomenological formula that will be of practical use to assess the presence of dynamical localization corrections in general situations.
To this aim we go back to Sec. III D, where we have shown that localization processes arise when the localization time is shorter than the velocity decorrelation time, i.e t loc τ d .
We observe that this typically occurs in a regime where the level of disorder is sizable on the scale of the band energy (cf. Fig. 2). In this regime, the localization time t loc scales with the semiclassical scattering rate τ within a numerical factor, and the two merge in the strong disorder limit, where t loc τ. Specifically, we have checked from the data of Fig. 2 that even in the isotropic 2D case, where localization corrections are the weakest, the ratio t loc /τ falls below 2 as soon as ∆J/J 0.4. In 1D, this is attained already at ∆J/J 0.2. It is then appealing to rewrite Eq. (9) by replacing the (quantum) localization time τ loc with the more familiar (semiclassical) scattering time τ. We therefore write with η a numerical prefactor. As a consistency check for Eq. (10), when τ d → ∞ we recover the known result that in the low density limit appropriate to non-degenerate semiconductors, localization corrections are always relevant when disorder is static 13 . To fix the value of η, in Fig. 3 we compare Eq. (10) with the rigorous condition obtained from Eq. (8), for two different values of the fluctuation time, τ d = 10/J and 20/J. We find that a very good match is obtained for η 4 (shaded area and dotted lines respectively). The criterion Eq. (10) can now be straightforwardly applied to a variety physical situations of interest, provided that τ d is identified with the timescale of the relevant degrees of freedom that couple to the electron motion. It is interesting to compare this result with the phenomenological Mott-Ioffe-Regel (MIR) criterion, which signals the disappearance of Bloch states. In its original formulation 39 , this is identified as the point where the semiclassical meanfree-path becomes comparable to the lattice spacing, i.e. a. Using 2 = v 2 τ 2 , this can be rewritten in terms of the scattering rate as This condition is illustrated in Fig. 3 as a dashed line, and it is located well above the crossover region. The fact that deviations from band theory arise before the MIR limit is reached is consistent with the fact that the buildup of quantum localization corrections requires the existence of coherently propagating Bloch waves to start with. The comparison of Eqs. (10) and (11) shows that this can only happen in systems where the disorder dynamics is sufficiently slow, as is the case here. Finally, we note that a slightly different formulation of the MIR criterion is often used 40,41 , which predicts the disappearance of Bloch states when the semiclassical scattering rate 1/τ reaches a fraction of the bandwidth, i.e.
with ξ again a numerical factor. The two conditions Eq. (11) and Eq. (12) become equivalent in the ultra-high temperature limit T J, in which case v 2 1/2 /a ∝ J/ independent on temperature (the average velocity entering in Eq. (12) instead recovers instead v 2 ∝ k B T/m * at low temperature). Comparing Eq. (12) and Eq. (10), we again conclude that for sufficiently slow vibrations, i.e. J ω 0 , the dynamical localization corrections emerge in a region where Bloch states are not yet suppressed.

Appendix A: Band transport theory
Here we start from the textbook equations of Bloch-Boltzmann band transport theory to evaluate the timedependent diffusivity D S C (t) and anticommutator correlation function C S C (t) needed in the main text, and provide useful analytical formulas for the mobility in different models of electron-phonon coupling that are relevant to high-mobility organic semiconductors.

Time-dependent diffusivity
In band transport theory, one starts with Bloch states in the periodic molecular lattice, having momentum k, energy k and velocity v k in a given direction. The diffusion constant of each band state is D k = v 2 k τ k , with τ k a decay time determined by scattering off disorder and lattice vibrations (see Section A 3 below). The mobility is then obtained as µ band = eβD S C = βe v 2 k τ k , where the symbol · · · indicates thermal averaging over the states and β = 1/k B T , with k B the Boltzmann constant.
The time-dependent diffusivity is related to the real part of the frequency-dependent conductivity σ(ω) per particle via the following expression 15 which is exact for non-degenerate carriers. The optical conductivity per band carrier is given by 42 with Z = k e −β k the partition function. The corresponding diffusivity is obtained via direct integration of Eq. (A1). As Eq. (A2) is a linear superposition, the contributions from individual states can be separated as The integral can be performed via contour integration, leading to where we have introduced the bosonic Matsubara frequencies ω n = 2πn/β. The result of the last summation has a closed expression in terms of the Hurwitz-Lerch transcendent function Φ. Denoting z = e −2πt/β , we have n>0 2e −ωn t 1−(ω n τ k ) 2 = [zβ/π(τ k ) 2 ](Φ(z, 1, 1 + ω 1 τ k ) − Φ(z, 1, 1 − ω 1 τ k )). The anticommutator velocity-velocity correlation function C S C (t) can be straightforwardly obtained from Eq. (A4) using the definition C S C (t) = 2dD S C /dt.
In the high temperature/weak scattering limit, T 1/τ k , the explicit sum over Matsubara frequencies drops out and the above expression simplifies to D k (t) v 2 k τ k [1 − e −t/τ k ] (this result can be obtained straightforwardly from Eq. (A3) by taking the classical limit for the detailed balance factor, tanh(βω/2) → βω/2). This form corresponds to a simple exponential decay of the velocity correlation function 15 . It describes ballistic behavior at short times, D k ∼ v 2 k t, followed by diffusive behavior D k → const = v 2 k τ k at times t τ k . The diffusivity in this case is a monotonically increasing function of time. According to the full expression Eq. (A4), however, an anomalous time dependence can arise when the scattering rate becomes much larger than the thermal energy, in which case deviations from the simple monotonic form above can appear. This happens because the ballistic velocity at short times becomes larger than v k . Due to this initial overshoot, D k (t) reaches values larger than the long-time diffusion constant, which is then attained after going through a maximum. Such non-monotonic behavior arises when 1/τ k 3.15T . Note that, because the scattering rate for classical vibrations increases with √ T due to the equipartition principle, a regime of anomalous diffusivity can in principle be attained at low temperature.

Average scattering time
The full electrodynamic response of non-degenerate carriers in the Bloch-Boltzmann approximation, as given by Eq. (A2), is a superposition of Lorentzians of widths 1/τ k , weighted by the corresponding thermal factors. This, in principle, deviates from a simple Lorentzian shape, as would be predicted instead within Drude theory. A simpler approximation for σ(ω) 42 , and therefore for D S C (t), is obtained by introducing a single, k-independent relaxation time τ. The latter is univocally determined from the knowledge of the long-time diffusivity, as which is the proper thermal average of the transport scattering time τ k , as defined in Ref. 42 . By construction, the above equation recovers the correct diffusivity D S C = v 2 k τ = v 2 k τ k and mobility µ band = eβ v 2 k τ. The corresponding optical conductivity then takes the simple Drude form In all cases studied the time-dependent diffusivity and optical conductivity obtained from the average τ are either very close to or indistinguishable from those obtained from the full k-dependent expressions. We therefore use the former simplified framework for the evaluation of the quantum corrections in the main text.

Calculations on specific models
Let us consider the scattering of a k-state off phonon modes of momentum q and frequency ω 0 , as described by the interaction Hamiltonian Here N is the number of molecules, c + k , c k the creation and annihilation operators for carriers, x q,δ the deformation mode corresponding to a given bond direction δ, and we set = 1. Straightforward algebra allows to write the interaction matrix elements for uncorrelated bond disorder as [α (δ) k,q ] 2 = 4α 2 δ [cos((k + q/2) · δ)] 2 , with δ the vectors connecting nearestneighbours as shown in Fig. 1a and α δ = dJ δ /dx δ the sensitivity of the transfer integrals to intermolecular deformations. Upon substituting this expression, Eq. (A7) becomes equivalent to Eq. (1) of the main text, now expressed in momentum space. The canonical 2nd quantization expression for the electron-phonon interaction is obtained by expressing the bond coordinate in terms of dimensionless bosonic operators as x q,δ = (ω 0 /2K) 1/2 (b + q,δ + b q,δ ) with K the spring constant, so that the electron-phonon coupling matrix element becomes g (δ) k,q = (ω 0 /2K) 1/2 α (δ) k,q , and correspondingly g δ = (ω 0 /2K) 1/2 α δ . As is customarily done, we introduce a set of dimensionless coupling parameters λ δ = α 2 δ /(4J δ K) = (g 2 δ /ω 0 )/2J δ . The classical (thermal) fluctuation of the transfer integrals can be written as (∆J δ ) 2 = α 2 δ T/K = 2λ δ J δ T using the equipartition principle. In general, one can define global parameters J 2 = δ J 2 δ and ∆J 2 = δ ∆J 2 δ . In Sec. III C we consider a model where the relative fluctuations in the different bond directions are all equal, i.e. ∆J δ /J δ = ∆J/J for all δ. This corresponds to the choice of an isotropic coupling λ δ ≡ λ, and leads to (∆J) 2 = 2λJT in the thermal fluctuation regime.
Other models of interest can be put in the form of Eq. (A7). In the case of fully correlated bond disorder, as studied in Refs. 16,31 , the matrix element reads [α (δ) k,q ] 2 = 4α 2 S S H [sin((k + q) · δ) − sin(k · δ)] 2 , with now (∆J S S H ) 2 = 2α 2 S S H T/K = 4λ S S H JT . The prefactor 4 instead of 2 arises from the fact that the fluctuation of the transfer integral now arises from two independent modes located on adjacent sites. Finally, diagonal (intra-molecular) electron-phonon interactions correspond to a constant α k,q = α H , where α H measures the variation of the local molecular energy level with respect to an intra-molecular deformation x, and correspondingly ∆ 2 = α 2 H T/K = 2λ H JT . The momentum scattering rate is evaluated in d dimensions as with n b the phonon population and f k+q the Fermi occupation of the final state k + q, which can be set to zero at low carrier densities, and we have introduced the compact notation g 2 k,k+q = δ [g (δ) k,k+q ] 2 . The two terms between brackets originate from phonon emission and absorption respectively. The geometric factor F k,k+q = 1 − v k · v k+q /v 2 k measures the loss of momentum occurring at each scattering event, thereby differentiating the transport scattering time from the quasiparticle scattering time, which is instead obtained by setting F k,k+q = 1. The form of F k,k+q used here is the proper generalization to generic band structures of the textbook expression −k · q/k 2 , which only applies to isotropic, parabolic band dispersions. Note that the factor F k,k+q is often omitted in practical calculations 17,22 , generally leading to quantitatively incorrect values for the mobility (see below).
In the quasi-elastic limit where the intermolecular vibration frequencies set the smallest energy scale in the problem, ω 0 T, J the scattering time simplifies to From this result and from the expressions of g 2 k,k+q given above it is clear that in all models considered here, the band mobility at any given temperature in the classical regime T ω 0 scales with disorder strength as µ ∝ ∆J −2 , which reflects the secondorder nature of the scattering process.

Analytical results in 1D
For charge carriers in one space dimension, most calculations can be performed analytically. In what follows we express distances in units of the lattice spacing a, we set = 1 and mobility units to µ 0 = ea 2 / . Taking J a = J and J b = J c = 0, we have k = −2J cos(k) and v k = 2J sin(k) = 4J 2 − 2 k . The density of states is ρ(ν) = 1/ √ 4J 2 − ν 2 /π and we have The last equation allows to determine the average scattering time from the knowledge of the diffusion constant D, via Eq. (A5). The scattering time, diffusion constant and charge mobility can be calculated from Eq. (A9) in the following cases. a. Uncorrelated bond fluctuations ∆J, i.e. off-diagonal thermal disorder: where the last expression is the quasiparticle scattering time, obtained by setting F k,k+q = 1 in Eq. (A9). By performing the momentum integrations we obtain where I 0 and 0 F 1 denote respectively the modified Bessel function of the first kind and the regularized hypergeometric function. In the last equation, we have introduced (∆J/J) 2 = 2λT/J using the definition of the dimensionless electronphonon coupling λ = α 2 δ /(2K J). b. Diagonal disorder ∆. In this case we have µ = µ 0 2βJ cosh(2βJ) − sinh(2βJ) (π/2)(βJ) 2 (∆/J) 2 I 0 (2βJ) (A20) where we have used (∆/J) 2 = 2λ H T/J. Note that the transport scattering rate is formally equivalent to that arising from uncorrelated bond disorder (this equivalence is instead lost for the quasiparticle scattering time). As a result, the mobility has exactly the same functional form as in Eqs. (A14-A16), but with a reduced prefactor: for a given energetic spread ∆, diagonal disorder appears to be 4 times less effective than offdiagonal disorder ∆J in limiting the charge mobility, which is therefore 4 times larger. This result, which is exact in one dimension, remains approximately true in all the twodimensional structures studied here. c. Correlated bond fluctuations ∆J. In this case we have with λ = α 2 δ /(2K J), so that (∆J/J) 2 = 4λT/J (see above). d. Summary. Fig. 4 illustrates the temperature dependence of the mobility obtained, in the quasi-static limit, in the different cases studied above. We first observe that the effect of the geometrical factor F k,k+q is rather small for uncorrelated bond disorder, especially at low temperature, and it is strictly irrelevant in the case of diagonal disorder, cf. Eq. (A18). This is in contrast with the case of correlated bond disorder, where its neglect leads to an overestimate of the mobility by a factor of 2 [cf. Eq. (A23)], as already reported in Ref. 31 . Second, the temperature dependence for uncorrelated bond disorder in one dimension (and, similarly, for diagonal disorder) is rather weak in the relevant temperature range around and below room temperature. At T/J = 0.25 the temperature exponent of the mobility µ ∝ T −γ is γ = 0.66, and it tends to γ = 0.5 when T J, cf. Eq. (A16). The exponent is slightly reduced, γ = 0.46, if the geometrical vertex factor is neglected, as reported in Fig. 3a of Ref. 17 , again recovering γ = 0.5 when T J. The temperature exponent is larger for correlated bond disorder, γ 1.5, cf. Eq. (A26), which holds for all T J. We note that the mobilities calculated for correlated and uncorrelated bond disorder have accidentally very similar values at room temperature.

Quantum phonons
In all cases reported above we have considered the quasielastic limit ω 0 k B T , by setting explicitly ω 0 → 0 in Eq. (A9). We illustrate here how the quantum nature of the phonons is restored in the specific case of correlated disorder, by using the full expression Eq. (A8). The calculations can be straightforwardly extended to the other cases considered above.
For each of the two processes described in Eq. (A8), i.e. phonon emission and absorption, the delta function yields two solutions The corresponding squared matrix element g 2 k,k+q is equal to 4g 2 times 1 + sin 2 k − k ± ω 0 2t 2 − 2 sin k 1 − k ± ω 0 2t 2 (A29) at q = q + , and 1 + sin 2 k − k ± ω 0 2t 2 + 2 sin k 1 − k ± ω 0 2t 2 (A30) at q = q − . In the absence of vertex corrections, i.e. setting F k,k+q = 1, the square root term cancels from the sum, leading to: Including the geometrical vertex F k,k+q = 1 − v k · v k+q /v 2 k restores the square root terms in Eqs. (A29) and (A30). After some elementary algebra the transport scattering time is then obtained as The calculation of the diffusion constant D and of the mobility µ now proceed as in the classical case. The energy integrals v 2 k τ k cannot be cast in closed analytical form and must be performed numerically.
We report for completeness the expressions obtained for uncorrelated disorder, i.e.
and for diagonal disorder: . The inclusion of quantum phonons has two effects: the first is that the amount of intermolecular fluctuations increases as compared to the classical value, due to the fact that 1 + 2n b > 2T/ω 0 , with the classical limit only being attained asymptotically for T ω 0 . In particular, intermolecular fluctuations do not vanish at zero temperature as predicted by the classical formula but rather saturate to (∆J/J) 2 = λ δ ω 0 /J. This effect tends to increase the scattering rate, and hence to decrease the mobility with respect to the classical value. The second effect is that, due to energy conservation, scattering by phonons is suppressed within a shell of ±ω 0 around k , which leads instead to an increase in mobility. As shown in Fig. 5(a), the latter always dominates and the mobility with quantum phonons is larger than what is predicted in the quasistatic limit. The correction is around 2% at T = 0.25J = 5ω 0 , validating the use of the quasi-static expression when investigating room temperature mobilities. The discrepancy between the quantum and quasi-static results however grows rapidly upon reducing the temperature, (it is about a factor of 2 already at T = ω 0 ) so that one should use the full quantum expression when addressing the low-temperature properties. The apparent power-law exponent γ also rapidly increases upon reducing the temperature. The bottom panel in Fig. 5 shows an analogous comparison in the case of correlated bond disorder [Eq. (A22) (orange), Eq. (A32) (blue)], showing a qualitatively similar behavior. The correction in this case is respectively 7% at T = 5ω 0 , and 60% at T = ω 0 . The results reported in Fig. 5 also indicate that the effect of the geometrical vertex can be slightly modified by the inclusion of phonon quantum fluctuations. This is especially true in the case of correlated disorder, where the exact reduction of a factor 2 in the mobility in the quasi-static phonon limit [cf. Eqs. (A22) and (A23)], is increased in the quantum phonon case as soon as T ω 0 , reaching a factor 3 when T ω 0 . For uncorrelated disorder, vertex corrections become negligible at low temperature both in the quasi-static and in the quantum phonon case.

Errata on previously reported results
Eq. (A25) for correlated bond disorder correctly appears in Refs. 31 and 16 . In both cases, however, there is a misprint in the low-temperature expansion, which appears with the numerical prefactor 1/8 instead of 1/16, which actually corresponds to the result calculated in the absence of the geometrical vertex F k,k+q (this result was used in the figures of these papers in order to provide a direct comparison with the results of the Kubo bubble apporximation, which does not contain vertex corrections). An analogous overestimate of the mobility is reported in Fig. 1(d) of Ref. 29 . The correct band theory result is the one reported in Fig. 1(a) above.
Similarly, the Boltzmann result reported as a dashed line in Fig. 3b of Ref. 17 for uncorrelated bond disorder was also overestimated by a factor of 2. The correct result therefore moves closer to the transient localization value (see also Fig. 1b of the present work), although the significant discrepancies pointed out in that work remain. The temperature exponent shown as a dashed line in Fig. 3a  To evaluate the optical conductivity we consider the Kubo formula where C − (ω) is the current-commutator correlation function (we use units in which e = = 1 and a unit lattice spacing). The detailed balance condition relates C − (ω) to the anticommutator correlation function C(ω) Combining Eqs. (B1,B2) yields From Eq.
(2) of the main text we obtain where F T indicates the Fourier transform and p = τ −1 d is the inverse of a velocity decorrelation time (see text and Appendix C below). Substituting this expression in the exact relations Eqs. (B1,B2) yields, for the real part of the optical conductivity: and the semiclassical term in Eq. (B5) is given by Eq. (A6), i.e.
where σ S C (0) and τ are respectively the DC conductivity and the relaxation time obtained within the semiclassical approximation 42 and evaluated in Section A.
The term in Eq. (B7) can be calculated via the the Lehman representation (Ref. 23 , Eq. (A12)) as Direct inspection of Eq. (B9) shows that lim p→0 σ S C (ω, p) = σ S C (ω), so that Eq. (B5) recovers the optical conductivity of the statically disordered Hamiltonian when p → 0. In the opposite limit, p → ∞, the contributions from Eqs. (B6,B7) vanish as can be verified explicitly from the form Eq. (B9). It is worth noting that at any finite p, σ S C (ω, p) is not a simple Lorentzian convolution of σ S C (ω).

Mobility
The quantum corrections to the mobility are explicitly derived in the main text. Here we show how these can be alternatively obtained as the ω → 0 limit of the AC conductivity. Using the relation 23 which is Eq. (5) in the main text. For practical calculations, we recognize that the second term in Eq. (B11) is the mobility µ T L obtained from transient localization theory, which can be calculated by standard methods 17,27 . To evaluate the remaining terms we perform explicitly the integral in Eq. (B9) and take the limit ω → 0, which yields the final expression denotes the integral representation of the z th harmonic number. A similar expression can be obtained for the full frequency-dependent conductivity.

Crossover condition
According to Eq. (B5), the second derivative of the optical conductivity at ω = 0 consists of three parts. The semiclassical terms are directly evaluated from the explicit expression Eq. (B8) The transient localization term is calculated from the exact eigenvectors (|n ) and eigenvalues (E n ) of the statically disordered Hamiltonian as where we have defined ω nm = E n − E m . The crossover from band transport to transient localization is found when the sum of the three terms appearing in Eqs. (B14,B15,B17) vanishes. Fig. 6 illustrates the crossover line obtained for p = 0.05 (τ d = 20 /J) using four different system sizes, showing that finitesize effects are well controlled. The results reported in Fig.  3 of the main text correspond to the largest clusters studied, consisting of 32x32 sites. Similar to Eq. (B13), a compact analytical expression for the term appearing in Eq. (B15) can be obtained in the large temperature limit βp 1, which corresponds to replacing Eq (B5) with the simpler form σ(ω) = σ S C (ω) + σ 0 (ω + ip) − σ S C (ω + ip) with σ 0 (ω) the optical conductivity of the the random static Hamiltonian. Taking the second derivative yields d 2 σ(ω) dω 2 | ω=0 = ∂ 2 σ 0 (ω, p) ∂ω 2 | ω=0 + 2σ S C (0)τ 2 1 − (1 + pτ) 3 (1 + pτ) 3 . (B19) This approximate expression gives results comparable to the original Eq. (B5) in all the cases studied in the main text.
Appendix C: Thouless argument and identification of the parameter τ d The effects of decorrelation due to dynamic degrees of freedom in disordered systems have been studied within the frame of the scaling theory of localization in the past 13 . The typical situation addressed in these works corresponds to disordered systems in which a (small) source of inelastic scattering is added, e.g. due to interaction of charge carriers with phonons. Here we show how the Thouless argument 24 , that was devised to deal with this situation, can be generalized to the case of purely dynamic disorder (i.e. where the lattice vibrations act both as the source of disorder and of decorrelation) provided that a correct identification of the decorrelation time is made.
According to the Thouless argument, in a disordered system in the presence of inelastic scattering of phononic or electronic origin, there exists a finite length scale L T h beyond which diffusive motion can take place. The Thouless length is related to the inelastic scattering time τ in via L 2 T h = D S C τ in , i.e. it can be interpreted as the diffusion length on the inelastic scattering timescale, which acts as a decorrelation time. According to the scaling theory of localization in two dimensions, this corresponds to a finite DC conductivity 13 σ 2D (L th ) = σ S C + δσ DC (L th ) (C1) where the two-dimensional correction to the conductivity due to weak localization is δσ DC (L th ) = e 2 2π 2 log with a microscopic length enforcing a lower cutoff in the scaling theory. The presence of a finite length (L 2 T h ) in the scaling term Eq. (C1) implies that the DC conductivity is nonvanishing, at variance with the case of a purely static disorder. We can interpret this result in terms of an infrared cutoff ω * in the quantum correction of the AC conductivity, by imposing δσ DC (L th ) = δσ AC (ω * ). (C3) To determine the corresponding ω * , we make use of the scaling theory expression for the AC conductivity 13 δσ AC (ω) = σ S C k F log |ωτ|.
Equating the above result with Eq. (C2) yields If τ in τ, the microscopic length can be identified with the semiclassical diffusion length, 2 D S C τ, and consequently ω * = 1/τ in .
Eq. (B18) allows us to directly compare our results with that coming from the previous argument. Let us now consider Eq. (B18). In the ω = 0 limit, the quantum corrections to the conductivity are δσ = Reσ 0 (ip) − Reσ S C (ip). (C6) Continuing Eq. (C4) to the complex plane and substituting into Eq. (C6) yields Comparing Eq. (C7) in the pτ 1 limit with Eq. (C3) we get p = ω * = 1/τ in . Since p = 1/τ d , we can therefore identify the inelastic time in the Thouless argument with the decorrelation time introduced in the main text. We note that in the standard situation where preexisting localization effects are decorrelated by the inclusion of phonons, the decorrelation time is expected to diverge when the electron-phonon coupling strength vanishes. In the case considered in this work, instead, where the disorder is itself of dynamical origin, the decorrelation time diverges when the disorder becomes static, i.e when the phonon frequency ω 0 vanishes. It is then natural to assume 25,26 p ∝ ω 0 . Beyond the analytical argument given here, the comparison with the exact QMC results shown in Fig. 1(a) for the correlated bond disorder model, and with the FTLM results of Fig. 1(b) for the Holstein model confirms the proprotionality between p and ω 0 , and indicates that the choice p 2.2ω 0 provides the most accurate quantitative results. This identification is also compatible with the values in the range p = 2 − 2.9ω 0 inferred from best fits of the Ehrenfest dynamics of Ref. 15 , performed at various values of λ and ω 0 in the strong disorder regime.
• for all other cases, i.e. for the calculation of the quantum corrections to the mobility and for the determination of the crossover from band transport to transient localization, we use the exact diagonalization method described in Ref. 27 . The maximum system size studied here, of 32x32 sites, is sufficient to reach convergence on the quantities of interest (see Fig. 6).