Introduction

The last two decades have witnessed the rapid development of high-mobility crystalline organic semiconductors1,2,3. The first proposition of using Marcus’ semiclassical hopping model coupled with density functional theory to design high mobility molecules has been very successful and popular4. It can be regarded as the strong local electron–phonon coupling (EPC) limit, which was later improved by considering quantum nuclear effect and delocalization effect5,6, through which, isotope effect is found to be always negative and the dynamic disorder does not play appreciable role with some experimental supports7,8,9,10. However, such local EPC picture is challenged by the recently established transient localization (TL) model which invokes nonlocal EPC11,12,13,14. In due course, the molecular design principles derived from TL such as suppressing intermolecular vibration, which are quite different from the local picture, have proved successful in a number of experiments15,16. The applicability of the TL picture, nevertheless, is restricted to the regime of moderate transfer integral (electronic coupling) V and strong nonlocal EPC. Since EPC is a complicated many-body problem, it is highly desirable to present a general transport picture taking both local and nonlocal EPC into considerations in a rigorous way, instead of uncontrolled approximations.

Many recent efforts have been devoted to developing approximate methods that are able to portray a broader parameter space, including the band-like (BL) conduction regime17,18. Besides, unlike in the case of Holstein model in which it is beyond doubt that local EPC represents an obstacle for carrier diffusion19,20,21, how nonlocal EPC affects mobility does not have a definitive answer and the interplay between the local and nonlocal EPC is unclear. Early theoretical treatments for carrier mobility in crystalline organic semiconductors such as the Munn-Silbey approach and the polaron transformation often reach the conclusion that the nonlocal EPC leads to phonon-assisted (PA) transport and enhances mobility22,23,24, in sheer contrast with the basic starting point of the TL scenario. The findings of several numerical studies on the carrier mobility of organic semiconductors also contradict with the TL theory8,25. In principle, PA, TL, and BL are all possible mechanisms for charge transport with nonlocal EPC, valid at their respective parameter regimes, yet a universal theoretical treatment for the role of nonlocal EPC is not available due to the complex many-body electron–phonon interaction.

In this work, we present a nearly exact study of the charge transport mechanism in the Holstein–Peierls model using the time-dependent finite temperature matrix product state (MPS) formalism26,27. By studying EPC effect on the carrier mobility, mean free path, optical conductivity, and one-particle spectral function, we have located the PA, TL, and BL regimes simultaneously on the transfer integral – nonlocal EPC strength plane. We have also identified an intermediate regime where none of the existing pictures is truly applicable, as a generalization of the hopping-band crossover in the Holstein model.

Results

System Hamiltonian and Kubo Formula

We take the following one-dimensional Holstien–Peierls model with nearest-neighbor interaction and periodic boundary condition:

$$\hat{H} =\; {\hat{H}}_{{\rm{e}}}+{\hat{H}}_{{\rm{ph}}}+{\hat{H}}_{{\rm{e-ph}}}\\ {\hat{H}}_{{\rm{e}}} =-V\mathop{\sum}\limits_{n}({c}_{n+1}^{\dagger }{c}_{n}+{c}_{n}^{\dagger }{c}_{n+1})\\ {\hat{H}}_{{\rm{ph}}} = \mathop{\sum}\limits_{n,m}{\omega }_{m}{b}_{n,m}^{\dagger }{b}_{n,m}+\mathop{\sum}\limits_{n}{\omega }_{\theta }{b}_{n,\theta }^{\dagger }{b}_{n,\theta }\\ {\hat{H}}_{{\rm{e-ph}}} = \mathop{\sum}\limits_{n,m}{g}_{m,{\rm{I}}}{\omega }_{m}({b}_{n,m}^{\dagger }+{b}_{n,m}){c}_{n}^{\dagger }{c}_{n}\\ +\mathop{\sum}\limits_{n}{g}_{\theta ,{\rm{II}}}{\omega }_{\theta }({b}_{n,\theta }^{\dagger }+{b}_{n,\theta })({c}_{n+1}^{\dagger }{c}_{n}+{c}_{n}^{\dagger }{c}_{n+1})$$
(1)

where c† (c) and b† (b) are the creation (annihilation) operator for electron and phonon respectively, and V is the intermolecular transfer integral. The electronic motion is limited to single-electron manifold. ωm and gm,I are the frequency and the dimensionless EPC constant of the mth intramolecular vibration mode. ωθ and gθ,II are the intermolecular vibration counterparts. ℏ is set to 1. The thermal fluctuation ΔV is related to intermolecular coupling constant gθ,II by28:

$${{\Delta }}V={g}_{\theta ,{\rm{II}}}{\omega }_{\theta }\sqrt{\coth \frac{{\omega }_{\theta }}{2{k}_{B}T}}$$
(2)

The one-dimensional model in Eq. (1) is an approximation to realistic organic semiconductors, which typically adopt two-dimensional transport network8,13,29. In the Supplementary Fig. 4 we demonstrate that this approximation is valid for anisotropic materials by comparing the simulated one-particle spectral function with experimental angle resolved ultraviolet photoemission spectra (ARUPS)30, and at the end of the section we go beyond the one-dimensional model to discuss the isotropy effect on the different transport regimes.

In order to elucidate how nonlocal EPC affects charge transport at different transport regimes we focus on the role of transfer integral V and nonlocal EPC constant gθ,II (or equivalently ΔV at a given T). Other parameters are fixed throughout this paper unless otherwise specified with values drawn from representative organic semiconductors. In organic semiconductors, the most common values of V and ΔV range from 10 meV to 150 meV and from 10 meV to 60 meV, respectively31,32,33. The intramolecular vibration frequency ωm and local EPC constant gm,I are taken from our previous DFT calculations for rubrene and the total 3N − 6 normal vibration modes are reduced to four effective modes6,21, namely ωm = 26 meV, 124 meV, 167 meV, and 198 meV, with the corresponding dimensionless gm,I 0.83, 0.26, 0.34, and 0.37, respectively. The intermolecular vibration frequency ωθ is set to be 50 cm−1 (6.2 meV) as commonly adopted in literature13,17,34. The system is translational-invariant and we do not consider the effect of both diagonal and off-diagonal static disorder here35.

The carrier mobility is obtained via the Kubo formula36:

$$\mu =\frac{1}{{k}_{B}T{e}_{0}}\int_{0}^{\infty }\left\langle \hat{j}(t)\hat{j}(0)\right\rangle {\mathrm{d}}t=\frac{1}{{k}_{B}T{e}_{0}}\int_{0}^{\infty }C(t){\mathrm{d}}t$$
(3)

where for the Holstein–Peierls Hamiltonian in Eq. (1) the current operator \(\hat{j}\) takes the form:

$$\hat{j}=\frac{{e}_{0}R}{i}\mathop{\sum}\limits_{n}\left[-V+{g}_{\theta ,{\rm{II}}}{\omega }_{\theta }({b}_{n,\theta }^{\dagger }+{b}_{n,\theta })\right]({c}_{n+1}^{\dagger }{c}_{n}-{c}_{n}^{\dagger }{c}_{n+1})$$
(4)

Here R is the intermolecular distance and is set to 7.2 Å as in the case of rubrene crystal. Although we have treated the model as a closed system, in our study the recurrence problem is not severe and C(t) in general rapidly decays to nearly zero, except when both V and ΔV are small. In such cases we resort to a more strict model with 10 modes in total for the convergence of C(t). Lying at the heart of our calculation is the evaluation of the current–current correlation function \(C(t)=\left\langle \hat{j}(t)\hat{j}(0)\right\rangle\), which is achieved by the time-dependent MPS formalism21,26,27. In most of our simulations, the number of molecules in the periodic one-dimensional chain is 21 and the virtual bond dimension is 80. More details on the model with 9 intramolecular modes as well as numerical convergence check on system size and MPS parameters are included in Supplementary Table 1, Supplementary Fig. 1, and Supplementary Fig. 2.

Carrier mobility

Firstly, we analyze the role of local and nonlocal EPC on different parameter regimes by comparing the mobility calculated based on the Holstein–Peierls model (μH-P) with the mobility calculated based on pure Holstein model (μH) and pure Peierls model (μP) at 300 K. We have also included the mobility derived from Matthiessen’s rule (1/μM = 1/μH + 1/μP), presumably valid in the BL regime, because in the BL regime both local and nonlocal EPC can be considered as independent scattering sources and the total scattering rate of the wave-like electronic motion is the sum of both scattering rates. The overall results are shown in Fig. 1. When V = 5 meV, μH-P is generally higher than μH, which implies that nonlocal EPC is beneficial to charge transport. This is considered to be a signature of the PA picture. However, this behavior quickly vanishes as V increases from 5 meV to 20 meV. At the intermediate range of V, e.g. from 45 meV to 120 meV, μH-P could be higher than μP for larger ΔV. That local EPC could enhance instead of reduce mobility is quite counter-intuitive. Such unusual behavior can be best understood by TL picture, in which the quantum coherent interference responsible for Anderson localization can be damaged or destroyed by local EPC as the dephasing noise. The mechanism is studied in detail by means of open system dynamics for systems with static disorder37,38. In addition, the “band width narrowing” caused by local EPC could allow the carrier to thermally access much delocalized states17. These lead to the local EPC enhanced mobility. Upon further increasing V to 150 meV, the TL scenario also becomes less significant. Instead, it is found that μM coincides with μH-P remarkably well, serving as a piece of evidence for band-like transport.

Fig. 1: Carrier mobility at 300 K calculated based on the Holstein–Peierls model (μH-P), Holstein model (μH), Peierls model (μP), and Matthiessen’s rule (μM) with various transfer integral V and transfer integral fluctuation ΔV.
figure 1

From a to i the transfer integrals are 5, 10, 20, 30, 45, 60, 90, 120, and 150 meV respectively. Other parameters relevant to the carrier mobility such as local EPC constants are fixed with values taken from rubrene. The parameters for the Holstein (Peierls) model is the same as that of Holstein–Peierls model except that the nonlocal (local) EPC is left out.

Full size image

Mean free path and optical conductivity

Although in the BL regime μM is expected to be a good approximation for μH-P, μM ≈ μH-P alone is not a sufficient condition for band-like transport, and we additionally employ the Mott-Ioffe-Regel limit for the determination of the BL regime. The carrier mean free path lmfp is estimated as lmfp = vτ with the group velocity v and relaxation time τ evaluated by39:

$$v =\sqrt{\left\langle \hat{j}(0)\hat{j}(0)\right\rangle }/{e}_{0}\\ \tau = \int_{0}^{\infty }\left|\frac{{\rm{Re}}C(t)}{{\rm{Re}}C(0)}\right|{\mathrm{d}}t.$$
(5)

And the calculated lmfp in the (V, ΔV) plane at 300 K is shown in Fig. 2a. The overall tendency of lmfp matches well with the carrier mobility of the Holstein–Peierls model μH-P in Fig. 1. The region where lmfp > R is colored with blue in Fig. 2a and it lies within the large V, small ΔV limit, and in agreement with common perception. Another possible criteria for BL conduction is the appearance of Drude-like peak in the per carrier optical conductivity:

$$\frac{\sigma (\omega )}{n{e}_{0}}=\frac{1-{e}^{-\omega /{k}_{B}T}}{\omega }\int_{0}^{\infty }C(t){e}^{i\omega t}{\mathrm{d}}t$$
(6)

which is illustrated in Fig. 2b. In the case of V = 150 meV without nonlocal EPC, a broad Drude-like peak appears near ω = 0. Upon adding a small amount of nonlocal EPC, the Drude-like peak becomes invisible. However, the optical conductivity is still significantly different from the ΔV = 60 meV cases, in which a localization peak at ω ≈ 200 meV characteristic for the TL regime18 is present.

Fig. 2: Further analysis of the transport regimes.
figure 2

a Carrier mean free path lmfp/R at 300 K for the Holstein-Peierls model at various transfer integral V and transfer integral fluctuation ΔV. In the blue region (bottom right) the carrier mean free path exceeds the lattice constant R. b Per carrier optical conductivity of the Holstein–Peierls model at various transfer integral V and transfer integral fluctuation ΔV. c, d Correlation functions obtained from our simulation with Holstein-Peierls model C(t)H-P, pure Holstein model C(t)H and pure Peiels model C(t)P as well as the correlation function obtained from phonon-assisted transport theory24C(t)PA and transient localization theory40C(t)TL for c V = 5 meV, ΔV = 5 meV and d V = 90 meV, ΔV = 40 meV.

Full size image

Connection with semi-analytical theories

The PA regime and the TL regime can be further confirmed by semi-analytical results. In Fig. 2c we compare C(t)H-P of our numerical simulation with C(t)PA of phonon-assisted charge transport theory24:

$$\mu =\frac{{e}_{0}{R}^{2}}{{k}_{B}T}\int_{-\infty }^{\infty }[{V}^{2}+{({g}_{\theta ,{\rm{II}}}{\omega }_{m})}^{2}{{{\Phi }}}_{\theta }(t)]{e}^{-{{\Gamma }}(t)}{\mathrm{d}}t\\ {{\Gamma }}(t) =2\mathop{\sum}\limits_{m}{g}_{m,{\rm{I}}}^{2}[1+2{N}_{m}-{{{\Phi }}}_{m}(t)]+4{g}_{\theta ,{\rm{II}}}^{2}[1+2{N}_{\theta }-{{{\Phi }}}_{m}(t)]\\ {{{\Phi }}}_{m} =(1+{N}_{m}){e}^{-i{\omega }_{m}t}+{N}_{m}{e}^{i{\omega }_{m}t}\\ {N}_{m} =\frac{1}{{e}^{{\omega }_{m}/{k}_{B}T}-1}$$
(7)

The parameters are V = 5 meV and ΔV = 5 meV. For both real and imaginary part the two curves are in excellent agreement, and show significant increase with respect to the correlation function with only local EPC C(t)H. Thus we can confidently conclude that in this parameter regime the transport mechanism can be understood as phonon-assisted transport. We note that the derivation of Eq. (7) employs narrow-band approximation, which is valid in the small V limit. When V = 90 meV and ΔV = 40 meV shown in Fig. 2d, the correlation function given by TL theory with relaxation time approximation40 C(t)TL is in agreement with our calculation based on pure Peierls model C(t)P. The observation implies that in this regime the TL theory can successfully account for the transport property of the pure Peierls model, from which the TL theory is derived. If the Holstein coupling is included, the correlation function C(t)H-P exhibits significant difference from C(t)P and C(t)TL, however, the integrated mobility turns out to be rather insensitive to Holstein coupling in this regime (Fig. 1g). We note that it is possible to integrate Holstein coupling in the transient localization theory if the intramolecular vibration frequency is much smaller than the transfer integral33,41, however such scheme is not employed in this work because in most cases ωm is at the same order with V. We believe it is suitable to ascribe the V = 90 meV and ΔV = 40 meV case as TL, because although the TL theory with relaxation time approximation may not correctly produce the correlation function for realistic materials with Holstein coupling, the picture provided by the theory serves as a nice starting point for further analysis.

Effect of local EPC strength

In order to investigate how local EPC strength will affect the results in Fig. 1, we have further calculated the carrier mobility of the Holstein–Peierls model with stronger local EPC. More specifically, the values of gm,I are multiplied by \(\sqrt{2}\) so that the respective reorganization energies \({g}_{m,{\rm{I}}}^{2}{\omega }_{m}\) are doubled. The results are illustrated in Fig. 3. In the small V limit shown in Fig. 3a, the PA mechanism prevails, in agreement with the results in Fig. 1. However, from Fig. 3b it can be seen that with enlarged local EPC strength the PA picture remains valid even if V becomes as large as 20 meV, in contrast to the V = 20 meV results presented in Fig. 1, indicating the expansion of the PA region. Accordingly, the TL region is diminished, as can be inferred from Fig. 3c, d by noting that the parameter space in which μH-P ≥ μP is satisfied is smaller than that of Fig. 1. In Fig. 4 we show the carrier mean free path lmfp/R with increased local EPC strength. When V is relatively large and ΔV is relatively small, namely in the BL regime, lmfp/R with increased local EPC strength is generally smaller than that of the original local EPC strength. This observation implies that the BL region in the (V, ΔV) plane moves toward the even larger V area (V > 150 meV). Our findings are in agreement with physical instinct because in the large local EPC limit the hopping mechanism dominates and Eq. (7) is a good approximation for mobility. Another factor that might affect charge transport is the distribution of intramolecular vibration frequencies with fixed total reorganization energy. The problem is equivalent to the isotope effect problem and recent studies on the rubrene molecule concluded negative isotope effect6,9,21.

Fig. 3: Carrier mobility at 300 K calculated based on the Holstein–Peierls model (μH-P), Holstein model (μH), Peierls model (μP), and Matthiessen’s rule (μM) with enlarged local EPC at various transfer integral V and transfer integral fluctuation ΔV.
figure 3

From a–d the transfer integrals are 5, 20, 60, and 150 meV respectively.

Full size image
Fig. 4: Carrier mean free path with enlarged local EPC.
figure 4

Carrier mean free path lmfp/R at 300 K for the Holstein–Peierls model with enlarged local EPC at various transfer integral V and transfer integral fluctuation ΔV.

Full size image

One-particle spectral function

To further analyze the charge transport properties in the regimes implied by Figs. 1 and 2, we calculated the momentum resolved one-particle spectral function:

$$A(k,\omega )=\frac{1}{N\pi }\mathop{\sum }\limits_{mn}^{N}{e}^{ikR(m-n)}\int_{0}^{\infty }\left\langle {c}_{m}(t){c}_{n}^{\dagger }(0)\right\rangle {e}^{i\omega t}{\mathrm{d}}t$$
(8)

for nine sets of representative parameters at 300 K in Fig. 5. A Lorentzian broadening with η = 5 meV is applied for a smooth spectra. When V = 5 meV and ΔV = 10 meV (Fig. 5a), the spectral function exhibits dispersionless bound states separated by the intramolecular vibration frequencies ωm, which marks the formation of small polaron. On the contrary, when V = 150 meV and ΔV = 10 meV (Fig. 5c) an intense quasiparticle peak is observed near k = 0 and the overall shape of the spectra resembles the dispersion of free electron \(E(k)=-2V\cos kR\). In the TL regime represented by V = 60 meV and ΔV = 60 meV (Fig. 5e) the signature of either small polaron or delocalized state is almost completely smeared out. In combination with the limited carrier mean free path in this regime, it can be deduced that the charge carrier is localized by nonlocal EPC instead of local EPC. The same “blurred” spectral function is observed for other sets of typical parameters in the TL regime (Fig. 5d, f). With moderate V and ΔV shown in Fig. 5g, the spectral function exhibits none of the typical features described above. Namely, although the spectral function does not manifest the formation of small polaron or delocalized states, the peak intensity is still strong enough to be discernible from the TL regimes cases (Fig. 5d–f). In the absence of the local EPC (Fig. 5h, i), the quasiparticle peak has more intensity, implying the disruption of quantum coherent with the addition of local EPC. Combined with the mobility data shown in Fig. 1, it can be inferred that in the TL regime, the effect of the disruption is to alleviate the localization caused by nonlocal EPC, leading to increased mobility (Fig. 5e, h), while in the BL regime, on the contrary, the disruption scatters charge carrier and reduces mobility (Fig. 5c, i). The horizontal peak at the center of the band in Fig. 5h is a result of the pure nonlocal EPC in the Peierls model and is expected to vanish upon the inclusion of infinitesimal local EPC42.

Fig. 5: Spectral function for the Holstein-Peierls model and pure Peierls model.
figure 5

Spectral function at 300 K for the Holstein–Peierls model (a–g) and the pure Peierls model (h, i). Panel h and i has the same parameter as panel e and c respectively except that in Panel h and i the local EPC is left out.

Full size image

The isotropy effect

In a number of recent works it is established that dimensionality plays an indispensable role in the charge transport process, especially when dynamic disorder is taken into consideration8,13,43. To study the isotropy effect beyond the one-dimensional model, we employ a quasi-two-dimensional ladder Holstein–Peierls Hamiltonian:

$$\hat{H} = \, {\hat{H}}_{{\rm{e}}}+{\hat{H}}_{{\rm{ph}}}+{\hat{H}}_{{\rm{e-ph}}}\\ {\hat{H}}_{{\rm{e}}} = -{V}_{1}\mathop{\sum}\limits_{l=1,2}\mathop{\sum}\limits_{n}({c}_{l,n+1}^{\dagger }{c}_{l,n}+{c}_{l,n}^{\dagger }{c}_{l,n+1})-{V}_{2}\mathop{\sum}\limits_{n}({c}_{0,n}^{\dagger }{c}_{1,n}+{c}_{1,n}^{\dagger }{c}_{0,n})\\ {\hat{H}}_{{\rm{ph}}} = \mathop{\sum}\limits_{l,n,m}{\omega }_{m}{b}_{l,n,m}^{\dagger }{b}_{l,n,m}+\mathop{\sum}\limits_{l,n}{\omega }_{\theta }{b}_{l,n,\theta }^{\dagger }{b}_{l,n,\theta }\\ {\hat{H}}_{{\rm{e-ph}}} = \mathop{\sum}\limits_{l,n,m}{g}_{m,{\rm{I}}}{\omega }_{m}({b}_{l,n,m}^{\dagger }+{b}_{l,n,m}){c}_{l,n}^{\dagger }{c}_{l,n}\\ +\mathop{\sum}\limits_{l,n}{g}_{\theta ,{\rm{II}}}{\omega }_{\theta }({b}_{l,n,\theta }^{\dagger }+{b}_{l,n,\theta })({c}_{l,n+1}^{\dagger }{c}_{l,n}+{c}_{l,n}^{\dagger }{c}_{l,n+1})$$
(9)

here V1 and V2 represent the electronic coupling at the high-mobility direction and the low-mobility direction respectively. The intermolecular vibration at the V2 direction is neglected for simplicity. The setup, while still approximate compared to a full-fledged two-dimensional model, is believed to be reasonable for anisotropic materials (V2 ≪ V1) and can at least partially capture the dimensionality effect. In Fig. 6 we present the computed correlation function C(t) and mobility μ based on the model for several typical values of V1 and ΔV1. In the hopping limit shown in Fig. 6a, e, it is found that carrier mobility is irrelevant to the isotropy effect, because in this limit V2 does not affect the hopping process at V1 direction. In the phonon-assisted transport regime shown in Fig. 6b, f, μ is rather insensitive to isotropy effect. In the band-like regime shown in Fig. 6c, g, we find that isotropy effect tends to slightly increase mobility. Lastly, in the transient localization regime shown in Fig. 6d, h, it is observed that carrier mobility is susceptible to the isotropy effect. By increasing V2 from 0 to 0.2V1, the mobility increases by ~40%. Such increase may appear difficult to understand as C(t) in Fig. 6d does not seem to vary much. This is because transient localization implies that during the integration of C(t) the positive and the negative part of C(t) are canceled out, and thus minor changes in C(t) will result in a big difference in mobility. Our result is generally in agreement with previous literatures8,13. Based on these findings we can conclude that when the isotropy of the system is increased, the band region tend to expand while the transient localization regime would tend to diminish18. Physically, the first conclusion can be understood by enlarged bandwidth in two dimension and the second conclusion can be understood by considering that Anderson localization length for two dimension is larger than that in one dimension44.

Fig. 6: The isotropy effect for several typical values of V1 and ΔV1.
figure 6

a–d are the correlation functions C(t) and e–h are the corresponding mobilities μ. For a and e, V1 = 5 meV and ΔV1 = 0 meV; For b and f, V1 = 5 meV and ΔV1 = 20 meV; For c and g, V1 = 90 meV and ΔV1 = 0 meV; For d and h, V1 = 90 meV and ΔV1 = 40 meV.

Full size image

General charge transport regime diagram

Based on the EPC effect on the carrier mobility, the mean free path, the optical conductivity, and the one-particle spectral function, we are able to sketch a schematic “regime diagram” for the charge transport mechanisms as shown in Fig. 7. The PA regime is determined by μH-P > μH, short lmfp, C(t)H-P ≈ C(t)PA and narrow bound state states in the spectral function. The TL regime is determined by μH-P ≥ μP, intermediate lmfp, localization peak in optical conductivity, C(t)P ≈ C(t)TL and a “smeared out” spectral function. The BL regime is determined by μH-P ≈ μM, large lmfp, Drude-like peak in optical conductivity and sharp quasiparticle peak in the spectral function. In Fig. 7 we use μH-P > μH, μH-P ≥ μP and lmfp > R for the boundaries of the PA regime, TL regime, and BL regime respectively, and using other indicators such as the appearance of Drude-like peak for the BL regime may shift the boundaries to some extent but the general picture remains intact. On this (V, ΔV) plane we are also able to identify an “intermediate” regime that lies among the PA regime, TL regime, and BL regime. In this regime, μH-P is significantly lower than both μH and μP, and the carrier mean free path is still less than the lattice constant, forbidding the band description. In fact, for the pure Holstein model case (ΔV = 0), the intermediate regime simply degenerates into the canonical hopping-band crossover. The crossover from the BL regime to the TL regime has also been reported by introducing transient localization correction to band transport18. The gray solid arrows and gray dashed arrows in Fig. 7 indicate the shift of the boundaries upon increasing local EPC and increasing electronic coupling isotropy respectively, based on Figs. 3, 4 and 6. To provide a rough intuition of the distribution of parameters for realistic organic semiconductors on this (V, ΔV) plane, in Fig. 7 we have also marked the value of V and ΔV for several common organic semiconductors using reported values from recent literature32. These materials are pMSB, pyrene, naphthalene, perylene, anthracene, DATT, rubrene, and pentacene from left to right. We note that the colors at the location of the markers do not represent a prediction of the charge transport mechanism for the corresponding materials because the materials do not necessarily share the same local EPC coupling strength, transport network and intermolecular vibration frequency with the parameters used in this work.

Fig. 7: A schematic “regime diagram” showing the phonon-assisted transport (PA) regime, transient localization (TL) regime, band-like (BL) regime, and intermediate regime on the (V, ΔV) plane for the carrier mobility of the Holstein–Peierls model.
figure 7

The hopping regime is achieved in the ΔV = 0 limit of the PA regime. Gray solid arrows show qualitatively the shift of the boundaries when local EPC increases. Gray dashed arrows show qualitatively the shift of the boundaries when transport network changes from one dimension to quasi-two-dimension or equivalently when electronic coupling isotropy increases. The green dots represent the V and ΔV of several common organic semiconductors.

Full size image

Discussion

In this work, we present a nearly exact theoretical study of the carrier mobility in Holstein–Peierls model with parameters relevant to organic semiconductors. By carefully investigating the effect of both local and nonlocal EPCs on the carrier mobility μ, mean free path lmfp, per carrier optical conductivity \(\frac{\sigma (\omega )}{n{e}_{0}}\), and one-particle spectral function A(k, ω), we are able to identify the PA regime, TL regime and BL regime on the (V, ΔV) plane. The PA regime features μH-P > μH, short lmfp, and narrow bound states in the spectral function. The TL regime features μH-P ≥ μP, intermediate lmfp, localization peak in optical conductivity and a “smeared out” spectral function. And the BL regime features μH-P ≈ μM, large lmfp, Drude-like peak in optical conductivity and sharp quasiparticle peak in the spectral function. The semiclassical Marcus hopping regime is found to be around the corner of small V and ΔV. Furthermore, some of the parameters in the (V, ΔV) plane are recognized to lie in an intermediate regime that does not exhibit the typical features described above, and this regime can be considered as a generalization of the hopping-band crossover regime in the Holstein model. We find that when increasing local EPC strength, the PA regime will expand while the TL and BL regime will diminish. When going from one dimension to quasi-two-dimension, the TL regime will diminish and the BL regime will expand. It should be noted that the localization effect due to static disorder is not considered here, which deserves further investigation.

Methods

Matrix product states

The evaluation of the current–current correlation function \(C(t)=\left\langle \hat{j}(t)\hat{j}(0)\right\rangle\) is performed by time-dependent matrix product states through imaginary and real time propagation. The matrix product states method represent the wavefunction of many-body system as the product of a series of matrices26:

$$\left|{{\Psi }}\right\rangle =\mathop{\sum}\limits_{\{a\},\{\sigma \}}{A}_{{a}_{1}}^{{\sigma }_{1}}{A}_{{a}_{1}{a}_{2}}^{{\sigma }_{2}}\cdots {A}_{{a}_{N-1}}^{{\sigma }_{N}}\left|{\sigma }_{1}{\sigma }_{2}\cdots {\sigma }_{N}\right\rangle$$
(10)

\(\left|{\sigma }_{i}\right\rangle\) is the basis for each degree of freedom. \({A}_{{a}_{i-1}{a}_{i}}^{{\sigma }_{i}}\) are matrices in the chain connected by indices ai. {⋅} in the summation represents the contraction of the respective connected indices, and N is the total number of degrees of freedom (DOFs) in the system. The dimension of ai is called (virtual) bond dimension, while the dimension of σi is called physical bond dimension. In principle, the time-dependent algorithms for MPS27 is able to solve the time-dependent Schrödinger equation in an exact manner if the bond dimension is infinite. In practice, the accuracy of the method can be systematically improved by using a larger bond dimension, until convergence of interested physical observables within arbitrary convergence criteria.

Finite temperature algorithm

The finite temperature effect is taken into account through thermal field dynamics, also known as the purification method26,45. The thermal equilibrium density matrix of any mixed state in physical space P can be expressed as a partial trace over an enlarged Hilbert space P ⊗ Q, where Q is an auxiliary space chosen to be a copy of P. The thermal equilibrium density operator can then expressed as a partial trace of the pure state Ψβ in the enlarged Hilbert space over the Q space:

$${\hat{\rho }}_{\beta }=\frac{{e}^{-\beta \hat{H}}}{Z}=\frac{{{\rm{Tr}}}_{Q}|{{{\Psi }}}_{\beta }\rangle \langle {{{\Psi }}}_{\beta }|}{{{\rm{Tr}}}_{PQ}|{{{\Psi }}}_{\beta }\rangle \langle {{{\Psi }}}_{\beta }|}$$
(11)

and the pure state \(|{{{\Psi }}}_{\beta }\rangle\) represented as an MPS is obtained by the imaginary time propagation from the locally maximally entangled state \(\left|I\right\rangle ={\sum }_{i}{\left|i\right\rangle }_{P}{\left|i\right\rangle }_{Q}\) to β/2 in the one electron manifold:

$$|{{{\Psi }}}_{\beta }\rangle ={e}^{-\beta \hat{H}/2}|I\rangle .$$
(12)

To calculate C(t), \(|{{{\Psi }}}_{\beta }\rangle\), and \(\hat{j}(0)|{{{\Psi }}}_{\beta }\rangle\) are propagated in real time to obtain \({e}^{-i\hat{H}t}|{{{\Psi }}}_{\beta }\rangle\) and \({e}^{-i\hat{H}t}\hat{j}(0)|{{{\Psi }}}_{\beta }\rangle\) and then C(t) is calculated by:

$$C(t)=\langle {{{\Psi }}}_{\beta }|{e}^{i\hat{H}t}\hat{j}(0){e}^{-i\hat{H}t}\hat{j}(0)|{{{\Psi }}}_{\beta }\rangle /Z.$$
(13)

Here the current operator \(\hat{j}(0)\) is represented as an MPO and inner-product for \(|{{{\Psi }}}_{\beta }\rangle\) includes tracing over both P space and Q space. The construction of the MPOs is performed in an automatic and optimal fashion through our recently proposed algorithm46. Note that different from the simulation of diffusion dynamics, the initial state of the formulation does not require electronic excitation from the zero electron manifold. In principle, both imaginary and real time propagation can be carried out by any time evolution methods available to matrix product states27. In this work, we use the time-dependent variational principle based projector splitting time evolution scheme47,48, which is found to be relatively efficient and accurate combined with graphic processing unit (GPU) in our recent work49.