Communications: A nonperturbative quantum master equation approach to charge carrier transport in organic molecular crystals

The problem of a quantum system that couples to a complex environment is ubiquitous in different disciplines of physics.^1,2 Typical examples include atoms or molecules in photonic band gap materials,³ molecular crystals and aggregates,^4–6 and light harvesting complexes in photosynthesis.^7–9 Standard theoretical methods to investigate the quantum dynamics in these systems, such as the Redfield theory¹⁰ or Fermi’s golden rule,² apply perturbative treatment to either the system-bath coupling, or the coupling between the quantum units. As these approximate methods have limited range of applications, general methods that treat the system-environment coupling full quantum mechanically are still lacking. In this letter, we present an exact nonperturbative quantum master equation to investigate dissipative dynamics in the above mentioned systems. The new method is applied to study charge carrier transport in organic molecular crystals (OMCs) in the presence of nonlocal electron-phonon interactions, while the model and theory presented in this study are rather general and applicable to many other problems such as spectra and excitation energy transfer in various types of quantum aggregates.

An intriguing problem in charge carrier transport in OMCs is the temperature dependence of mobility at high temperatures.^11–13 At room temperature, the mobility in many OMCs decreases with rising temperature, showing a “thermal deactivated” bandlike behavior.^11,12 However, the traditional band transport model results in a too small mean free path of the charge carrier that is inconsistent with the delocalized band picture.¹⁴ At the same time, nearly temperature-independent mobility has also been observed in weak electronic coupling crystallographic directions, indicating the existence of thermal activated mobility.¹¹ An important finding in recent theoretical and computational studies is that, in many OMCs, the electronic couplings between neighboring molecules fluctuate strongly due to the intermolecular thermal motion.^15–18 This effect of dynamic disorder in the intermolecular electronic couplings, which is also termed as nonlocal electron-phonon coupling, may be an important factor to understand the above mentioned temperature dependence of intrinsic charge carrier mobility in OMCs, and is the subject of Refs. 15, 16, and 18.

The effects of nonlocal electron-phonon coupling have been investigated by a number of authors.^{15,16,18–20} In Ref. 19, Munn and Silbey proposed a method to calculate the effects of nonlocal coupling on the polaron bandwidth and mobility. More recently, using parameters obtained from first-principle calculations, contributions from both local and nonlocal electron-phonon couplings to charge carrier transport were evaluated.^15,18 However, the above methods are based on perturbative treatments of the effective electronic coupling terms in the small polaron picture, which may not be valid when the intermolecular electronic coupling is large. Troisi and co-workers^16,20 have solved the time-dependent Schrödinger equation of a semiclassical Hamiltonian to investigate the effect of nonlocal electron-phonon coupling. However, the classical treatment of intermolecular vibrational modes is only valid at high temperatures for slow vibrational motions.

In contrast to the above mentioned works, we consider a nonperturbative quantum master equation approach for a full quantum description of the electron-phonon coupling, which is based on the recently developed Liouville space hierarchical equations of motion (HEOM) approach.^21–24 To this end, we consider a one-dimensional (1D) system described by a tight-binding model that couples to a bath of harmonic oscillators. The total Hamiltonian is written as

Here, the electronic Hamiltonian $He$ is given by

where $|N+1⟩=|1⟩$ and $|m⟩$ describes a state where the charge carrier is located on the $mth$ site. The site energy $ϵm$ is assumed to be the same for all sites.

$Hph$ describes the phonon bath. It is assumed that the transfer integral between sites $m$ and $m+1$ couples independently to its own harmonic bath, such that

where $Nbm$ is the number of bath modes. The electron-phonon coupling $He-ph$ causes fluctuations in the transfer integral between neighboring sites, and is written as

where $Fm=∑j=1Nbmcmjxmj$ and $Qm=|m⟩⟨m+1|+|m+1⟩⟨m|$⁠.

The model of independent electronic coupling fluctuations described in Eq. (4) is different from that used by Troisi and co-workers,^16,20 where fluctuations of the intermolecular electronic couplings are correlated due to the shared vibrational modes. However, since the electron-phonon couplings in Refs. 16 and 20 can be written in a form similar to Eq. (4), extending the following derivations and calculations to those different models is straightforward.

The spectral density $Jm(ω)$ is defined as

$Jm(ω)$ is assumed to be the same for all sites, and the Debye–Drude-type spectral density is used,

The spectral density in Eq. (6) leads to the following time correlation function for the bath collective coordinate $Fm$ (⁠$ℏ=1$ is used throughout this paper, and $β=1/kBT$⁠):

where $γ0=γ$⁠, and $γk=2kπ/β(k≥1)$ are the Matsubara frequencies,

In the high temperature limit,

Equation (9) indicates that the parameters $η$ and $γ$ describe the amplitude and time scale of the fluctuations in the intermolecular electronic couplings, respectively.

The exact and nonperturbative HEOM for the model Hamiltonian can be written as^21–25 

where $Lρ=[He,ρ]$⁠. The subscript $n={{n10,n11,…},…,{nN0,nN1,…}}$⁠, and $nmk±$ differs from $n$ only by changing the specified $nmk$ to $nmk±1$⁠. $ρS=ρ0$ with $0={0}={{0,0,…,},…,{0,0,…}}$ is the system reduced density operator (RDO), while the others, $ρn$⁠, are auxiliary density operators (ADOs). Equation (10) consists of a set of hierarchically coupled equations for the RDO and ADOs; the level of hierarchy for each ADO is denoted as $L=∑m=1N∑knmk$⁠. While the HEOM approach is computationally more demanding than the approximate methods, recent algorithm developments have greatly improved its efficiency.^22,24,25

We now briefly discuss the connection between the exact HEOM equation and the approximate methods of stochastic Liouville equation (SLE) and generalized master equation (GME),^4,26 which have been widely applied to charge carrier (and exciton) transport. By applying the high temperature approximation [Eq. (9)] and neglecting the imaginary part of $⟨Fm(t)Fm(0)⟩$⁠, Eq. (10) reduces to a SLE with a Gaussian–Markovian bath.²² Further Markovian approximation to the bath correlation functions leads to the original form of SLE.^22,26 Truncating the HEOM at the first level $L=1$ results in a special type of time-nonlocal generalized quantum master equation (GQME) that treats the electron-phonon coupling perturbatively to the second order.^23,27 Thus the exact HEOM can be viewed as an extension of the approximate theories of SLE and GME.

Simulations of model systems were performed using the following parameters: $γ−1=130fs$ corresponding to a bath characteristic frequency of $41cm−1$ and $η=323cm−1$⁠. Estimated from Eq. (9), the root mean square deviation of the transfer integral fluctuation is $150cm−1$ at 100 K, and $260cm−1$ at 300 K. Intermolecular electronic coupling strengths up to $J=300cm−1$ were used in the simulations, which seems to be small compared to some large electronic coupling strengths $(≈1000cm−1)$ found in quantum chemistry calculations.¹⁷ However, in many OMCs, there also exists strong local electron-phonon coupling to high frequency intramolecular modes that can renormalize the electronic coupling significantly. The effective electronic coupling is given by $Jeff=Jexp[−g2coth(βω/2)]$⁠, where $ω$ is the frequency of high frequency intramolecular mode and $g$ is the coupling constant. Since we did not consider the high frequency intramolecular modes explicitly in our model Hamiltonian, such choice of nearest site electronic coupling is reasonable.

We first investigate equilibrium properties of the model systems by starting from an initial state and propagating the HEOM equation until thermal equilibrium. As all the sites are equivalent, the equilibrium populations are evenly distributed. However, the off-diagonal terms of the reduced density matrix, which describe the coherence between different sites, contain important information to characterize the extent of charge delocalization. To describe the intersite coherence, we define the quantity $gmn=ρmn/ρmmρnn$ using the matrix elements of the system RDO. Figure 1 shows $gmn$ as a function of the distance between two sites $|m−n|$⁠, for two different average electronic couplings $J=300cm−1$ and $J=50cm−1$ at various temperatures. We can see that in the case $J=300cm−1$⁠, the charge delocalization effect is significant with a charge carrier shared by several sites. For small electronic coupling $J=50cm−1$⁠, large intersite coherence only exists for the nearest neighbors. The extent of charge carrier delocalization is also seen to be more significant at low temperatures. In related areas of quantum transport, the intersite coherence was found to be of fundamental importance in electronic energy transfer in light harvesting complexes²⁸ and conjugated polymers.^29,30 It would be interesting if the intersite coherence can also be measured experimentally for charge carrier transport in organic materials.

The diffusive motion of charge carriers is then calculated. An initial state with the charge density shared equally on three central sites of the 1D chain was used and simulations were performed with the number of sites $N$ up to 50. Figure 2 shows the mean square distances $Δ⟨(Δr)2⟩(t)≡⟨(Δr(t))2⟩−⟨(Δr(0))2⟩$ as a function of time, where $r=∑m=1Nm|m⟩⟨m|$ is the position operator, and $(Δr)2≡Tr(r2ρS)−[Tr(rρS)]2$⁠. The result of dynamics without dissipation is also plotted for comparison, where $Δ⟨(Δr)2⟩(t)∝t2$ is the result of pure coherent dynamics. It can be seen clearly from Fig. 2 that two different charge carrier diffusion mechanisms occur in the case of large and small average electronic couplings, respectively. In Fig. 2(a), where $J=300cm−1$⁠, the dynamics initially follows the coherent dynamics at short time, but turns into diffusive motion $(Δ⟨(Δr)2(t)⟩∝t)$ at longer times. The overall picture is that the coherent electronic motion is dissipated by the fluctuations of intermolecular electronic coupling, which leads to diffusive motion of the charge carrier. Increasing the temperature will further destroy the intermolecular coherence and leads to smaller diffusion constants. In this case, the charge carrier diffusion mechanism is consistent with the previous findings by Troisi and co-workers,^16,20 although our result is based on the numerical exact simulation of the model Hamitonian [Eq. (1)]. In Fig. 2(b), where $J=50cm−1$⁠, the initial dynamics is much faster than the coherent dynamics without electron-phonon coupling. Also in contrast to Fig. 2(a), the diffusion constant is larger at higher temperatures. The main reason is that since the average of the intermolecular couplings is much smaller than their thermal fluctuations, the large fluctuating part of the transfer integral becomes the main driving force for charge carrier transport, and leads to thermal assisted diffusion. Comparing with the results shown in Fig. 1 indicates that thermal assisted transport exists when the charge carrier is sufficiently localized, and the extents of intersite coherence have important implications for the charge carrier transport mechanism in organic materials.

The diffusion constants were obtained by using the equation $Δ⟨(Δr)2⟩(t)≈2Dt$ at long times, and the charge carrier mobilities were obtained using the Einstein relation $μ=eD/kBT$ and are shown in Fig. 3. The distance between neighboring sites is assumed to be $d=4Å$⁠. It can be seen that for $J=300cm−1$⁠, the mobility $μ$ decreases as the temperature increases, and follows an approximate power law $μ(T)∼1/Tα$ with the index $α=1.49$ fitted from the simulation data (dashed line). This power law behavior for large average intermolecular coupling agrees with the experimental observation of bandlike transport.^11,12 While for $J=50cm−1$⁠, the thermal assisted diffusion leads to nearly temperature-independent mobility. This new finding may also provide an explanation to the experimentally observed mobility in weak electronic coupling crystallographic directions.¹¹ 

Finally, we compare the exact HEOM result with those from perturbative GQMEs to investigate their convergence with the perturbation order. As stated above, HEOM [Eq. (10)] truncated at level $L$ leads to a time-convolution GQME that treats the system-bath coupling to order $N=2L$⁠.^23,27 Figure 4 shows $Δ⟨(Δr)2⟩(t)$ as a function of time calculated with HEOM and perturbative GQMEs to the eighth order for $J=300cm−1$ and $T=300K$⁠. It can be seen that low order (second and fourth) perturbative GQMEs result in large deviations in the diffusive dynamics, and the results converge to the HEOM result with increasing $L$⁠.

In summary, the nonperturbative quantum master equation provides a full quantum treatment of the electron-phonon couplings, and is applied to study the charge carrier transport in OMCs. Both the coherent and thermal assisted transport mechanisms are characterized within a unified theory, and the intersite coherence is found to have important implications for the charge carrier transport mechanism. The new method could also find wide applications in areas such as charge transport in DNA,³¹ and exciton dynamics in conjugated polymers,^29,30 as well as many other problems mentioned in the introduction.

This work is supported by NSFC (Grant Nos. 20733006 and 20873157), the Chinese Academy of Sciences (Grant No. KJCX2.YW.H17 and the Hundred Talents Project).