THE JOURNAL OF CHEMICAL PHYSICS 139, 174109 (2013)

Quantized Hamiltonian dynamics captures the low-temperature regime of charge transport in molecular crystals Linjun Wang,1,a) Alexey V. Akimov,1,2 Liping Chen,1 and Oleg V. Prezhdo1,a) 1 2

Department of Chemistry, University of Rochester, Rochester, New York 14627, USA Chemistry Department, Brookhaven National Laboratory, Upton, New York 11973-5000, USA

(Received 30 August 2013; accepted 22 October 2013; published online 7 November 2013) The quantized Hamiltonian dynamics (QHD) theory provides a hierarchy of approximations to quantum dynamics in the Heisenberg representation. We apply the first-order QHD to study charge transport in molecular crystals and find that the obtained equations of motion coincide with the Ehrenfest theory, which is the most widely used mixed quantum-classical approach. Quantum initial conditions required for the QHD variables make the dynamics surpass Ehrenfest. Most importantly, the first-order QHD already captures the low-temperature regime of charge transport, as observed experimentally. We expect that simple extensions to higher-order QHDs can efficiently represent other quantum effects, such as phonon zero-point energy and loss of coherence in the electronic subsystem caused by phonons. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4828863] I. INTRODUCTION

Charge transport is fundamental to all organic electronic devices, including organic field-effect transistors (OFETs), light-emitting diodes (OLEDs), and photovoltaic cells (OPVs).1–4 The carrier mobility, characterizing the charge transport efficiency of organic semiconducting functional layers, strongly influences the overall device performance. The temperature dependence of mobility marks the charge transport mechanism and has been extensively studied in literature in the past decades.5–11 Two distinct pictures are usually evoked when modeling charge transport in molecular crystals, that is, band and hopping models.8, 9 The signature of each regime is conveyed by a characteristic temperature dependence of mobility. The band model applies when the electron-phonon coupling can be regarded as a perturbation. Under such circumstances, the charge is delocalized over the whole crystal, and the mobility decreases with temperature, since higher temperature results in more disorder that destroys the translational symmetry of the electronic Hamiltonian.12–15 In contrast, the hopping model works when the electronic couplings, also called transfer integrals, are treated as a perturbation. Then, the charge is self-trapped on individual molecules due to polaronic effects, and the mobility increases with temperature since there exists a barrier for charge transfer due to nuclear reorganization.16–20 Both charge carrier mobility regimes can be described analytically at the perturbation level. The analytic expressions have been applied to different classes of organic materials.12–20 In principle, the band and hopping regimes of charge transport may coexist in the same material, because the effective electron-phonon coupling strength relies on the magnitude of nuclear vibrations, which is temperature dependent. At low temperatures, nuclear coordinate fluctuations are a) Electronic addresses: [email protected] and oleg.prezhdo@

rochester.edu

0021-9606/2013/139(17)/174109/10/$30.00

rather weak, and thus can be considered a perturbation, as in the band model. On the contrary, at high temperatures, large amplitude vibrations break the electronic symmetry, the charge becomes localized, and the hopping model takes over. This phenomenon is known as the band-to-hopping transition.21–25 The transition between the two regimes is not instantaneous and covers a finite temperature range, due to involvement of phonons with different frequencies, disorder, and other effects. At extremely low temperatures, quantum rather than thermal effects govern nuclear fluctuations, leading to a special charge transport behavior. Thus, according to the nature of the nuclear vibrations coupled to the charge carriers, charge transport can be classified into three regimes, as schematically shown in Fig. 1. The band-to-hopping transitions in molecular crystals have been extensively studied in the past decades under the framework of the polaron theory,22–24, 26–30 where the polaronic couplings are generally treated in a perturbative manner. In recent years, mixed quantum-classical dynamics (MQCD) has developed into a powerful tool to study charge transport in organic materials.31–38 As a nonperturbative approach, MQCD treats nuclear vibrations classically and electron dynamics quantum mechanically. It has two major flavors: mean field (MF), also known as Ehrenfest,39–41 and surface hopping (SH).42–44 In both cases, the electronic state is driven by the classical external field created by the nuclear trajectory. The Schrödinger representation of quantum dynamics is usually used. The main difference between the Ehrenfest and SH approaches lies in the description of the electronic back-reaction – how the classical nuclei respond to the evolution of electronic degrees of freedom. The nuclear evolution in the Ehrenfest approach is governed by a single potential energy surface (PES) averaged over all states according to the electronic probability amplitudes. In the SH strategy, nuclei always move on a distinct, typically adiabatic PES, but incorporate nonadiabatic transitions, or hops, between different PESs. The probability to stay in each PES is

139, 174109-1

© 2013 AIP Publishing LLC

174109-2

Wang et al.

FIG. 1. Schematic representation of the three regimes for charge transport: quantum regime, classical perturbative regime, and regime beyond perturbation. Typical temperature dependence of charge mobility in the log-log plot and applicable temperature range are indicated. ω is the phonon frequency, J represents the intermolecualr electronic coupling, K is the force constant of nuclear vibrations, and α (β) is the local (nonlocal) electron-phonon coupling.

determined from the corresponding quantum mechanical amplitude. By treating nuclear motions classically, MQCD neglects such quantum effects as zero-point energy, tunneling and loss of coherence. The latter in particular greatly affects quantum transition rates in condensed phase systems and has to be incorporated as a semiclassical correction.45 In the Ehrenfest approach, nuclear vibrations do not induce destruction of the electronic superposition states, indicating that the coherence time is infinite. Thus, Ehrenfest should be optimal in the limit of weak electron-phonon coupling, as requested by the band transport. The SH approach implies a finite decoherence time,44–49 corresponding to a relatively large electron-phonon coupling, and thus more likely falls into the hopping regime. Using analytical results as a benchmark, it has been recently confirmed by numerical calculations that SH is able to reproduce the hopping transport, while MF gives completely wrong trends in this regime.38 A proper implementation of SH works very well in the intermediate regime, where the transport is already band-like,37 showing a significantly wider applicability range compared to Ehrenfest. In spite of the great successes of SH, the Ehrenfest technique remains the most widely used approach, due to its simple formulation. Maybe more importantly, the most actively studied molecular crystals, such as pentacene and rubrene, generally show large mobility with a band-like temperature dependence,50–52 which favors the MF description.31, 33 The Ehrenfest studies have shown that the feedback from the electron dynamics to the nuclear motion can be safely neglected when room-temperature mobility is higher than a critical value. The value is ∼4.8 cm2 /Vs and ∼0.14 cm2 /Vs for oneand two-dimensional crystals, respectively.36 Such negligible importance of the feedback, which is a second-order effect of electron-phonon coupling, implies that a weak electronphonon coupling regime and thus band-like transport can be realized up to a certain mobility. Similar observations have been made in the SH studies, where charge transport becomes band-like when carrier mobility exceeds ∼1 cm2 /Vs.37 There-

J. Chem. Phys. 139, 174109 (2013)

fore, when the charge mobility is higher than these critical values, which is the case for many organic systems, such as oligoacene crystals and graphene-based carbon materials,50–55 MF should be a reliable approach to study charge transport. Instead of the Schrödinger representation, charge transport dynamics can be also studied within the Heisenberg representation. By treating quantum mechanically not only the electron creation and annihilation operators, but also the nuclear position and momentum variables, the time derivative of the expectation value for any system observable of interest can be achieved from the Heisenberg equation. An interesting phenomenon exists in this method: The original operator becomes coupled to higher-order operators, resulting in a chain of equations. The resulting quantized Hamiltonian dynamics (QHD) approach can provide a hierarchy of approximations to quantum dynamics.56–59 The approximations are achieved through termination of the chain with a closure that expresses the expectation values of the higher-order operators in terms of products of the expectations of the lower-order operators. At this point, it would be very interesting to obtain the exact physical interpretation of different levels of QHD approximations, and to link QHD to the conventional MQCD approaches, especially the Ehrenfest theory. Besides, following the Heisenberg formulation of quantum mechanics, the initial conditions for the QHD variables, including nuclear positions and momenta, should fulfill the quantum distributions. In this way, some quantum effects of nuclear vibrations60, 61 can be naturally incorporated. This issue is particularly important when studying the temperature dependence of mobility, especially the low temperature behavior of systems with high frequency modes.18, 62 In this study, we apply the first-order QHD to charge transport in molecular crystals. Higher-order QHDs will be considered in further investigations. We derive the equations of motion for both electrons and nuclei based on model Hamiltonians with nuclear vibrations coupled to electron site energies and intermolecular electronic couplings, and compare with the Ehrenfest expressions. We implement both classical and quantum initial conditions, and reveal the importance of the nuclear quantum effects in describing the lowtemperature transport behavior. We consider variations of the parameters of the Hamiltonians in order to establish their influence on the temperature dependence of mobility. Finally, we investigate the dependence of the results on the number of vibrations per molecular site.

II. METHODOLOGY A. Local and nonlocal Hamiltonians

In order to investigate charge transport in molecular crystals, we adopt a one-dimensional model system containing N molecular sites with periodic boundary conditions. All nearest neighbors are equally spaced by a distance L. Each molecule k has one harmonic intramolecular vibrational degree of freedom with position xk and momentum pk , and the force constant K and effective mass m independent of k. Although real molecular crystals are generally assembled layer by layer and the molecules are arranged in a herringbone

174109-3

Wang et al.

J. Chem. Phys. 139, 174109 (2013)

pattern within each layer,63, 64 the chosen one-dimensional model with one effective vibrational mode is still widely used in literature.31–38 Despite its simplicity, the model reflects most characteristics of charge transport in bulk materials. Note that the methodologies described throughout this contribution can be easily extended to higher dimensions. For each molecule k, there is one associated molecular orbital |k, representing the highest occupied molecular orbital for hole transport and the lowest unoccupied molecular orbital for electron transport. At equilibrium geometry, the electronic coupling between the nearest neighbors is J, and the onsite energies of all molecules are the same and are therefore set to zero for simplicity. The electron-phonon couplings have two flavors: (1) the local electron-phonon coupling, where the onsite energy of any molecule k is linearly modulated by xk with the local electron-phonon coupling constant α and (2) the nonlocal electron-phonon coupling, where the electronic coupling between molecule k and its neighbor k + 1 is linearly modulated by xk + 1 − xk with the nonlocal electronphonon coupling constant β. We consider both Hamiltonians as follows:36   Hlocal = J (Pk,k+1 + Pk+1,k ) + αxk Pk,k k k   1 pk2 2 + Kxk + , (1) k 2 m  [J + β(xk+1 − xk )](Pk,k+1 + Pk+1,k ) Hnonlocal = k   1 pk2 2 Kxk + , (2) + k 2 m where Pk, l ≡ |kl| are the electronic projection operators. B. QHD equations of motion

Within the QHD framework, the time derivative of a ˆ of interest can be achieved from the system observable A Heisenberg equation of motion (EOM),56–59 i¯

d ˆ ˆ Hˆ ], A = [A, dt

(3)

where Hˆ is the system Hamiltonian and the brackets denote a quantum average over the initial state. In order to derive the first-order QHD (QHD-1) EOMs for charge transport, we need to regard the nuclear positions and momenta in Eqs. (1) and (2) as quantum operators, and take into account the canonical commutation relations between position and momentum, namely, [xk , pl ] = i¯δ kl . Applying the local Hamiltonian in Eqs. (1)–(3), we can get dPij  J = (Pi,j +1  + Pi,j −1  − Pi+1,j  − Pi−1,j ) dt i¯ α + (xj Pij  − xi Pij ), (4) i¯ dxi  1 = pi , dt m

(5)

dpi  = −Kxi  − αPii . dt

(6)

Using similar procedures for the nonlocal Hamiltonian in Eq. (2), one obtains dPij  J = (Pi,j +1  + Pi,j −1  − Pi+1,j  − Pi−1,j ) dt i¯ β + (xj +1 Pi,j +1  − xj Pi,j +1  i¯ + xj Pi,j −1  − xj −1 Pi,j −1 ) β (xi+1 Pi+1,j  − xi Pi+1,j  i¯ + xi Pi−1,j  − xi−1 Pi−1,j ), −

dxi  1 = pi , dt m dpi  = −Kxi  + β(Pi,i+1  + Pi+1,i  dt − Pi,i−1  − Pi−1,i ).

(7)

(8)

(9)

Hereby, we can find that Eqs. (5), (6), (8), and (9) depend only on first-order variables, while Eqs. (4) and (7) rely also on second-order correlations between nuclear positions and electronic projections. Principally, one needs to move on to include the second-order correlations xi Pkl  inside the EOMs. However, the simplest treatment is to neglect these correlations and make a closure at the first-order using the approximation, xi Pkl  ≈ xi Pkl . As a result, Eq. (4) goes to dPij  J = (Pi,j +1  + Pi,j −1  − Pi+1,j  − Pi−1,j ) dt i¯ α (10) + (xj  − xi )Pij , i¯ and Eq. (7) becomes dPij  J = (Pi,j +1  + Pi,j −1  − Pi+1,j  − Pi−1,j ) dt i¯ β + [(xj +1  − xj )Pi,j +1  i¯ + (xj  − xj −1 )Pi,j −1 ] β [(xi+1  − xi )Pi+1,j  i¯ + (xi  − xi−1 )Pi−1,j ]. −

(11)

If one compares the final EOMs for the local Hamiltonian, i.e., Eqs. (5), (6), and (10), with the Ehrenfest formulations,36 one can immediately find that the equations for QHD-1 are exactly identical to that of the Ehrenfest approach. Similar observations can be also made for the nonlocal Hamiltonian. Thus, we can conclude that the QHD reduces to Ehrenfest at the first-order in case of no higher-order correlations between electron and nuclei.65 C. Classical and quantum initial conditions

The initial settings of nuclear positions and momenta in the EOMs are critical for the charge transport dynamics. The classical initial conditions, where both nuclear positions and

174109-4

Wang et al.

J. Chem. Phys. 139, 174109 (2013)

D. Charge carrier mobility

Using classical or quantum initial conditions, the expectation values for the electronic projections can be obtained through solving the EOMs numerically with the standard fourth-order Runge-Kutta (RK4) algorithm.70 The meansquare displacement (MSD) can then be obtained through  2  MSD(t) = Pkk (t)(kL)2 − Pkk (t)kL . k k (15) In the normal diffusive regime, there is a linear relationship between MSD and the simulation time.7 Then the diffusion coefficient can be evaluated through the time derivative of the MSD as FIG. 2. Schematic representation of the nuclear coordinate and momentum variances as functions of temperature. The trends at low and high temperatures are indicated as blue and red dashed lines, respectively, to guide the eye.

momenta are randomly chosen from Boltzmann distributions, are widely used in literature when implementing the Ehrenfest dynamics.31–38 This approach works well for high temperature investigations but fails completely to describe low temperature behaviors because of the absence of quantum effect. Due to the quantum nature of QHD, one needs to make use of the quantum distribution of nuclear positions and momenta. For a harmonic oscillator, it is well-known that nuclear position follows a Gaussian distribution,66 1 2 2 ρ(x) = √ e−x /ξ . ξ π

(12)

 ¯ coth( 2k¯ω ), where kB is the Boltzmann conHere, ξ = mω BT stant,√T is the temperature, and the vibrational frequency ω = K/m. Thus the variance of the nuclear position reads   ¯ ¯ω . σx2 = √ coth 2kB T 2 Km

(13)

In Fig. 2, we show a schematic representation of the variance as a function of temperature. In the low temperature limit, Eq. (13) goes to 2√¯Km , which is independent of temperature, representing the coordinate uncertainty of ground state nuclear vibrations. And in the high temperature limit, Eq. (13) becomes kBKT , which is actually the classical nuclear coordinate distribution. Similarly, the variance for nuclear momentum is expressed as σp2

√   ¯ω ¯ Km coth . = 2 2kB T

(14)

Therefore, using Eqs. (13) and (14) as initial conditions, we should have a better description of charge transport at low temperatures by QHD-1. Note that similar quantum initial conditions based on Wigner distributions have been widely used in nonequilibrium, nonadiabatic dynamics of finite systems.67–69

D=

dMSD 1 lim . 2 t→∞ dt

(16)

Finally, the charge carrier mobility can be calculated by means of the Einstein relation,71, 72 e D. (17) μ= kB T In the log μ vs. log T plot, an inverse power-law temperature dependence, μ ∝ T−n , becomes a straight line, where the slope n ≡ −d(log μ)/d(log T). This slope varies along different lattice directions for different molecular crystals.73 Thus, we will focus on the temperature dependence of mobility, μ(T), and this slope, n(T), in the following discussions. E. Computational details

Unless otherwise specified, we adopt the following set of reference parameters: J = 300 cm−1 , α = 3500 cm−1 /Å, β = 995 cm−1 /Å, K = 14500 amu/ps2 , and m = 250 amu, based on extensive earlier studies of charge transport in molecular crystals.31, 36–38 The investigated temperatures fall into the range of 5–500 K. The intermolecular distance, L, is set to be 5 Å. The time interval is fixed to be 0.025 fs. The total simulation time (at least 1 ps) and the system size (at least 100 sites) are adjusted to ensure that equilibrium diffusion can be obtained within the chosen simulation time without biased by boundary effects. For each calculation, 2000 realizations are carried out to get a smooth time dependence of the meansquare displacement. III. RESULTS AND DISCUSSIONS A. Critical temperatures for charge transport

As the only source of disorder to the electronic Hamiltonian, the nuclear coordinate fluctuations determine entirely the temperature dependence of mobility. As shown in Fig. 2, the nuclear coordinate variance (σx2 ) at extreme low temperatures is a constant value due to the quantum effect of nuclear vibrations. Then the time-dependent electron dynamics remains the same and the diffusion coefficient is invariant of temperature. In this regime, the temperature contribution to mobility solely comes from the Einstein relation (see Eq. (17)), resulting in a μ ∝ T−1 power-law dependence (see Fig. 1). When temperature is higher than a critical value,

174109-5

Wang et al.

Tc1 = ¯ω/kB , nuclear coordinates start to follow classical distributions, and the coordinate variance can be written as a linear relationship with temperature, i.e., σx2 = kB T /K. If the magnitude of the electronic Hamiltonian √ fluctuations, that is, ασ x for the local Hamiltonian and 2βσx for the nonlocal Hamiltonian, is small in comparison to the intermolecular transfer integral, J, the electron dynamics can be dealt with the perturbation theory. In this regime, the diffusion coefficient decreases with temperature for a linear response, and then the temperature dependence of mobility follows μ ∝ T −n1 with n1 > 1. When temperature exceeds the critical value, Tc2 = J2 K/α 2 kB and Tc2 = J2 K/2β 2 kB for local and nonlocal Hamiltonians, respectively, the charge transport should be described beyond the perturbation theory. There, the charge becomes more localized, and thus the charge transport can become thermal activated or possess a weaker power-law temperature dependence of mobility. As a result, we can make a more quantitative description of the three charge transport regimes shown in Fig. 1. Note that if the classical instead of the quantum initial conditions are used, the first regime is dismissed and the temperature range, T < Tc2 , is dominant purely by the second regime. Besides, if Tc2 < Tc1 , the second regime does not apply. Analyzing the two critical temperatures, Tc1 and Tc2 , is very helpful for understanding the overall temperature dependence of mobility. B. Role of quantum initial conditions

In Fig. 3(a), we show μ(T) based on the Hamiltonian given by Eq. (1) with only local electron-phonon couplings. Using classical initial conditions, the Ehrenfest dynamics gives an almost perfect linear relationship of μ(T) in the log-

J. Chem. Phys. 139, 174109 (2013)

log plot, which implies that the carrier mobility follows an overall power-law dependence. With quantum initial conditions instead, the temperature dependence obtained by the QHD-1 remains the same above 50 K, but the mobility is strongly reduced at lower temperatures. According to the reference parameters used in this study, the vibrational frequency is calculated to be about 40 cm−1 , corresponding to a critical temperature Tc1 of 58 K. The numerical calculations here thus confirmed the interpretations discussed above: a strong change in the transport behavior happens around Tc1 because of the nuclear quantum effect. Due to the same reason, for the Hamiltonian given by Eq. (2) with only nonlocal electronphonon couplings, the role of quantum initial conditions remains basically the same (see Fig. 3(b)). The above observations can be viewed more clearly when we extract n(T) from μ(T). As shown in Figs. 3(c) and 3(d), for both the local and nonlocal Hamiltonians, the n(T) predicted by the Ehrenfest approach using classical initial conditions are more or less constant at low temperatures. This is basically the perturbation regime shown in Fig. 1. However, they experience a reduction at temperatures higher than 150 K. This is related to the transition around Tc2 , which is calculated to be 128 K for the local Hamiltonian. The decrease of n is even more significant for the nonlocal Hamiltonian due to the difference in nature of the local and nonlocal electronphonon couplings. For the local Hamiltonian, since the site energies at equilibrium geometries are zero, the fluctuation of site energies due to nuclear vibrations always result in a enhancement of electronic disorder and thus a reduction of transport efficiency. However, this is not always the case for nonlocal Hamiltonians. There, the charge transport depends strongly on the absolute values of transfer integrals. When the

FIG. 3. Temperature dependence of (a) and (b) mobility, μ, and (c) and (d) its derivative, −d(log μ)/d(log T), for the Ehrenfest dynamics with classical initial conditions and the QHD-1 dynamics with quantum initial conditions. (a) and (c) are based on the local Hamiltonian, Eq. (1), while (b) and (d) are based on the nonlocal Hamiltonian, Eq. (2). In all cases, the reference parameter values are used: J = 300 cm−1 , α = 3500 cm−1 /Å, β = 995 cm−1 /Å, K = 14500 amu/ps2 , and m = 250 amu.

174109-6

Wang et al.

J. Chem. Phys. 139, 174109 (2013)

FIG. 4. Temperature dependence of mobility for QHD-1 calculations based on the local Hamiltonian, Eq. (1). For each calculation, either (a) the electronic coupling, J, (b) the local electron-phonon coupling, α, (c) the force constant, K, or (d) the effective mass, m, is changed, while keeping the other parameters to their reference values: J = 300 cm−1 , α = 3500 cm−1 /Å, K = 14 500 amu/ps2 , and m = 250 amu.

temperature is high enough, the fluctuation of the transfer integral can be of opposite sign of the transfer integral at equilibrium geometries, and thus transport can be relatively enhanced,19, 20 resulting in a weaker temperature dependence. For the QHD-1 approach using quantum initial conditions, n(T) is close to unity at extreme low temperatures, indicating the signatures of quantum regime shown in Fig. 1. The high temperature behaviors are similar to the Ehrenfest results using classical initial conditions. C. Role of parameters in the Hamiltonian

We start with the local Hamiltonian in Eq. (1). Based on the reference parameters, we vary each fundamental parameter within the Hamiltonian one by one, namely, the intermolecular transfer integral, J, the local electron-phonon coupling, α, the force constant, K, and the effective mass, m. The temperature dependences, μ(T) and n(T), obtained by the QHD-1 approach using quantum initial conditions are shown in Figs. 4 and 5, respectively. At temperatures below 50 K, μ(T) is simply shifted when using different J and α (see Figs. 4(a) and 4(b)), and thus n(T) keeps the same (see Figs. 5(a) and 5(b)). This is because the transition between the quantum regime and the classical perturbative regime in Fig. 1 is unrelated to J and α. However, at high temperatures, we can observe a lower n for smaller J or larger α. As we mentioned above, beyond the perturbation regime, the stronger the dynamic disorder strength the weaker the temperature dependence of mobility. When reducing J or enhancing α, the critical temperature Tc2 gets smaller, resulting in a lower n. Comparatively, K and m have a much stronger impact on the temperature behavior of charge transport (see Figs. 4(c)

and 4(d)), because they together √ determine the frequency of the nuclear vibrations, ω = K/m, which plays the central role in the quantum effect. With the increase of K or decrease of m, ω is increased, then the quantum effect plays a role up to a higher critical temperature Tc1 (see Figs. 5(c) and 5(d)). The major difference between the roles of K and m happens at high temperatures. An increase of K results in a larger Tc2 , while changing m has no impact on Tc2 . Thereby, different K can end up with similar n at high temperatures, because both Tc1 and Tc2 changes with K, and these high temperature calculations fall in the same classical perturbative regime. Smaller m results in smaller n at high temperatures since the perturbation regime is compressed, and its characteristic temperature dependence cannot be fully reflected. The results with only nonlocal electron-phonon couplings (as shown in Figs. 6 and 7) are very similar to those obtained using only local electron-phonon couplings. Recent surface hopping studies have shown that the bandlike charge carrier mobility decreases with increasing nonlocal electron-phonon coupling strength in the presence of local electron-phonon couplings.37 Thereby, we expect that combined investigations with both local and nonlocal electronphonon couplings should yield similar conclusions in the present Ehrenfest-like dynamics, and thus are not carried out. In line with the observations in Fig. 3, the major difference between the Hamiltonians with local and nonlocal electronphonon couplings happens at high temperatures. A much stronger reduction in n can be found, because strong nonlocal electron-phonon couplings at high temperatures tend to increase the effective electronic couplings between neighboring molecules, counteracting the decrease of mobility due to larger disorder.19, 20

174109-7

Wang et al.

J. Chem. Phys. 139, 174109 (2013)

FIG. 5. Temperature dependence of −d(log μ)/d(log T) for the same calculations as in Figure 4.

D. Role of multiple phonons per molecule

In real molecular crystals, there exist a lot of vibrational modes with different frequencies reacting on the same molecule.19, 23 Thereby, we study further the role of multiple phonons, using the local Hamiltonian with two uncorrelated vibrations per molecule as an example. We keep J and K to be unchanged. The electron-phonon coupling strengths are re-

√ duced to 1/ 2 times the reference  value to conserve the total reorganization energy λ = i αi2 /Ki ,36 where i covers all considered vibrational modes. We choose m1 = 250 amu with the vibrational frequency as 40 cm−1 , and m2 = m1 /16 so that the frequency is 4-fold larger. As shown in Fig. 8, the temperature dependence of mobility using two phonons per molecule falls exactly between the results based on one single phonon

FIG. 6. Temperature dependence of mobility for QHD-1 calculations based on the nonlocal Hamiltonian, Eq. (2). For each calculation, either (a) the electronic coupling, J, (b) the nonlocal electron-phonon coupling, β, (c) the force constant, K, or (d) the effective mass, m, is changed, while keeping the other parameters to their reference values: J = 300 cm−1 , β = 995 cm−1 /Å, K = 14500 amu/ps2 , and m = 250 amu.

174109-8

Wang et al.

J. Chem. Phys. 139, 174109 (2013)

FIG. 7. Temperature dependence of −d(log μ)/d(log T) for the same calculations as in Figure 6.

per molecule and proper electron-phonon couplings yielding the same reorganization energy. Thus, the major charge transport mechanism under multiple phonons per molecule can be fully understood with one single vibrational mode as systematically investigated above. E. Comparisons with experiment

Due to the great difficulty in fabricating large size and ultra pure molecular crystals for time of flight measurements,52 reliable experimental data on the temperature dependence of intrinsic mobility over a wide temperature range is not widely available in literature. The work done by Karl on naphthalene

FIG. 8. Temperature dependence of mobility with either one nuclear vibration (α = 3500 cm−1 /Å, m = 250 amu or m = 15.625 amu) or two nuclear vibrations (α 1 = α 2 = 2475 cm−1 /Å, m1 = 250 amu, and m2 = 15.625 amu) coupled to each molecule. The QHD-1 approach based on the local Hamiltonian, Eq. (1), is used.

crystals73 is widely used in comparison with theories.22–24 The reference parameters adopted in the present study are based on pentacene crystal,31 which has a slightly larger intermolecular transfer integral and a smaller electron-phonon coupling strength than naphthalene crystal.74 By varying these reference parameters (i.e., Figs. 4–7) we are able to represent the major transport behavior of other related systems, including naphthalene crystal. The transition in the slope of μ(T) observed at about 30 K, arises from the quantum effect of nuclear vibrations discussed above. A power-law temperature dependence of mobility is generally observed for higher temperatures. The slope n is around 1.5 for most transport directions, with the exception of hole transport along the a-axis, exhibiting n = 2.9.73 As shown in Figs. 4–7, the value of n calculated at the QHD-1 level falls in the range between 1 and 2, agreeing with these measurements on naphthalene. Recently, it has been found that thermal expansion of molecular crystals has a strong impact on the temperature dependence of mobility.24 With the increase of temperature, the larger spacing between neighboring molecules induce smaller intermolecular transfer integrals and also larger electronphonon couplings. When these temperature-dependent structure modifications are taken into account, a larger n is expected and thus the large n in experiment can be explained. Note that acoustic phonons, which have not been considered in the present models, may also have a strong impact on the temperature dependence of mobility at low temperatures.22 Besides, we notice that there is another transition in the experimental temperature dependence of mobility along the b-axis of naphthalene crystal at around 100 K. It can be attributed to our second transition from the perturbative regime to the nonperturbative regime (see Fig. 1). In our calculations, such transition happens at around 100–200 K, depending strongly on the intermolecular transfer integrals (see Figs. 5(a) and 7(a)).

174109-9

Wang et al.

Finally, it is important to emphasize the application range of QHD-1 for charge transport. We have proved that the Ehrenfest approach is a special case of QHD, in the absence of higher-order correlations and quantum effects. Thus, similarly to Ehrenfest, QHD-1 does not work well in the hopping regime of charge transport.38 In the regime beyond perturbation shown in Fig. 1, higher-order QHD or other advanced approaches like the flexible surface hopping37 should be used to yield more accurate descriptions. As a result, QHD-1 works properly for charge transport in a relatively low temperature range, where quantum initial conditions should be considered to achieve correct description of the nuclear quantum effects. IV. CONCLUSIONS

In summary, we have derived the equations of motion for charge transport in molecular crystals using the first-order quantized Hamiltonian dynamics theory, and found that the QHD equations of motion are identical to the Ehrenfest formalism. The difference between the Ehrenfest and QHD-1 methods arises in the initial conditions. Ehrenfest uses classical initial conditions, while the initial conditions of QHD are quantum mechanical. Neglect of correlations between different degrees of freedom constitutes the main approximation of QHD-1. The use of quantum initial conditions for the nuclear position and momentum variables allowed us to reproduce the low-temperature regime of charge transport, as observed experimentally. We suggested a three-regime picture of charge transport in order to understand the temperature dependence of mobility at the QHD-1 level. We investigated the effect of all parameters in the Hamiltonian on charge transport for cases including both single and multiple phonons per molecular site. The developed description allowed us to explain the experimental results on naphthalene crystal. ACKNOWLEDGMENTS

This work is supported by the US National Science Foundation Grant No. CHE-1300118. A.V.A. was funded by the Computational Materials and Chemical Sciences Network (CMCSN) project at Brookhaven National Laboratory under Contract No. DE-AC02-98CH10886 with the U.S. Department of Energy and supported by its Division of Chemical Sciences, Geosciences & Biosciences, Office of Basic Energy Sciences. O.V.P. acknowledges financial support of the U.S. Department of Energy, Grant No. DE-SC0006527. 1 Y.

Shirota and H. Kageyama, Chem. Rev. 107, 953 (2007). R. Murphy and J. M. J. Fréchet, Chem. Rev. 107, 1066 (2007). 3 A. P. Kulkarni, C. J. Tonzola, A. Babel, and S. A. Jenekhe, Chem. Mater. 16, 4556 (2004). 4 J.-M. Nunzi, C. R. Phys. 3, 523 (2002). 5 J.-L. Brédas, D. Beljonne, V. Coropceanu, and J. Cornil, Chem. Rev. 104, 4971 (2004). 6 V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey, and J.-L. Brédas, Chem. Rev. 107, 926 (2007). 7 L. J. Wang, G. J. Nan, X. D. Yang, Q. Peng, Q. K. Li, and Z. Shuai, Chem. Soc. Rev. 39, 423 (2010). 8 F. C. Grozema and L. D. A. Siebbeles, Int. Rev. Phys. Chem. 27, 87 (2008). 9 Z. G. Shuai, L. J. Wang, and Q. K. Li, Adv. Mater. 23, 1145 (2011). 10 A. Troisi, Chem. Soc. Rev. 40, 2347 (2011). 2 A.

J. Chem. Phys. 139, 174109 (2013) 11 F.

Ortmann, F. Bechstedt, and K. Hannewald, Phys. Status Solidi B 248, 511 (2011). 12 Y. C. Cheng, R. J. Silbey, D. A. da Silva Filho, J. P. Calbert, J. Cornil, and J.-L. Brédas, J. Chem. Phys. 118, 3764 (2003). 13 M. Q. Long, L. Tang, D. Wang, L. J. Wang, and Z. G. Shuai, J. Am. Chem. Soc. 131, 17728 (2009). 14 L. Tang, M. Q. Long, D. Wang, and Z. G. Shuai, Sci. China, Ser. B: Chem. 52, 1646 (2009). 15 L. P. Chen, L. J. Wang, Z. G. Shuai, and D. Beljonne, J. Phys. Chem. Lett. 4, 2158 (2013). 16 Y. Olivier, V. Lemaur, J.-L. Brédas, and J. Cornil, J. Phys. Chem. A 110, 6356 (2006). 17 X. D. Yang, L. J. Wang, C. L. Wang, W. Long, and Z. G. Shuai, Chem. Mater. 20, 3205 (2008). 18 G. J. Nan, X. D. Yang, L. J. Wang, Z. Shuai, and Y. Zhao, Phys. Rev. B 79, 115203 (2009). 19 L. J. Wang, Q. K. Li, Z. Shuai, L. P. Chen, and Q. Shi, Phys. Chem. Chem. Phys. 12, 3309 (2010). 20 H. Geng, Q. Peng, L. J. Wang, H. J. Li, Y. Liao, Z. Y. Ma, and Z. G. Shuai, Adv. Mater. 24, 3568 (2012). 21 L. B. Schein, C. B. Duke, and A. R. McGhie, Phys. Rev. Lett. 40, 197 (1978). 22 K. Hannewald and P. A. Bobbert, Appl. Phys. Lett. 85, 1535 (2004). 23 L. J. Wang, Q. Peng, Q. K. Li, and Z. G. Shuai, J. Chem. Phys. 127, 044506 (2007). 24 L. J. Wang, Q. K. Li, and Z. G. Shuai, J. Chem. Phys. 128, 194706 (2008). 25 S. Fratini and S. Ciuchi, Phys. Rev. Lett. 103, 266601 (2009). 26 T. Holstein, Ann. Phys. (N.Y.) 8, 343 (1959). 27 H. Sumi, J. Chem. Phys. 70, 3775 (1979). 28 R. Silbey and R. W. Munn, J. Chem. Phys. 72, 2763 (1980). 29 V. M. Kenkre, J. D. Andersen, D. H. Dunlap, and C. B. Duke, Phys. Rev. Lett. 62, 1165 (1989). 30 Y.-C. Cheng and R. J. Silbey, J. Chem. Phys. 128, 114713 (2008). 31 A. Troisi and G. Orlandi, Phys. Rev. Lett. 96, 086601 (2006). 32 M. Hultell and S. Stafström, Chem. Phys. Lett. 428, 446 (2006). 33 A. Troisi, Adv. Mater. 19, 2000 (2007). 34 A. Troisi, D. L. Cheung, and D. Andrienko, Phys. Rev. Lett. 102, 116602 (2009). 35 A. Troisi and D. L. Cheung, J. Chem. Phys. 131, 014703 (2009). 36 L. J. Wang, D. Beljonne, L. P. Chen, and Q. Shi, J. Chem. Phys. 134, 244116 (2011). 37 L. J. Wang and D. Beljonne, J. Phys. Chem. Lett. 4, 1888 (2013). 38 L. J. Wang and D. Beljonne, J. Chem. Phys. 139, 064316 (2013). 39 P. Ehrenfest, Z. Phys. 45, 455 (1927). 40 F. A. Bornemann, P. Nettesheim, and C. Schütte, J. Chem. Phys. 105, 1074 (1996). 41 O. V. Prezhdo and V. V. Kisil, Phys. Rev. A 56, 162 (1997). 42 J. C. Tully and R. K. Preston, J. Chem. Phys. 55, 562 (1971). 43 J. C. Tully, J. Chem. Phys. 93, 1061 (1990). 44 H. M. Jaeger, S. Fischer, and O. V. Prezhdo, J. Chem. Phys. 137, 22A545 (2012). 45 O. V. Prezhdo and P. J. Rossky, J. Chem. Phys. 107, 5863 (1997). 46 O. V. Prezhdo, J. Chem. Phys. 111, 8366 (1999). 47 J. C. Tully, Classical and Quantum Dynamics in Condensed Phase Simulations, edited by B. Berne, G. Ciccotti, and D. F. Coker (World Scientific, Singapore, 1998). 48 B. R. Landry and J. E. Subotnik, J. Chem. Phys. 135, 191101 (2011). 49 B. R. Landry and J. E. Subotnik, J. Chem. Phys. 137, 22A513 (2012). 50 T. Minari, T. Nemoto, and S. Isoda, J. Appl. Phys. 99, 034506 (2006). 51 V. Podzorov, E. Menard, J. A. Rogers, and M. E. Gershenson, Phys. Rev. Lett. 95, 226601 (2005). 52 M. E. Gershenson, V. Podzorov, and A. F. Morpurgo, Rev. Mod. Phys. 78, 973 (2006). 53 J. Zaumseil and H. Sirringhaus, Chem. Rev. 107, 1296 (2007). 54 J. E. Anthony, Chem. Rev. 106, 5028 (2006). 55 S. V. Morozov, K. S. Novoselov, M. I. Katsnelson, F. Schedin, D. C. Elias, J. A. Jaszczak, and A. K. Geim, Phys. Rev. Lett. 100, 016602 (2008). 56 O. V. Prezhdo and Y. V. Pereverzev, J. Chem. Phys. 113, 6557 (2000). 57 D. S. Kilin, Y. V. Pereversev, and O. V. Prezhdo, J. Chem. Phys. 120, 11209 (2004). 58 O. V. Prezhdo, Theor. Chem. Acc. 116, 206 (2006). 59 A. V. Akimov and O. V. Prezhdo, J. Chem. Phys. 137, 224115 (2012). 60 A. Donoso and C. C. Martens, Phys. Rev. Lett. 87, 223202 (2001).

174109-10

Wang et al.

J. Chem. Phys. 139, 174109 (2013)

61 J. M. Riga, E. Fredj, and C. C. Martens, J. Chem. Phys. 122, 174107 (2005).

68 H.

62 V.

69 T. C. Berkelbach, D. R. Reichman, and T. E. Markland, J. Chem. Phys. 136,

Coropceanu, R. S. Sánchez-Carrera, P. Paramonov, G. M. Day, and J.-L. Brédas, J. Phys. Chem. C 113, 4679 (2009). 63 F. R. Ahmed and D. W. J. Cruickshank, Acta Crystallogr. 5, 852 (1952). 64 R. Hajlaoui, D. Fichou, G. Horowitz, B. Nessakh, M. Constant, and F. Garnier, Adv. Mater. 9, 557 (1997). 65 C. Brooksby and O. V. Prezhdo, Chem. Phys. Lett. 346, 463 (2001). 66 C. Cohen-Tannoudji, B. Diu, and F. Laloë, Quantum Mechanics (WileyInterscience, New York, 1977). 67 X. Sun, H. B. Wang, and W. H. Miller, J. Chem. Phys. 109, 7064 (1998).

B. Wang, M. Thoss, and W. H. Miller, J. Chem. Phys. 115, 2979 (2001).

034113 (2012). H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes, 2nd ed. (Cambridge University Press, Cambridge, 1992). 71 A. Einstein, Ann. Phys. 322, 549 (1905). 72 M. von Smoluchowski, Ann. Phys. 326, 756 (1906). 73 N. Karl, Synth. Met. 133–134, 649 (2003). 74 K. Hannewald, V. M. Stojanovi´ c, J. M. T. Schellekens, P. A. Bobbert, G. Kresse, and J. Hafner, Phys. Rev. B 69, 075211 (2004). 70 W.

The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp