A unification of the Holstein polaron and dynamic disorder pictures of charge transport in organic semiconductors

We present a unified and nonperturbative method for calculating spectral and transport properties of Hamiltonians with simultaneous Holstein (diagonal) and Peierls (off-diagonal) electron-phonon coupling. Our approach is motivated by the separation of energy scales in semiconducting organic molecular cystals, in which electrons couple to high-frequency intramolecular Holstein modes and to low-frequency intermolecular Peierls modes. We treat Peierls modes as quasi-classical dynamic disorder, while Holstein modes are included with a Lang-Firsov polaron transformation and no narrow-band approximation. Our method reduces to the popular polaron picture due to Holstein coupling and the dynamic disorder picture due to Peierls coupling. We derive an expression for efficient numerical evaluation of the frequency-resolved optical conductivity based on the Kubo formula and obtain the DC mobility from its zero-frequency component. We also use our method to calculate the electron-addition Green's function corresponding to the inverse photoemission spectrum. For realistic parameters, temperature-dependent DC mobility is largely determined by the Peierls-induced dynamic disorder with minor quantitative corrections due to polaronic band-narrowing, and an activated regime is not observed at relevant temperatures. In contrast, for frequency-resolved observables, a quantum mechanical treatment of the Holstein coupling is qualitatively important for capturing the phonon replica satellite structure.


I. INTRODUCTION
Organic semiconductors are a promising class of soft materials with applications in photovoltaics, display technologies and plastic electronics. 1-3 In particular, ultrapure organic molecular crystals (OMCs) provide a rich test-bed for the fundamental mechanisms and limitations of charge transport in soft materials. 4,5 Exceptional room-temperature mobilities in the tens of cm 2 V −1 s −1 have been measured in acene single crystals such as rubrene and pentacene. [6][7][8] These large mobilities are generally accompanied by an apparently "band-like" power-law temperature dependence, µ ∝ T −γ (γ > 0), suggesting significant carrier delocalization over many unit cells; however controversy still exists over the microscopic mechanisms leading to this behavior, with both experimental and theoretical evidence suggesting an interplay between coherent dynamics, localization, and incoherent hopping. 9,10 Theoretical efforts to understand electronic dynamics in OMCs has led to the development of several competing and complementary pictures that capture different limiting behaviors. For example, in the limit of weak electronphonon coupling, Boltzmann transport theory is quite successful, [11][12][13][14][15] similar to its use in simple inorganic semiconductors. However, many of the most interesting materials exhibit strong electron-phonon coupling, producing localized electronic states and invalidating the use of Boltzmann transport theory. 13 In actuality, many of the material parameters that govern charge transport in most OMCs -the electronic bandwidth, phonon energy, and electron-phonon coupling strength -exist on energy scales that are comparable to one another and to k B T , precluding any simple perturbative treatment.
In this state of affairs, two theoretical pictures have emerged, which focus on two distinct forms of electronphonon coupling. The first picture is that of the Holstein polaron formed by strong local coupling to intramolecular vibrations . 16 In this picture, electronic motion is generally thought of as incoherent, and nuclear tunnelling quantum effects are non-negligible. [17][18][19] The second picture is that of dynamic disorder or transient localization due to strong nonlocal coupling to intermolecular vibrations. [20][21][22][23] In this case, the delocalized coherent evolution of the electronic wavefunction cannot be neglected, but nuclear quantum effects are considered insignificant. In many real OMCs, both forms of electronphonon coupling are present. However, the different physical ingredients underlying these two theories prevents their straightforward unification. In this work, we achieve this goal and present a unified nonperturbative approach to calculating dynamical observables in OMCs. Importantly, our method reduces to the two powerful theories described above in their respective limits of validity, and gives insight into the relative importance of each energy scale and the crossover from one regime to another. Furthermore, the method is generalizable and computationally affordable, suggesting its straightforward future application to real materials beyond model Hamiltonians.
The layout of this manuscript is as follows. In Sec. II, we describe the Hamiltonian and theoretical and computational tools used to simulate dynamical properties of interest, in particular the conductivity and the spectral function. We further discuss the relation to previous works and the validity of our approximations, referring to an Appendix for a comparison to numerically exact results that we generate. In Sec. III, we present the results of our simulations as a function of temperature, analyzing spectral properties and the temperature dependence of the mobility. Our results are presented at various values of the strength of both types of electron-phonon coupling. Enabled by our unified theory, we identify a mechanism by which an increase in the electron-phonon coupling strength counterintuitively yields a more delocalized wavefunction. In Sec. IV, we summarize our results and discuss future directions.

A. Hamiltonian
We study the spinless Holstein-Peierls model with the Hamiltonian H = H el + H ph + H el−ph with where a † i (a i ) and b † i (b i ) are the creation (annihilation) operators on lattice site i for electrons and phonons, and i j indicates nearest neighbors. Here and throughout, we use = 1 for notational simplicity, but all results will be given with proper physical units. In Eq. (1a), h i j are one-electron Hamiltonian matrix elements. Due to particle-hole symmetry our analysis applies equally to electron and hole doping, and in this work we use parameters inspired by hole transport in rubrene. [24][25][26][27] Although our ensuing formalism can be applied to an arbitrary electronic Hamiltonian, here we will study a nearest-neighbor tight-binding model with h i j = −τ for i and j nearest neighbors and zero otherwise. This model represents one molecular orbital per lattice site and we consider two vibrational modes on each site with frequencies ω H and ω P . The first is a high-frequency intramolecular vibration whose dimensionless displacement X i,H ≡ (b † i,H + b i,H ) couples to the on-site electron density via the Holstein mechanism with dimensionless strength g H . The second is a low-frequency lattice mode corresponding to the position of the molecule in the unit cell; the difference in displacements of neighboring molecules, (X i,P − X j,P ), couples to the kinetic energy via the the Peierls (or Su-Schrieffer-Heeger 28 ) mechanism with dimensionless strength g P . These two electron-phonon coupling mechanisms are commonly referred to as being local and nonlocal, respectively.
The Hamiltonian parameters relevant for molecular organic crystals presents a challenge to theoretical analysis. For example, the parameters we will use in this work are τ = 100 meV, ω H = 150 meV, and ω P = 6 meV. These exemplify the typical behavior ω H /τ > 1 and ω P /τ < 1, i.e. the intramolecular dynamics occur on timescale that is faster than the electronic dynamics and the intermolecular vibrational dynamics occur on a timescale that is slower than the electronic dynamics. Furthermore, the dimensionless electron-phonon coupling strengths g H and g P are both on the order of 1 and the temperature ranges from 1 to 40 meV. Collectively, these parameters preclude an obvious perturbative treatment of the combined electron-phonon interactions. In the remainder of this section, we will discuss our nonperturbative approach to simulating the dynamics of this Hamiltonian with the knowledge that the Holstein and Peierls components describe fun-damentally different physical processes.

B. Conductivity and mobility
Our principle object of interest is the temperaturedependent AC conductivity σ(ω), which is obtained from the current autocorrelation function, where Z = Tre −βH and L is the number of lattice sites. Although in general, the conductivity is a tensor quantity with respect to the lattice vectors, for the sake of simplicity we will focus on the one-dimensional case. The DC conductivity, σ DC ≡ σ(ω → 0), gives access to the mobility via µ = σ DC /ne 0 , where n = N/L is the number density of carriers and e 0 is the fundamental electric charge. The current operator is (3) and a is the lattice constant and J i j is an operator in the Peierls phonon subspace. We work on a 1D lattice with periodic boundary conditions in the canonical ensemble with N = 1 electron, which is appropriate for the low density limit achievable in most semiconductors through doping or photoexcitation.
We treat the Peierls phonons quasi-classically, which is justified based on their relatively low frequency. Although approximate, this treatment is nonperturbative in g P . Under this approximation, the Peierls phonon displacements are scalar variables that obey classical equations of motion X i,P (t) = X i,P (0) cos(ω P t) + P i,P (0) sin(ω P t) with local Hamiltonian H i,P = (ω P /4) P 2 i,P + X 2 i,P . This leads to a Holstein Hamiltonian with time-dependent electronic matrix elements, and a time-dependent current operator where The current autocorrelation function is now approximately given by C JJ (t) = dX P dP P P(X P , P P ) where P(X P , P P ) is the phase-space distribution of the classical Peierls variables and U H (t, 0) is the time-ordered evolution operator, In order to approximately respect quantum statistics at low temperature, especially zero-point motion, the phase-space distribution P is a product of Gaussian Wigner distributions with X 2 i,P = tanh(βω P /2) and P 2 i,P = 1/[tanh(βω P /2)]. The integral is then evaluated by Monte Carlo sampling of X P and P P , and these variables are then evolve in time according to classical dynamics.
To treat the Holstein electron-phonon coupling, we use the Lang-Firsov polaron transformation, 29 defined for operators O asÕ ≡ e S Oe −S with Applying the transformation to H H (t) and J(t) gives where H ω H is the polaron self-energy. We now realize a separable, finite-temperature mean-field Hamiltonian by replacing the phonon operator e C i e −C j by its thermal average, This approximation assumes that the Holstein phonons equilibrate much faster than the electrons, which is reasonable because of their high frequency ω H > τ. We do not make any approximation in the polaronic treatment ofJ.
Combining the above procedures, the total Hamiltonian has now been approximated as a sum of separable and quadratic electron and phonon terms, where the electronic Hamiltoniañ h(t) is time-and temperature-dependent, with the evolution operator The transfer integral is given bỹ and reflects Peierls-induced dynamic disorder and Holsteininduced band narrowing. Finally, the correlation function becomes where J i jkl (t) and F i jkl (t) are purely electronic and phononic correlation functions, × Tr el ũ(0, t)a † i a jũ (t, 0)a † k a l e −βh(0) /Z el (18) and where The time dependence of F i jkl (t) is responsible for incoherent electron-phonon scattering and suggests a separation into static coherent and dynamical incoherent contributions, 30,31 Retaining only F coh produces a theory of dynamically disordered transport with Holstein-induced band narrowing of the matrix elements in both the Hamiltonian and the current operator.
We have implemented an algorithm for Eq. (18) based on the time evolution of single-particle eigenstates which, for a general electronic Hamiltonian, scales as L 5 N t per trajectory where L is the number of lattice sites and N t is the number of timesteps. If the electronic Hamiltonian has only finiterange interactions, such as the nearest-neighbor interaction used here, then the scaling is reduced to L 3 N t . Finally, if the incoherent phonon term is neglected, the nearest-neighbor case has a scaling of L 2 N t , which is the typical cost of the nearest-neighbor semiclassical dynamic disorder approach at finite temperature. C. One-particle spectral function The above approach can be straightforwardly extended to many other correlation functions. In particular, we will also present results for the inverse photoemission spectrum of the undoped parent semiconductor corresponding to a measurement of the conduction band's many-body density of states. 11 This is calculated via the electron addition Green's function of the model in the ground state with N = 0 electrons, using the same approximations introduced above. From this, we calculate the momentum-resolved spectral function and the density of states ρ(ω) = L −1 k A(k, ω).

D. Limiting behaviors and relation to previous works
In the limit g P = 0, i.e. without Peierls electron-phonon coupling, our approach reduces to that presented by Ortmann et al. 30 In this limit, the Hamiltonian has no time dependence and the dynamics can be re-written exactly in terms of eigenstates ofh. At low temperature, this approach reduces to Boltzmann transport theory with a constant scattering time t s and band-narrowed transfer integralτ. 32 The theory also reproduces the narrow-band and Marcus-Levich mobility in the limit of large polaron binding energy (large g H or ω H ), and ultimately Marcus theory in the limit of large polaron binding energy and high T . 33 A major distinguishing feature of our approach, in comparison to previous studies that have combined Holstein polaron approach with dynamic disorder, 19 is that we make no narrow-band approximation, which means that we recover the exact unitary electronic dynamics in the absence of Holstein coupling.
In the limit that g H = 0, i.e. without Holstein electronphonon coupling, this approach is analogous to the "dynamic disorder" picture of Troisi and co-workers 20 or the transient localization scenario of Ciuchi, Fratini and Mayou. 22,23 In previous work, 20,21,24,34 the evolution operator Eq. (9) is used to propagate a localized electronic state and the mean squared displacement (MSD) is calculated; in the long-time limit, the slope of the MSD is used to determine the diffusion coefficient and the electron mobility. This approach commonly includes feedback of the electronic system on the classical phonons at a mean-field level (Ehrenfest dynamics). We neglect this force, resulting in the analytical solution for the phonon dynamics in Eq. (4); this has been shown to be a safe approximation in the medium-to high-mobility regime studied here. 21 Ignoring this back-reaction, our Kubo approach has three advantages. First, quasiclassical dynamics deteriorates in the longtime limit, approaching a behavior that is similar to that at infinite temperature. 23,35 Compared to the MSD, the current autocorrelation function is less sensitive to the long-time dynamics because it decays to zero. Second, our approach naturally gives access to the full frequency-resolved conductivity rather than just the DC conductivity (the mobility). Third, the incorporation of a fully quantum Holstein electron-phonon interaction and its treatment by a similarity transformation is natural for the correlation function but much harder for the wavefunction dynamics.

E. Discussion of approximations
The accuracy of our approach ultimately relies on a separation of timescales between the two types of phonons and the electronic system. A semiclassical treatment of lowfrequency Peierls modes in molecular crystals has been shown to be well-justified and highly accurate. 23,36 The ratio of energy scales, ω P /τ = 0.06, indicates a strongly adiabatic regime, while the temperature k B T ω P in typical experimental conditions further justifies the use of classical phonons.
With respect to the Holstein modes, while the criterion of ω H /τ > 1 is satisfied, this ratio (1.5 in our case) indicates that the system is only weakly anti-adiabatic, which would limit the accuracy of our approach; many OMCs will have a smaller bandwidth than considered here, and thus ω H /τ 1 may be more strongly satisfied. To the best of our knowledge, no systematic and detailed study has been performed that compares the Lang-Firsov approach to numerically exact results for the Holstein model in this regime. Therefore, in the Appendix, we give a detailed analysis of the accuracy of the Lang-Firsov treatment compared to numerically exact calculations obtained via exact diagonalization over the variational Hilbert space (EDVHS). 37,38 The comparison is performed at zero temperature, which is useful to assess accuracy even at finite temperatures, since ω H k B T for the temperature range of interest. We find that despite application to a challenging parameter regime, the approximate Lang-Firsov approach captures the essential physics of the Holstein polaron in the form of coherent band-narrowing and the incoherent electron-phonon interaction. The accuracy of the Lang-Firsov approach holds in the limits of both weak and strong electron-phonon coupling, demonstrating its nonperturbative nature. Crucially, the physics enabled by our approach go well beyond those with perturbation theory and the common narrow-band approximation.
The idea of combining classical slow phonons with quantum high-frequency modes has been successful in the past in studies of nonadiabatic dynamics [39][40][41][42][43][44][45] and to a limited extent in OMC transport studies; however, such an approach generally involves additional approximations such as static lowfrequency phonons 46 or a zero-bandwidth perturbative treatment of high-frequency phonons. 19

A. Simulation Details
Numerical simulations were performed for a 1D periodic lattice with L = 64 sites for the full theory and L = 128 sites for the pure Holstein polaron and Peierls theories. Phonons frequencies are ω P = 6 meV and ω H = 150 meV and the electronic nearest-neighbor transfer integral is τ = 100 meV. These parameters were chosen to simulate the high-mobility axis of an anisotropic organic single crystal such as rubrene or pentacene, particularly when using a coupling strength of g H = 1 and g P = 1 − 2. [24][25][26][27] In keeping with this, we also use the lattice constant for the b-axis of rubrene, 7.2 Å. In addition to being representative of OMCs, these parameters satisfy the conditions necessary to ensure the accuracy of the method presented in this work, ω H > τ > ω P , as discussed above. Numerical propagation was performed using the fourth-order Runge-Kutta algorithm with a timestep set to be a factor of ten smaller than the shortest timescale in the problem, dt = 0.1 * min( /(4τ), 1/ω H ), which in practice varies between 0.05 and 0.4 fs. For DC transport properties, to ensure convergence at zero frequency in a finite system, we apply a Gaussian damping to the correlation functions, simulating a weak dispersionless Drude scatterer such as neutral impurity. 32 We use a scattering time of /t s = 0.25 meV, which is smaller than any energy scale of the problem by more than an order of magnitude. The choice of t s affects the absolute magnitude of the mobility in cases where the current autocorrelation does not decay to zero naturally due to other coupling, but our choice has no effect on the overall temperature dependence. While the spectrally-resolved quantities are obtained from the same calculation, we instead apply an exponential damping in time corresponding to a Lorentzian broadening in energy with half-width at half maximum η = 2.5 meV to produce smooth spectra. Convergence with respect to the ensemble sampling of the classical Peierls modes required averaging of up to 10, 000 trajectories.

B. Spectral functions
Before exploring transport, it is informative to look at the one-particle spectral properties of the Holstein-Peierls Hamiltonian within our approach, as described in Sec. II C. In Fig. 1, we show the momentum-resolved spectral function A(k, ω) at 300 K, with and without Peierls dynamic disorder and with various Holstein coupling values. In the absence of any electron-phonon coupling (top left), the spectral function has sharp quasiparticle peaks following the expected noninteracting dispersion relation E(k) = −2τ cos(ka). When moderate Holstein coupling is added (top center), we see a mixture of coherent and incoherent features in the spectral function. The zero-phonon transition is shifted by the selfenergy Σ H = −g 2 H ω H and has a dispersion whose bandwidth is narrowed to 4τ. Furthermore, satellite peaks now appear with a spacing of ω H , as seen in previous results. [47][48][49] With strong Holstein coupling (top right), the renormalized electronic bandwidth is very small, leading to a dispersionless vibronic progression with a spacing of ω H and intensities determined by g H and T .
In the case with Peierls coupling only (lower left), we see a broadened version of the purely electronic band structure, but with an anomalous peak at the band center. This unphysical peak arises in models with off-diagonal disorder. 50,51 Although the peaks acquire a linewidth, the quasiparticle picture is maintained. Dynamic disorder also yields states outside of the undisordered band and the onset of the band edge is softened. We recall that in the static (Anderson localization) limit, all states are localized in one dimension; however the localization lengths of states near the band edge are much smaller than those in the band center. 13 We will return to this point in Sec. III D when we analyze our mobility results in terms of localization properties.
For the case of simultaneous Peierls and moderate Holstein coupling (bottom center), we see that the dispersion and satellite structure largely depend on the Holstein coupling g H , whereas the Peierls dynamic disorder contributes a broadening in frequency and momentum (as well as the spurious peak at the band center, which is now present in the phonon replica structures). These effects of dynamic disorder become less pronounced as the band narrows, and in the case with renormalized bandwidth approaching zero (bottom right), the Peierls disorder has no discernible effect.

C. Optical conductivity
We now turn our attention to the primary focus of this work, the conductivity. In Fig. 2 we present the AC conductivity at different temperatures for the pure (disorder-free) Holstein model (g P = 0, left column) and for the full Holstein-Peierls model with strong dynamic disorder (g P = 2, right column). For the pure Holstein model, due to the aforementioned condition ω H k B T , we see little qualitative difference between the T = 100 K and T = 300 K cases, Fig. 2(a) and (c), and indeed little difference from the zero-temperature result in Fig. 7(a). In order to see qualitative changes in the pure Holstein conductivity, we need to go to much higher temperatures (T = 1000 K), at which point the phonon peaks associated with the incoherent terms in Eq. (21b) become more pronounced.
The full Holstein-Peierls model shows significantly differ- ent behavior. As also seen in the one-particle spectral function, spectra with strong Holstein coupling (g H = 2, green curves) are largely unaffected by Peierls disorder. The only qualitative change is the appearance of weak side peaks at ±ω P at T = 1000 K. In all cases with g H = 1 (blue curves), the addition of dynamic disorder significantly broadens and suppresses the zero-frequency peak, while shifting weight toward finite frequency peaks, particularly enhancing those located at odd multiples of ω H . The absence of even-numbered peaks is a coincidental destructive interference effect, because these phonon frequencies coincide with odd multiples of τ. In the presence of dynamic disorder with no Holstein coupling (red curves in Fig. 2(b), (d) and (f)), structure is observed at ω = 2τ and 4τ, while the ω = 0 peak is suppressed, in agreement with previous studies. 52 Summarizing our results on the finite-temperture AC conductivity of the Holstein-Peierls Hamiltonian, the lowfrequency features, which determine transport, are most strongly affected by Peierls-induced dynamic disorder and are virtually unaffected by the presence of moderate Holstein coupling. However, the Holstein coupling significantly modifies the high-frequency structure, which can be probed spectroscopically.

D. DC mobility
We now turn to a study of the DC mobility, which is the observable of greatest practical interest for device performance, presented in Fig. 3. First, we recall the single-carrier (N = 1) mobility in the absence of any electron-phonon coupling, corresponding to the Boltzmann transport equation 11,30,32 where t s is a Drude scattering time and v k = ∇ k ε(k). At low temperatures over which the band is effectively parabolic, a single electron in the canonical ensemble obeys ideal statistics v 2 = k B T /m * , producing a constant value which is the Drude mobility. At high temperatures, k B T τ, the average over the finite bandwidth is constant, J 2 ≈ 4πe 2 0 τ 2 a 2 , giving i.e. decaying like T −1 . The crossover between these two behaviors clearly occurs at k B T ≈ τ, as can be seen in Fig. 3(a), with g H , g P = 0 (red data).
With the addition of Holstein coupling, our theory of the mobility reduces to that of Ortmann et al. 30 and our analysis is correspondingly similar. The coherent contribution to the mobility has the same form as above, but with the temperaturedependent renormalizedτ(T ). Thus the low-temperature constant value is reduced by a factor of e −g 2 H , the crossover temperature is reduced, and the high-temperature decay is faster than T −1 . The mobility has an additional incoherent contribution that yields activated behavior at intermediate temperature, where the activation temperature decreases with increasing g H . At very high temperatures k B T g 2 ω H , we see the well-known T −3/2 power-law dependence of Marcus theory. 53 In Fig. 3(a), we extend our results to unrealistically high temperatures T ∼ 10 5 K only to show this asymptotic behavior.
Aside from the effects of band-narrowing -an overall reduction of mobility and faster than T −1 decay -the major signature of a high-frequency Holstein phonon is the presence of the activated regime. However, our numerical results indicate that for parameters relevant to a high-mobility OMC, the activated regime is only observed at temperatures exceeding 1000 K. This is due to the high frequency of the intramolecular vibrations, on the order of 150 − 200 meV. Therefore, we conclude that an intrinsic activated regime associated with polaron formation may not be observable in OMCs. However, this does not indicate the absence nor irrelevance of strong intramolecular electron-phonon coupling.
Lastly, we analyze the DC mobility in the presence of both Holstein and Peierls electron-phonon couplings, which is only enabled by the new theory developed here. In comparing Fig.  3(b) and (c) (moderate and strong Peierls coupling, g P = 1 and g P = 2) to (a) (no Peierls coupling), it is clear that the temperature dependence of the mobility is dramatically affected by dynamic disorder, except for the case of very strong Holstein coupling (green data). The reason for this latter insensitivity is illustrated in Fig. 4(a), which shows the mean value (solid lines) and standard deviation of the renormalized transfer integral under moderate (g P = 1) and strong (g P = 2) Peierls disorder (dark and light shaded regions, respectively). We see that for g H = 2, the mean value is nearly zero and, concomitantly, the fluctuations with these values of g P are negligible.
Having seen that Peierls dynamic disorder has little effect in the large g H limit, we will mainly compare the g H = 0 and g H = 1 cases. As previously stated, g H = 1 corresponds most closely to the estimated coupling strength in rubrene, [24][25][26][27] whereas the g H = 0 case reduces to the dynamic disorder/transient localization picture that has so successfully captured the transport behavior of high-mobility OMCs in the past. 20,23 In the temperature range of interest (10 2 -10 3 K), the g H = 0 and g H = 1 Holstein-Peierls mobilities are remarkably similar, exhibiting the well-known µ ∝ T −γ behavior, with γ ≈ 2 characteristic of dynamic disorder-induced transient localization. 13,22,23 For g P = 2 and g H = 0 or 1, our extracted  4. (a) Band-narrowed nearest-neighbor hoppingτ (solid lines), for which the bare value is τ = 100 meV. Shaded regions represent one standard deviation in the distribution of values thatτ can take under Peierls disorder with moderate coupling (g P = 1, darker shading) and strong coupling (g P = 2, lighter shading). (b) Energy-dependent localization length (solid lines) for zero Holstein coupling (g H = 0, red) and moderate Holstein coupling (g H = 1, blue) and thermally occupied density of states (shaded regions). Results are shown for T = 300 K and Peierls coupling g P = 2. (c) Temperature-dependent average localization length with Peierls coupling g P = 2. All other parameters are the same as in Fig. 1. value of γ ≈ 1.8 is slightly smaller than that determined from the slope of the mean-squared displacement in Refs. 20 and 24. This effect is potentially related to the spurious heating of quasiclassical dynamics 22,35 and which is significantly less influential in our formulation based on the Kubo formalism. The only major qualitative change induced by Holstein coupling is the activated behavior, which again occurs at unrealistically high temperatures exceeding 1000 K. In the dotted lines of Fig. 3, we show the mobility calculated in the absence of incoherent phonon effects, which is completely sufficient for experimentally relevant temperatures. This observation argues that the mobility of many OMCs should be simulated by a computational approach where the transfer integrals undergo a static but temperature-dependent renormalization due to the Holstein coupling to high-frequency intramolecular vibrations. These parameters are then used in a quasiclassical, time-dependent simulation of the electronic dynamics, following the usual prescriptions of dynamic disorder 20,21,24 or transient localization. 13,22,23 In Fig. 3(c), we compare our theoretical results to experimental Hall mobility data of the rubrene b-axis via Ref. 7. Both the g H = 0 and g H = 1 data show excellent agreement in terms of the absolute magnitude and temperature dependence of the mobility.
Before concluding, we seek to understand the origin of the mobility's insensitivity to the Holstein coupling strength g H , in cases where the Peierls coupling is nonnegligible. Figure 4(a) clearly shows that the magnitude of the renormalized transfer integral is significantly affected by g H , which would suggest a strong dependence. However, Fig. 4(a) also demonstrates that the band-narrowing reduces the variance of the transfer integrals, which lessens the degree of localization. To test the proposal, we analyze the localization properties of the electronic wavefunction averaged over disorder realizations and with statically renormalized transfer integrals. The energy-resolved coherent localization length, L coh (ω) can be obtained via the Thouless formalism, 54,55 whereε α =ε α (X P , P P ) are eigenvalues of the band-narrowed and statically disordered electronic Hamiltonian. The above quantities L coh (ω) and ρ coh (ω) are thermally averaged over the phonon degrees of freedom by Monte Carlo sampling as done in Eq. (8). In Fig. 4(b), we show this disorder-averaged energydependent localization length L coh (ω) without and with Holstein electron-phonon coupling (g H = 0 and g H = 1) at an example temperature of T = 300 K. We observe three main features: highly localized states at the band edges, states with a constant localization length of a few lattice constants in the middle of the band, and an unphysical delocalized state at the band center due to the purely off-diagonal nature of the disorder. Crucially, in the case with Holstein electron-phonon coupling (g H = 1), all of these features are compressed into a narrower energy spacing. The shaded region of Fig. 4(b) shows the thermally occupied density of states n(ω)ρ coh (ω), where n(ω) is the thermal occupancy. We see that the g H = 1 occupied density of states extends much further into the band, giving added weight to the more delocalized states. This yields a larger average localization length, which is defined by L coh = dωn(ω)ρ coh (ω)L coh (ω) and plotted in Fig. 4(c). Indeed, at all temperatures we see that the average localization length for a system with g H = 1 is larger than that with g H = 0. In other words, the addition of Holstein electronphonon coupling yields electronic states that are, on average, more delocalized. This parameter-specific effect increases the mobility and partially compensates for the decreased magnitude of the renormalized transfer integrals, producing a mobility that is relatively insensitive to the value of g H , as seen in Fig. 3.
More broadly, at all values of g H , we observe a maximum delocalization at low temperature and exponential decay with increasing temperature, as previously observed by Ciuchi and Fratini. 13 In our model, maximum delocalization does not occur at zero temperature because we include zeropoint energy in the Peierls phonons by sampling from the Wigner distribution, which allows for localized states even at zero temperature. Interestingly, at these low temperatures k B T < ω P (with Wigner sampling), the degree of disorder is temperature-independent. The temperature dependence of the localization length only comes from the thermal average over electronic states, where higher energy states have a larger localization length, causing the average localization length to increase with increasing temperature up to k B T ≈ ω P . To summarize, we have seen that the temperature-dependent mobility is determined by a subtle competition between a Holsteininduced renormalization of electronic parameters and their variance, the finite-temperature dynamics of quasiclassical Peierls modes, and the thermal occupancy of electronic states with different localization lengths.

IV. CONCLUSIONS
We have introduced a new approach to solving the Holstein-Peierls Hamiltonian for electron-phonon coupling in organic molecular crystals. We exploit the quasi-adiabatic nature of the intermolecular Peierls modes and treat them semiclassically. The lattice motions create a dynamically disordered landscape with broken translational symmetry that localizes the polaronic wavefunction. The intramolecular Holstein modes are accounted for by performing a Lang-Firsov polaron transformation, and separating the electron and Holstein phonon degrees of freedom through an anti-adiabatic finitetemperature mean-field approximation. The resulting dynamics are nonperturbative in both the electronic and electronphonon interactions. The polaronic dynamics exactly capture the limit of strong electron-phonon coupling, while reducing to the exact unitary electronic dynamics in the weak-coupling limit. We calculated the frequency-resolved spectral function and conductivity, as well as the zero-frequency DC mobility.
We found that finite-frequency observables such as the spectral function and optical conductivity are strongly affected by Holstein phonons even at moderate coupling strengths, especially at frequencies corresponding to multiples of the phonon energy. These features correspond to incoherent electron-phonon interactions, and are only enabled by a quantum mechanical treatment such as that presented here. Peierls disorder contributes peak-broadening and shifts spectral weight toward incoherent peaks. DC transport, in contrast to finite-frequency observables, is largely determined by dynamic disorder-induced localization of the coherent wavefunction. [20][21][22][23] The DC mobility is only weakly affected by the Holstein interaction in the form of temperature-dependent band-narrowing, i.e. effective mass renormalization. An incoherent activated regime due to Holstein phonons is observed, but only at experimentally irrelevant temperatures, at least for the frequency of intramolecular vibrations considered here and relevant for most OMCs. We further support this conclusion by showing that including only Holstein band-narrowing, and excluding incoherent features, can reproduce the total mobility to high accuracy up to temperatures of around 1000 K.
Interestingly, we find that coherent band-narrowing due to Holstein coupling has less of an effect on overall transport than might be expected, because the bandwidth reduction is accompanied by a reduction in the variance and thus a reduction in the nonlocal dynamic disorder. Such an observation is only made possible by the unified theory presented in this work. This finding further indicates that the success of theories with purely nonlocal electron-phonon coupling may be partially due to a serendipitous cancellation of effects which makes local coupling to intramolecular vibrations appear less significant. Previous studies have shown that dynamic disorder enhances hopping transport in 1D, but the mechanism was due to completely incoherent pathways and therefore unrelated to localization. 19 Because localization properties are strongly dependent on dimensionality and potential anisotropy, 56 the interplay of localization and bandrenormalization may be significantly different in higher dimensions, which we consider an interesting topic for future study.
Perhaps most importantly, the relative simplicity of our approach enables application to more complex models and realistic materials, where numerically exact techniques struggle. Because the theory leads to a time-dependent Hamiltonian with separable electronic and phononic degrees of freedom, we imagine the straightforward inclusion of electronic interactions, treated using conventional methods of electronic structure theory in the absence of phonons. This extension would be important in semiconductors at higher doping densities or in metals with strong electron-phonon interactions. Furthermore, a more sophisticated treatment of the electron-phonon interaction, for example via variational optimization of the Lang-Firsov transformation, is also possible. [57][58][59] Finally, we emphasize that the quasiclassical treatment of low-frequency nuclear motion is not limited to the study of linear electronphonon coupling or to harmonic nuclear degrees of freedom. Fully atomistic and anharmonic nuclear dynamics can be straightforwardly included; anharmonic effects are most important for low-frequency motions, which is precisely where the approach is valid. For example, in addition to OMCs, 60 nonlinear electron-phonon coupling and anharmonic effects have been shown to be important in the electronic properties of lead-halide perovskites 61,62 and for "phononic" control of nonequilibrium material properties. 63,64 Work along all of these directions is underway and will enable nonperturbative treatment of electron-nuclear interactions in realistic, complex materials.
In this appendix, we seek to benchmark the accuracy of the approximate Lang-Firsov treatment of the Holstein model by comparison with numerically exact results, using the same material parameters as considered in the main text. A large number of works have been devoted to the numerically exact treatment of electron-phonon problems, especially the Holstein model. In particular, variational techniques, 65 the density matrix renormalization group, 66,67 exact diagonalization of small clusters 37,38 and diagrammatic techniques 68 have been used to predict highly accurate or exact spectral properties for modestly sized 1D systems. Higher dimensional extensions have been explored within exact diagonalization 69 and dynamical mean field theory. 70 The later has been extended to finite doping and symmetry-broken phases. 71,72 Here, we employ exact diagonalization over the variational Hilbert space (EDVHS), which is described in the next section.

Exact diagonalization over variational Hilbert space
EDVHS can be used to obtain a numerically exact solution of the Holstein polaron ground state, 37 response functions 38 and nonequilibrium dynamics. 73,74 The construction of the variational space starts from an electron in a translationally invariant state with no phonon excitations. We generate new parent states by applying the Hamiltonian, Eqs. (1a) and (1c). By choosing the number of Hamiltonian operations N H , we control the maximal number of phonon quanta and the spatial extent of the polaron, which can be increased for convergence. After solving the ground state problem, the dynamical properties are evaluated using the continued fraction method. 75 All EDVHS calculations are converged with N H = 20 generations creating a Hilbert space with over 3 × 10 6 basis states. Lorentzian broadening with η = 10 meV for spectral functions and 15 meV for optical conductivity spectra were used in the presentation of spectra.

Spectral functions
The momentum-resolved spectral function for the Holstein model at zero temperature is presented in Fig. 5. We compare the numerically exact EDVHS results (top row) to the approximate Lang-Firsov treatment (bottom row). For intermediate coupling (g H = 1) we see good agreement, particularly in the dispersion of the quasiparticle peak and in the satellite structure near k = 0. These features are most relevant to the transport properties, which rely on the k = 0 contribution of the two-particle Green's function. 11 The satellite peaks, which are strictly dispersionless in the Lang-Firsov treatment, inherit some of the electronic dispersion in the exact results. At stronger electron-phonon coupling (g = 2), we see better agreement in the vibronic replica structure at all values of k; however the exact result again has slight dispersion in the peaks that is absent in the approximate result.
The positions of the spectral peaks are slightly shifted in the approximate result, which is due to an inaccurate groundstate energy of the N = 1 polaron problem. We examine this ground-state energy as a function of coupling strength g H in Fig. 6(a), which shows that the Lang-Firsov approach is accurate and both weak and strong coupling; the intermediatecoupling regime of g H = 1 − 2 is most challenging, but is qualitatively captured. We also note that a constant energy shift is not relevant for the conductivity, which depends only on energy differences at fixed electron number.
More important for transport is the quasiparticle effective mass at the bottom of the band. In Fig. 6(b) we look at the effective mass enhancement which, within the Lang-Firsov approximation, is the inverse of the band-narrowing factor; at zero temperature, m * /m 0 = τ/τ = e g H 2 . Again we see that the Lang-Firsov approach is accurate at small and large values of the coupling constant but slightly overestimates the mass enhancement at intermediate coupling.
Overall, we see excellent agreement between the approximate and exact spectral functions even in the weakly antiadiabatic, intermediate-coupling regime relevant to highmobility organic semiconductors. The approximate method is particularly accurate for the low-momentum, low-energy features that are relevant for conductivity and DC transport at low density.

Optical conductivity
Next we examine the frequency-dependent optical conductivity, σ(ω), which is presented in the main panels of Fig. 7. The conductivity can be separated into a zero-frequency Drude contribution Dδ(ω) and a regular finite-frequency contribution σ reg (ω); the latter is shown in the insets. Overall, the structure of the conductivity is well-reproduced by the Lang-Firsov treatment, including peak locations and lineshapes. The exact results show a complex mixture of coherent and incoherent features that contributes fine structure, which is absent in the approximate spectrum. The approximate Lang-Firsov approach predicts a structure that is dominated by regularly spaced vibronic replicas similar to that seen in the single-site limit. 11 The low-frequency behavior, which determines transport properties, is qualitatively reproduced by the approximate treatment. For the case with moderate coupling (g H = 1), the Drude weights are in satisfactory agreement: D = 0.68 (exact) and D = 0.37 (approximate). At strong coupling (g H = 2), the agreement is slightly worse: D = 0.054 (exact) and D = 0.018 (approximate). At low density, the Drude weight is related to the effective mass via D ∝ 1/m * and the underestimation of the Drude weight is consistent with the overestimation of the mass enhancement seen in Fig. 6(b).
To summarize our comparison of the conductivity at zero temperature, the low-frequency features show reasonable accuracy, including the zero-frequency Drude contribution and the locations and lineshapes of the first few excited-state peaks. Crucially, the nonperturbative Lang-Firsov approach captures exactly the weak-and strong-coupling limits, enabling far richer physics than those weak-coupling perturbation theory or the historically prevalent narrow-band approximation.