Coherent quantum transport in disordered systems: II. Temperature dependence of carrier diffusion coefficients from the time-dependent wavepacket diffusion method

The time-dependent wavepacket diffusion method for carrier quantum dynamics (Zhong and Zhao 2013 J. Chem. Phys. 138 014111), a truncated version of the stochastic Schrödinger equation/wavefunction approach that approximately satisfies the detailed balance principle and scales well with the size of the system, is applied to investigate the carrier transport in one-dimensional systems including both the static and dynamic disorders on site energies. The predicted diffusion coefficients with respect to temperature successfully bridge from band-like to hopping-type transport. As demonstrated in paper I (Moix et al 2013 New J. Phys. 15 085010), the static disorder tends to localize the carrier, whereas the dynamic disorder induces carrier dynamics. For the weak dynamic disorder, the diffusion coefficients are temperature-independent (band-like property) at low temperatures, which is consistent with the prediction from the Redfield equation, and a linear dependence of the coefficient on temperature (hopping-type property) only appears at high temperatures. In the intermediate regime of dynamic disorder, the transition from band-like to hopping-type transport can be easily observed at relatively low temperatures as the static disorder increases. When the dynamic disorder becomes strong, the carrier motion can follow the hopping-type mechanism even without static disorder. Furthermore, it is found that the memory time of dynamic disorder is an important factor in controlling the transition from the band-like to hopping-type motions.

Although a rigorous theoretical prediction of transport dynamics is still lacking, under certain limits, the transport properties are well understood. For example, in the case of static disorder, Anderson localization occurs in one-or two-dimensional systems for any finite amount of disorder, leading to a complete lack of transport in the long-time limit [24][25][26]. For homogeneous systems with strong dynamic disorder, the transport is assumed to follow the hopping mechanism where the carrier hops incoherently between adjacent molecular sites [27]. In this case, Fermiʼs golden rule (FGR) (small polaron theory) [21,28] was originally proposed to describe the hopping rate in organic crystals, and the well-known Marcus formula [29,30] corresponds to its high-temperature limit. The solvent dynamics can be further incorporated as the electronic coupling becomes strong (see, for instance, [17], and reference therein). In homogeneous systems with small dynamic disorder, however, the transport exhibits a coherent and delocalized behavior, and the well-established band-like theory is applicable.
When both the static and dynamic disorders are present, the transport properties become complicated. Several analytic expressions and numerical simulations [31][32][33][34][35][36][37] have revealed that while any finite amount of static disorder leads to a lack of diffusion in one-or two-dimensional systems, adding a dynamic disorder can be sufficient to allow for transport to occur by destroying the phase coherence responsible for Anderson localization. In our preceding paper [38], referred to as paper I, we carried out numerically exact calculations of infinite one-dimensional disordered systems over the entire regime of dephasing (proportional to temperature and reorganization energy) within the Haken-Strol-Reineker (HSR) model [39,40]. If the dephasing rate caused by the dynamic disorder is sufficiently large such that any coherence created during the course of the time evolution is quickly destroyed, then the diffusion coefficient displays a non-monotonic dependent on temperature. When the temperature is higher than the transition value, the diffusion constant is proportional to the inverse of temperature. In the regime of weak dephasing, it is demonstrated that the diffusion constant is proportional to the temperature and the square of the Anderson localization length of the disordered system Hamiltonian. Although the results are insightful, the HSR essentially considers the dynamic disorder classically, which may lead to an incorrect transport limit. For instance, the HSR predicts zero diffusion coefficient at zero temperature, which is not consistent with the quantum prediction at very low temperatures, and results in the asymptotic equal-population of the energy levels of quantum systems, which occurs for arbitrary energy differences.
In the present work, we use the time-dependent wavepacket diffusion (TDWPD) method [41] to investigate the diffusion coefficient of the one-dimensional molecular-crystal model when both the static and dynamic disorders are present. Formally, TDWPD can be understood as an approximation to the exact stochastic Schrödinger equation (SSE) [42][43][44][45][46] or stochastic wavefunction methods [47]. The TDWPD has a similar structure to the HSR; except we start from the coherent-state representation of the phonon, and the dynamic disorders are then incorporated by stochastic complex-valued forces, which are quantitatively generated from their quantum correlation functions. The TDWPD method essentially overcomes the deficiency of HSR, i.e. it incorporates the quantum tunneling effect and approximately satisfies the detailed balance principle. Moreover, its demand on computational time is similar to that of HSR, and it can thus be applied to nanoscale systems. To compare with existing theories, we also present the results from Redfield theory [38,[48][49][50], which is a standard approach to describe the coherent motions of carriers. The diffusion coefficients from the (FGR) and Marcus formula are further illustrated to investigate the validity of the TDWPD in the strong dynamic disorder limit.
The paper is arranged as follows. In section 2, we outline the TDWPD, Redfield theory and Marcus theory used in the present work. Section 3 shows the results of numerical simulations. The concluding remarks are presented in section 4.

TDWPD method
In our previous work, we proposed a TDWPD approach that deals with the charge carrier dynamics of homogeneous systems [51]. In order to make TDWPD satisfy the detailed balance principle, we further incorporated the memory effect and obtained a non-Markovian SSE [41]. Numerical simulations have shown that the non-Markovian term begins to play a role only for the case with long memory time. In the typical organic semiconductor with static disorder, this effect can be safely neglected. Therefore, in the present work, we use the Markovian limit of the SSE. In this case, the SSE becomes similar to our original TDWPD method. Here, we still name it the TDWPD method, and the corresponding formulation is outlined as follows. The total Hamiltonian with static and dynamic disorders for charge and energy transfer is given as  with the expectation ϵ 〈 〉 nn . In a similar fashion, one can incorporate the static disorder in the electronic coupling. In the present work, we only focus on static disorders on the site energies (diagonal terms), thus the off-diagonal couplings ϵ σ′ ( ) nm are set to constants. The dynamic disorder is described by the Hamiltonian H ph of molecular vibrational motions and the electron-phonon coupling − H e ph . They are given by where N ph n is the number of phonon modes in the n-th site, + a nmj and a nmj are the creation and annihilation operators of the jth phonon mode with frequency ω nmj , respectively (the mass scaled coordinates and units with  = 1 are used throughout the paper), C nmj represents the strength of electron-phonon interactions which cause the random fluctuations of site energies (n = m) and electronic couplings ( ≠ n m) (for the abbreviation, we use the subscript j instead of nmj in the following). It is noted that C j are related to the displacements ΔQ j from equilibrium configurations by the relationship . The corresponding mode-specific reorganization energies are ω ΔQ (1/2) j j 2 2 , in which ω j and ΔQ j can be easily obtained for the realistic molecules from electronic structure calculations [12,52].
In the following, we transform the total Hamiltonian equation (1) into the interaction representation with respect to the phonon Hamiltonian. In this case, the total Hamiltonian becomes time-dependent, and is expressed by and has a formal solution , and its explicit expression is , and ψ α t ( ) represents the component of the total wavefunction on the coherent state α . With the assumption that at the initial time the phonons are in the thermal equilibrium distribution ρ ph T and the carrier is located at the ψ (0) state, the reduced density operator of the carrier can be then written as , and β denotes the coherent state of phonon distribution at the initial time. With the use of creation (annihilation) operators and coherent states, we can obtain the following time evolution equation for the carrier wavefunction . Here, it must be pointed out that the stochastic forces of forward and backward propagation are not explicitly related through one stochastic equation, and the non-Markovian term thus turns up to incorporate the correlation such that equation (11) can be used to obtain the correct carrier dynamics and is equivalent to SSEs with a functional derivative [42]. Alternatively, Moix and Cao have recently proposed a hybrid stochastic hierarchical equation of motion to treat the non-Markovian characteristics of the bath exactly [47]. For simplicity of numerical simulations, we make several approximations of equation (11) (see appendix A or [41] for more details) and can finally obtain represents the probability of the wavefunction starting from coherent state β at initial time to the α state at time t. The pure carrier dynamics can be obtained from the sum over α and β indices with the thermal weight factor. In the numerical implementation, one alteratively considers α and β as stochastic quantities, and the resulting correlation function of F(t) for the diagonal terms should satisfy where β = ( ) , T is temperature, and ω J ( ) is the spectral density function. α t ( ) in equation (12) is the zero temperature limit of C(t) and does not depend on the stochastic quantities. The similar correlation functions for off-diagonal terms can be derived if the dynamic disorders are incorporated in the electronic couplings, i.e. if C nmj are not zero for ≠ n m.
It is noted that the Markovian integral part in equation (12) is obtained from the first-order approximation [41], and the more accurate expression can be derived as the higher order terms are incorporated in the memory integral. Nevertheless, the numerical simulations have approved that equation (12) can reproduce the exact results of the rigorous quantum approaches for systems with not too small a frequency and not too large a reorganization energy at low and high temperatures, and the possible reason for this consistence is that most of the memory effect is already incorporated from stochastic forces via equation (13).
From equation (13), the stochastic dynamic disorder F(t) can be generated by and ω max is the upper cutoff frequency of the spectral density. ϕ k are the independent random phases which are uniformly distributed over the interval π (0, 2 ). Now, we can solve equation (12) at a given F(t) (i.e. at the given β and α) by equation (14). We must emphasize that coherent states are used in the derivation of the TDWPD approach but are not used in the simulations. As the carrier wavefunction is expanded in the site where d is the dimension of the system and 〈 〉 q t ( ) 2 is the mean-squared displacement of the carrier, which can be calculated according to here, l is the distance between two adjacent sites [38].

Secular Redfield equations
When the electron-phonon interaction is weak, the perturbation theory leads to the well-known Redfield equations. Under the secular approximation [48][49][50] , , Here, μν μν R , and μμ κκ R , represent the coherent dephasing and population relaxation rates, respectively, and can be expressed as where m denote the site eigenstates, the state vectors labeled with Greek characters, such as μ and ν , are the eigenstates of the pure carrier Hamiltonian, and ω = − μν μ ν E E is the eigenstate energy difference. Once the reduced density matrix ρ μν t ( ) is obtained, the meansquared displacement of the carrier can be easily calculated by where μ T n is the transform matrix element between the site and eigenstate energy level representation [38].

Marcus formula and FGR
In the limits of weak electronic coupling and strong electron-phonon interactions (large reorganization energies), the coherent motion of the carrier among sites is negligible and the dynamics of the carrier essentially obey the hopping mechanism. In this situation, the wellknown Marcus formula and FGR can be used to estimate carrier diffusion coefficients based on the assumption that the carrier only jumps to its nearest neighbor site and that successive jumps are uncorrelated. The carrier diffusion coefficient is thus given by [53,54] 2 Here, k(T) is the carrier transfer rate from the molecule n to + n 1, and all the nearest neighbor transfer rates are assumed to be the same. At the high-temperature limit, k(T) is given by the Marcus formula is the driving force. As the nuclear tunneling is incorporated, the Marcus formula in equation (22) is replaced by the following FGR expression , and λ ω = S i i i is the Huang-Rhys factor. In the practical calculations, one commonly uses the approximation of a one-effective mode [55,56]. In this case, equation (23) derived from FGR can be analytically carried out and the carrier transfer rate is given by , and I p is the modified Bessel function. This single mode method has been demonstrated to successfully predict the rigorous charge transfer rate calculated from equation (23), especially for the mode having a high effective frequency ω eff [56].

Numerical results
To model the carrier dynamics, we use a simple one-dimensional chain, although the multidimensional systems can in principle be handled by the present approaches. In the regime of parameters used in the simulations, it is found that 200 sites (N = 200) are enough to reproduce the real time dynamics, i.e. the population never arrives at the two edges of the chain during the propagation time as the initial carrier is localized at the center. For the phonon motions, the spectral density of phonons is taken as a Debye-Drude-type form where η is the electron-phonon interaction strength, the corresponding reorganization energy λ is equal to η, and ω c is the characteristic frequency of the phonon bath. Since the reorganization energy between two sites is η for the Debye-Drude-type spectral density, it allows us to investigate the effects of dynamic disorder and bath memory length independently. In the simulations, all the parameters are given in the unit of electronic coupling, i.e. ϵ 〈 〉 = + 1.0 n n , 1 . The carrier dynamics are thus determined by four parameters: σ for static disorder, ω c , η and the temperature T (β T ) for dynamic disorder. In the following, we will demonstrate the effect of dynamic disorder in the weak, intermediate and strong damping regimes on the carrier dynamics as well as its joint effect with static disorder.

Carrier dynamics in the weak dynamic disorder regime
In the weak dynamic disorder regime, the reorganization energy is much smaller than the electronic coupling, as in conventional inorganic semiconductors. In the extreme limit of the zero dynamic and static disorders, the carrier motion is ballistic and the diffusion coefficient does not exist. In the other limit with only static disorder, the diffusion coefficient becomes zero because of the Anderson localization. To reproduce the diffusive behavior, both the static and dynamic disorders are required. To show this, we have calculated the carrier dynamics and determined the diffusion coefficients as the dynamic disorder is turned on. Figure 1 plots the temperature dependence of diffusion coefficients at a fixed static disorder (σ = 1.0 2 ) for several small values of η with ω = 1.0 c . Although the dynamic disorder parameter η changes from 0.001 to 0.01, the behaviors of diffusion coefficients with respect to temperature are similar. At low temperatures, the diffusion coefficients are nearly temperature independent. With the increasing of temperature, the diffusion coefficients begin to linearly increase with temperature at first, and then decrease after they reach to maximum values at the optimal temperature. This non-monotonic property is obviously different from the band-like theory prediction that the mobilities decrease with increasing temperature, but are quite consistent with many experimental observations [5,6]. In such a weak dynamic disorder regime, Redfield theory is expected to work quite well. At low temperatures, the diffusion coefficients from Redfield theory and the TDWPD method are indeed consistent, although the results from Redfield theory become divergent at high temperatures because of the excitonic basis set in equation (17). This agreement suggests that the carrier transport is governed by the thermal hopping between delocalized quantum states, which is the signature of quantum transport as discussed in paper I [38].
The detailed analysis for the different dynamic disorders reveals that the diffusion coefficient becomes larger for a larger value of η. This can be explained by phonon-assistant carrier transport. More interestingly, the optimal temperatures become lower for larger η, but the maximum values of the coefficients are almost the same, at very high temperatures, and the diffusion coefficient is nearly proportional to / T 1 . This behavior further demonstrates that the carrier transport follows the thermal hopping-type mechanism in these high-temperature limits because it can be approximately predicted from classical Marcus formula.
To investigate the influence of static disorder on carrier diffusion coefficients, we fix the dynamic disorder parameters η to 0.002 and ω c to 1.0, and calculate the diffusion coefficients with several values of static disorder σ. The obtained results are shown in figure 2. Although the temperature dependence of the diffusion coefficients is similar to that of figure 1, i.e. there are maximum diffusion coefficients with respect to temperature, the diffusion coefficients decrease with increasing strength of static disorder, which is opposite to the property of the dynamic disorder. At very high temperatures, moreover, the values of the diffusion coefficients become consistent for three different σ 2 .
One possible reason for the above properties of diffusion coefficients may arise from the ratio of strength between the dynamic and static disorder. At low temperatures, the static disorder dominates the carrier dynamics, and the tunneling effect of dynamic disorder breaks down the Anderson localization, resulting in small and temperature-independent diffusion coefficients. With the increasing temperature, the dynamic disorder begins to play an important role. Due to the static disorder, the carrier motion follows the hopping-type mechanism; therefore, phonon-assistant carrier transport occurs. At high temperatures, dynamic disorder becomes quite strong compared to static disorder, and the carrier dynamics are mainly determined by the dynamic ones, which also explains the failure of the Redfield equation at high temperatures.

Carrier dynamics in the intermediate dynamic disorder regime
In most organic semiconductors, the reorganization energy generally has the same order of magnitude as the electronic coupling. In this case, the carrier dynamics are diffusive even without the static disorder, as shown in our previous work [41,51]. To investigate the effect of static disorder, we display in figure 3 the temperature-dependent diffusion coefficients for the given dynamic disorder parameters η = 0.5 and ω = 1.0 c , and for several different values of static disorder σ. For the zero static disorder (σ = 0 2 ), the diffusion coefficient has a band-like property, i.e. it monotonously decreases with increasing temperature, consistent with our previous investigations. The implicit mechanism has been explained by the coherent and tunneling effects of carrier motions. With the increasing static disorder, the diffusion coefficient always decreases in the complete temperature regime. As σ 2 increases to a certain value, the optimal temperature appears, and the temperature dependence of the diffusion coefficients becomes similar to that in figures 1 and 2, suggesting that the carrier dynamics begin to follow the hopping mechanism. This can be easily understood by the fact that the action of static disorder always destroys the carrier coherence, and the carrier dynamics thus change gradually from coherent to hopping with the increasing strength of the static disorder. Figure 3 also shows that the stronger static disorder leads to the higher optimal temperature. But these optimal temperatures are much smaller than those in the case of weak dynamic disorder (see figures 1 and 2). Therefore, one expects that the optimal temperature is related to both the static and dynamic disorders: the larger dynamic disorder and smaller static disorder will result in lower optimal temperature. Indeed, we have found that the transition temperatures for η = 2.0 are obviously smaller than those in figure 3 (not shown). At sufficiently high temperatures, the dynamic disorder effect should be more dominant than the static disorder one. Figure 3 clearly shows this behavior because the diffusion coefficients for different static disorders become the same and independent of σ 2 .

Carrier dynamics in the strong dynamic disorder regime
For the strong dynamic disorder, the carrier commonly follows the hopping motion because of the high-energy barrier from the large reorganization energy, except at low temperatures where the coherent motion may still be present. In this deep tunneling regime, although the FGR is accurate for the hopping rate between two sites, it neglects the coherence of carrier motion among multi-sites (super-exchange). An accurate description of carrier dynamics with coherence and deep tunneling is still a challenge. The approximation in the present TDWPD expression requires that the dynamic disorder is not very strong. However, we have numerically demonstrated that the TDWPD method may still work even when the reorganization energy is about 10 times larger than the electronic coupling (λ ϵ ≃ 〈 〉 + 10 n n , 1 ) in a spin-boson model. Therefore, the TDWPD method should be suitable for most organic semiconductors or lightharvesting systems. In the simulations, we set λ η = = 4.0 and ω = 1.0 c , and the results are shown in figure 4 with two different static disorder strengths.
It is clearly seen from figure 4 that the optimal temperature appears in the system even without static disorder (σ = 0 2 ), which is significantly different from the results in the intermediate dynamic disorder regime. Furthermore, the incorporation of static disorder has little effect on the carrier motions although smaller diffusion coefficients are obtained at not too high temperatures. We thus consider that the carrier transport is dominated by hopping motion in the strong damping regime.
For hopping dynamics and weak electronic coupling, the FGR or its high-temperature limit, the well-known Marcus formula, should work. Therefore, we also calculate carrier diffusion coefficients with the Marcus formula and FGR. In the FGR calculation, one needs to know the effective frequency of phonon motions, which heavily depends on ω c . However, the Debye spectral density predicts a divergent value for this effective frequency, although the TDWPD method can straightforwardly use different ω c . Figure 5 displays the results with different ω c from TDWPD, the FGR and Marcus formula, respectively. It is found that ω c is a crucial factor in the strong dynamic disorder regime, as it controls the transition of carrier dynamics from band-like to hopping-type. At small ω c values, the carrier follows the hopping-type motion and the optimal temperature appears. This temperature does not depend on ω c , as shown in figure 5, because the strengths of dynamic disorder remain constant. As ω c increases, the diffusion coefficient monotonically decreases with temperature. In this case, the coherence and tunneling of carrier motion at low temperatures dominate the diffusion coefficient. Although those explicit differences exist at low temperatures, the diffusion coefficients become consistent for different values of ω c at high temperatures, where the hopping motion controls the diffusion coefficients.  In comparison with the TDWPD method, the Marcus formula, which is independent of ω c , predicts a hopping-type motion as expected, because it is suitable for the high-temperature limit and neglects the tunneling effect. However, the FGR predicts a similar property for the TDWPD method: as the effective frequency increases, the coherent motion together with the tunneling becomes dominant. Nevertheless, consistent results from FGR and TDWPD are not obtained. The possible reason for this may arise from the difference between two neighboring site simulations in the FGR and the multi-site dynamics in the TDWPD method, where the quantum correlation function of stochastic forces leads to carrier coherent motion.

Conclusion
In summary, the TDWPD approach is used to investigate the transport properties in a onedimensional molecular-crystal model incorporating both static and dynamic disorders. In order to clarify the effects of these disorders on carrier dynamics, three regimes from the weak, intermediate and strong dynamic disorder strengths are analyzed in the present paper. For the purpose of comparison, numerical simulations from the Redfield equation for weak dynamic disorder and the Marcus formula for strong dynamic disorder are also shown. At the weak dynamic disorder, the carrier diffusion coefficient is nearly independent of temperature at low temperatures, and then linearly increases with temperature ( ∝ D T); however, it becomes proportional to T 1/ at very high temperatures, which can be roughly predicted from the classical Marcus formula. Although the Redfield equation is applicable in this case, the numerical simulation shows that it can only apply to the relatively low temperature regime. In the intermediate dynamic disorder regime, the band-like transport is observed with the weak static disorder. As the static disorder increases to a certain value, however, an optimal temperature appears and the temperature dependence of the diffusion coefficients becomes similar to that in the weak dynamic disorder regime. For the strong dynamic disorder, the carrier dynamics can follow the hopping-type mechanism even without static disorder, which is significantly different from the results in the intermediate dynamic disorder regime. Furthermore, it is found that the memory time of dynamic disorder is also an important factor in controlling the transition of carrier motion from band-like to hopping-type, especially in the strong dynamic disorder regime, where the partially coherent motions of the carrier can still be maintained with enough long time memory, leading to a larger diffusion coefficient than that from conventional perturbation theory. , we can then obtain a truncated version of the exact SSE or equation (11)  which is used as equation (12).