- Email: [email protected]
Polaron dynamics in highly ordered molecular crystals
Chemical Physics Letters 428 (2006) 446–450 www.elsevier.com/locate/cplett
Polaron dynamics in highly ordered molecular crystals M. Hultell *, S. Stafstro¨m Department of Physics, Chemistry and Biology, Linko¨ping University, S-58183 Linko¨ping, Sweden Received 6 February 2006; in final form 24 May 2006 Available online 21 July 2006
Abstract From a numerical solution of the time-dependent Schro¨dinger equation and the lattice equation of motion we obtain a microscopic picture of polaron transport in highly ordered molecular crystals in the presence of an external electric field. We have chosen the pentacene single crystal as a model system, but study the transport as a function of the intermolecular interaction strength, J. We observe a smooth transition from a nonadiabatic to an adiabatic polaronic drift process over the regime 20 < J < 120 meV. For intermolecular interaction strengths above 120 meV the polaron is no longer stable and the transport becomes band like. 2006 Elsevier B.V. All rights reserved.
1. Introduction Charge transport in synthesized organic systems has been a subject of intensive research over the last three decades. In the past, conjugated polymers and charge transfer complexes were the main focus of these studies. Today, a lot of attention is also drawn to highly ordered molecular crystals, of particular interest for their high mobilities and their potential use in molecular electronic devices [1]. Recent observations of extremely high mobilities in rubrene [2] and several reports of high mobility in pentacene [3–5] showed that it is indeed possible to compete with inorganic semiconductors in this aspect. However, the high mobility state can easily by destroyed in the presence of disorder, which reduces the electronic coupling strength between neighboring molecules, induces localization of the charge carriers and/or creates traps for the charge carriers. These effects lead to hopping transport and low mobility [6]. Thus, a high degree of structural order as well as low impurity concentration is absolutely essential in order to achieve these high mobilities.
*
Corresponding author. Fax: +46 13 137568. E-mail addresses: [email protected] (M. Hultell), [email protected] (S. Stafstro¨m). 0009-2614/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.07.042
In addition to extrinsic effects there are also intrinsic properties that are important in the context of transport. The electron–phonon coupling present in most carbon based systems can give rise to self-localization of the charge carrier, i.e., polaron formation [7,8]. The constituent molecules undergo relaxation of atomic positions in order to lower the lattice energy by an amount El, defined here as the difference in lattice energy between the charged and the neutral geometrical ground states of a single molecule. In molecular crystals, this effect is competing with the electronic intermolecular interaction J, which favors charge delocalization. Depending on the relative strengths of these two quantities, we can describe intrinsic transport in three different categories [9]. In the limit of weak intermolecular interactions (J El) the small Holstein polaron localized to a single molecule is stable [7]. Since the electronic state is localized, transport occurs via hopping from one molecule to the next. This is a nonadiabatic process for which the carrier motion is slow compared to the molecular relaxation times and the mobility depends quadratically on the intermolecular interaction strength [7,10]. For intermediate intermolecular interactions (J El) the polaron is delocalized over several molecular units and the electron probability density can sample an even larger region in space. Under the influence of an external electric field the charge density associated with the polaron moves
M. Hultell, S. Stafstro¨m / Chemical Physics Letters 428 (2006) 446–450
along the direction of the electric field and changes the molecular geometry as the carrier drifts along the system [11]. Thus, in this regime, in the absence of disorder, the activation energy barrier is small and the mobility is weakly dependent on the intermolecular interaction strength. Finally, for J El the lattice relaxation can be neglected and the transport can be described by textbook band theory [12,13]. To the best of our knowledge, no single theory [7–11] of today can be applied to the full range of relations between J and El. In particular both the small polaron model by Holstein [11] and the Marcus theory [10] are based on perturbation theory and hence are valid only in the limit of weak intermolecular interaction. In this Letter, we introduce a new approach, capable of covering the intermediate range of intermolecular interactions, within which a quantum mechanical description of the charge carrier dynamics is readily derived from the simultaneous numerical solution to the time-dependent Schro¨dinger equation _ b el jWðtÞi; i hjWðtÞi ¼H
ð1Þ
ð2Þ
Here, a viscous damping term has been introduced in the lattice equation of motion in order to simulate thermal contact between the system and the surrounding medium. Since both J and El will be uniquely defined by the Hamilb , the relationship between J and El, and possibly tonian H also the external field strength, will determine if these solutions result in band transport or in the formation of a polaron, and if the polaron has an adiabatic or nonadiabatic dynamical behavior. 2. Methodology b of the system we invoke In deriving the Hamiltonian H the r–p electron separability and treat only the p-electrons associated with the carbon 2pz orbitals quantum mechanically, assuming a tight-binding model X bp ¼ H tnn0 ðtÞð^cyn^cn0 þ ^cyn0 ^cn Þ; ð3Þ nn0
with hopping integrals ( t0 aðrnn0 ðtÞ r0 Þ intra-molecular; tnn0 ðtÞ ¼ ð0Þ inter-molecular; kS nn0
low for the contribution of an external electric field E(t), such that the total system Hamiltonian reads X b latt b ¼H b p þ jej rn EðtÞ ^cy ^cn 1 þ H H n
X
n
b latt ; ^cyn ^hnn0 ^cn0 þ H
ð5Þ
nn0 2 where EðtÞ ¼ E0 exp ½ðt tc Þ =t2w ^e for t < tc, and E0^e otherwise (for the definition of ^e, see below). Note that tc = 75 fs and tw = 25 fs for the simulations in Section 3. The exact equation of motion for the monomer displacement may now be derived by differentiating the total energy of the system with respect to the atomic coordinates rn, unraveling the interdependence with Eq. (1) through the density matrix elements qnn0 ðtÞ. If we then make the ansatz P qnn0 ðtÞ ¼ p wnp ðtÞfp wn0 p ðtÞ, where p is the molecular orbital (MO) index, wnp(t) the time-dependent MO, and fp 2 [0,1,2] the time-independent occupation number of the pth MO, {wnp(t)} will be solutions to the time-dependent Schro¨dinger equation X ihw_ np ðtÞ ¼ hnn0 ðtÞwn0 p ðtÞ: ð6Þ n0
and the lattice equation of motion b jWðtÞi k_rn : M€rn ¼ rrn hWðtÞj H
447
ð4Þ
a being the electron-lattice displacement coupling constant ð0Þ and S nn0 the initial state overlap integrals between 2pz Slater-type atomic orbitals on site n and n 0 on adjacent molecules [14,15]. For the small anticipated lattice distortions we adopt a classical description and expand the r bonding energy to second order around the undimerized state. Supplemented with the constraint of fixed molecular b latt . Finally we allength [16] we denote this contribution H
Expanding the time-dependent MO’s wnp in a basis of instantaneous eigenfunctions [17] X wnp ðtÞ ¼ unp0 ðtÞap0 p ðtÞ; ð7Þ p0
P defined by n0 hnn0 ðtÞun0 p ðtÞ ¼ unp ðtÞp ðtÞ, we obtain, additionally, the time-dependent occupation number of the eigenstate as X fp0 japp0 ðtÞj2 : ð8Þ np ðtÞ ¼ p0
In ‘adiabatic dynamics’ the system stays mostly on the same Born–Oppenheimer state, as the occupation numbers are constant or change slowly [18]. This is not sufficient if the states come close to each other in energy as will be the case in our approach. We therefore rather solve the coupled differential equations (2) and (6) simultaneously, allowing for a time-dependent occupation of instantaneous eigenstates. In what follows we will use the time evolution of np(t) as a signature of the adiabacity at hand and hold the intermolecular overlap integrals constant in time and uniform in space. Furthermore, we will treat k in Eq. (4) as a parameter, providing us with the means to scale the value of the total electronic interaction strength J, here defined as half the energy level splitting of the lowest unoccupied molecular orbital (LUMO) of a molecular dimer. We thus have a straightforward procedure to investigate the transport dynamics in general and the adiabacity in particular as a function of J. Using the pentacene single crystal as a model system for molecular crystals, we aim to obtain a basic understanding of the transition from nonadiabatic to adiabatic transport. This is the main objective in this work.
448
M. Hultell, S. Stafstro¨m / Chemical Physics Letters 428 (2006) 446–450
3. Results As a starting point for our dynamics simulations we first determined the ground state of a neutral pentacene molecule for a large number of Hamiltonian parameter sets, and compared our results against those obtained from first principles density functional calculations. The best overall fit (as detailed in Table 1) displays an error margin in the average molecular bond-length deviations of less than ˚ compared to the geometry obtained from the ab ini0.02 A tio calculations. The localization energy is El = 97 meV for this parameter set. One may then readily construct any arbitrarily large single crystal pentacene structure using available experimental data on the triclinic unit cell structure [19]. Since these crystals are highly anisotropic with the most important propagation routes found along the nearest-neighbor directions, we use in our simulations an array ofi 30 pentacene molecules stacked along ^e ¼ h 1 ^a; 1 ^ b; 0 , with ^ a and ^ b defined in Ref. [19]. The system 2
2
is finite and no periodic boundary conditions are applied. Having derived the initial state geometry {rn(0)} of the system, we proceed to simulate charge transport by introducing an additional electron into the LUMO level p = 331 (such that fp = 2 for p 6 330 and fp = 1 for p = 331) of the system at rest ð_rn ¼ 0Þ, together with a momentary initial potential well of 50 meV placed on the first three molecules in the array, the presence of which will localize the charge to this end of the system. At the first timestep the potential well is removed and the accompanying lattice displacements start to develop. The transient response to the localization of charge to the end of the system is a strong distortion of molecular geometries in this region. The distortions decay however quite rapidly during the onset of the electric field (t < tc = 75 fs) and the propagating charge in a system where J = 50 meV and E0 = 5.0 · 104 V/cm finds its steady state localization after 170 fs. In Fig. 1 are shown the consecutive curves for the time evolution of: (a) the total bond-length deviations per molecule relative to the initial state DR and (b) the net charge per molecule, DQ, for molecule . . .5, 6, 7, . . . along the array of 30 pentacene molecules. At each instant of time, roughly five molecules are involved in the coupled charge-lattice Table 1 Hamiltonian parameter set Parameter
Unit
Reference hopping integral: t0 = 2.66 Reference bond-length: r0 = 1.412 Reference bond angle: h0 = 2p/3 Reference torsion angel: /0 = 0 Electron–phonon coupling constant: a = 6.8 Harmonic distance spring constant: K1 = 72.0 Harmonic bond angle spring constant: K2 = 70.0 Harmonic torsion angle spring constant: K3 = 200.0 Damping constant: k = 1.0 · 105 Mass of a CH-group: M = 1349312601.0
[eV] ˚] [A [rad] [rad] ˚] [eV/A ˚ 2] [eV/A [eV/rad2] [eV/rad2] ˚ 2] [eV as/A 2 ˚2 [eV(as) /A ]
distortions, indicating a localized charge carrier. Since there is an almost instantaneous response in DR towards changes in DQ, the charge carrier is clearly of polaronic type. The total time during which the charge carrier affects a single pentacene molecule is above 100 fs. During this time the center of charge has moved more than 5 molecular distances. From this behavior we conclude that the motion of the polaron is best described as a drift process rather than hopping. In this situation the localization energy plays a less important role than in the case of polaron hopping [10,11]. Since the electron density associated with the polaron grows gradually on the molecules in the forward direction of polaron motion there is in fact no barrier for this type of transport. The situation is similar to that of polaron motion along a polymer chain [20]. The lack of an energy barrier also implies that the process is independent on the temperature of the system, as will be discussed further below. A detailed study of the time elapsed as the center of the polaron moves from one molecule to the next through the system reveals that after a short acceleration phase, the polaron continues to drift through the system at virtually constant velocity v. The behaviour of v as a function of J for E0 = 5.0 · 104 V/cm is depicted in Fig. 2. It should be stressed that computations become increasingly demanding in the lower region of intermolecular interactions where the time of acceleration grows rapidly with decrements in J. We therefore adopt a slightly different approach in the dynamic simulations for J 6 40 meV and reduce the size of the systems to twenty molecules. The value of J is initialized at 50 meV and the charge carrier is allowed to propagate through the system for 100 fs. Then J is set to the desired value of 20, 30, and 40 meV, respectively. This procedure results in a substantial reduction in CPU time for the simulations. In the limit of weak intermolecular interactions we find charge carriers to be immobile for J 6 20 meV. We also observe an upper boundary of v(J) above which the width of the polaron becomes so large that the concept of a localized charge carrier breaks down. In this regime, i.e. above J 120 meV, charge carrier propagation is best described in terms of band transport. This is in agreement with the discussion above since these values of J exceed the localization energy (El = 97 meV). Finally, in the intermediate regime we find that v(J), although strictly increasing, is a function that is concave down on the interval of 20 < J < 120 meV (see Fig. 2). This regime most probably spans the value of J relevant for pentacene crystals [13]. We therefore conclude that the mobilities between 22 cm2/V s at J = 30 meV and 130 cm2/V s at J = 100 meV (for E0 = 5.0 · 104 V/cm) set an upper limit on the polaron mobility in pentacene single crystals, since the simulations are performed for a defect free system. Since we are dealing with systems at low temperatures it could be argued that a discretization of phonons imposed by a quantum mechanical description of the lattice is crucial for an accurate description of charge transport.
M. Hultell, S. Stafstro¨m / Chemical Physics Letters 428 (2006) 446–450
DR [A]
0.06
449
a
0.04
5
6
7
8
9
10
11
12
13
5
6
7
8
9
10
11
12
13
0.02 0
b
DQ [qe]
0.3 0.2 0.1 0 100
150
200
250
300
Time [fs] Fig. 1. The time evolution of: (a) total molecular bond-length deviation (DR) from the initial state geometry and (b) net charge (DQ) per molecule for J = 50 meV and E0 = 5.0 · 104 V/cm. Only the curves of molecules 5–13 are highlighted. 0.7 0.6
v [Å/fs]
0.5 0.4 0.3 0.2 0.1 0 0
20
40
60
80
100
120
J [meV]
Fig. 2. The dependence of J on carrier velocity v for applied field strength E0 = 5.0 · 104 V/cm.
This is true in the nonadiabatic limit of weak intermolecular interactions and may explain the observed immobility of charge carriers below J 6 20 meV. However, increasing the lattice kinetic energy using nonzero initial atomic velocities ð_rn Þ so as to simulate a temperature T ¼ hM r_ 2n i=3k B ¼ 100 K only produced negligible corrections to the functional dependence of v(J) on the interval 20 < J < 120 meV1. This behaviour is consistent with the drift process with very low activation energy barriers for polaronic motion previously noted and indicate that a classical description of the lattice is indeed sufficient for the previously derived results to be reliable. In this context though it is important to point out that J is kept constant and uniform throughout our simulations. This treatment is used to steady the intrinsic nature of polaron transport in ordered molecular crystals. In reality, 1
Note that in this case the value of k is set to zero in order to preserve the temperature by avoiding heat to dissipate from the system. Comparative simulations show that this will not effect the transport dynamics.
however, molecular motion modulation of the intermolecular transfer integrals results in a power law temperature dependence of the mobility [21,22]. In the future we aim to include also time-dependent intermolecular interactions in our model to study the transport dynamics of heated systems. Despite the described situation of a polaron drift rather than polaron hopping depicted in Fig. 1, transport in the region of low J can be described as nonadiabatic, i.e., the dynamics requires multiple electronic states. We can obtain a very direct view of the nonadiabatic behavior of the polaron dynamics by studying the time evolution of the occupation number of the instantaneous eigenstates (see Eq. 8 above). Fig. 3a shows the occupation numbers for J = 50 meV and E0 = 5.0 · 104 V/cm. The occupancy after the field has reached its maximum value E0 is in the higher region of the unoccupied spectrum of the manifold of pentacene LUMO levels. Due to the weak intermolecular interaction relative to the localization energy, the electronic states are localized. The states that localize to the high (low) potential region of the system are Stark shifted upwards (downwards) in the spectrum. The electron is injected into the high potential energy region and will consequently occupy the higher lying states. Thus, the electron is initially in an excited state in the instantaneous picture (e.g. at t = 100 as in Fig. 3a). The motion towards the low potential region has a driving force in terms of the change in the occupancy of the instantaneous eigenstates towards lower energy as a function of time. This is clearly seen from Fig. 3a. The process is driven by the external electric field and is essentially independent on the lattice vibrations since transport occurs as a result of polaron drift rather than polaron hopping. The lattice kinetic energy is instead generated as a result of the geometrical response of the molecule to the presence of electric charge. Fig. 3b shows the occupancy of the instantaneous eigenstates in the case of J = 100 meV and E0 = 5.0 · 104 V/cm. Clearly, the behavior is different from that shown in Fig. 3a.
450
M. Hultell, S. Stafstro¨m / Chemical Physics Letters 428 (2006) 446–450
strength J and localization energy El. Using the pentacene single crystal as a model system, we show that in the regime of intermediate J (i.e., J El) there is a smooth transition from a nonadiabatic to an adiabatic polaronic drift process on the region of 20 < J < 120 meV. For intermolecular interaction strength above 120 meV the polaron is no longer stable and the transport becomes band like. We thus cover a regime between J El, where text book band theory applies, and J El, where the polaron localize to a single molecule and transport is due to a nonadiabatic hopping process. Acknowledgement Financial support from the Center of Organic Electronics (COE), Swedish Foundation of Strategic Research, is gratefully acknowledged. References
Fig. 3. The time evolution of the occupation number at E = 5.0 · 104 V/ cm, (a) J = 50 meV and (b) J = 100 meV.
The occupancy is now almost constant over time indicating that the transport is adiabatic [18]. Still, however, we are dealing with a localized charge carrier. The difference as compared to the low J case is that all other states are delocalized and therefore do not exhibit such strong Stark shift as in the previous case. The only state which is really sensing the external field is the polaron state which itself drifts along the system towards the low potential energy regime. Thus, the dynamics is in this case best described as a change in the spatial location of one and the same occupied state as opposed to the nonadiabatic case in which the electron performs a sequence of hopping events between states that are localized to different regions in space. 4. Conclusions In summary, we have introduced an approach to study the transport of charge carriers in perfectly ordered molecular crystals as a function of intermolecular interaction
[1] P. Batail (Ed.), Molecular Conductors of Chemical Reviews, vol. 104, ACS Publications, 2004. [2] V. Podzorov, E. Menard, A. Borissov, V. Kiryukhin, J.A. Rogers, M.E. Gershenson, Phys. Rev. Lett. 93 (8) (2004) 086602. [3] W. Warta, L.B. Stehle, N. Karl, Appl. Phys. A 36 (1985) 163. [4] H. Klauk, M. Halik, U. Zschieschang, G. Schmid, W. Radlik, W. Weber, J. Appl. Phys. 92 (2002) 5259. [5] O.D. Jurchescu, T. Baas, T.T.M. Palstra, Appl. Phys. Lett. 84 (2004) 3061. [6] A. Miller, E. Abrahams, Phys. Rev. 120 (1960) 745. [7] T. Holstein, Ann. Phys. 8 (1959) 706. [8] D. Emin, Phys. Rev. B 61 (21) (2000) 14543. [9] V. May, O. Kuhn, Charge and Energy Transfer Dynamics in Molecular Systems, Wiley-VCH, 2004, Chapter 6. [10] R.A. Marcus, Rev. Mod. Phys. 65 (1993) 599. [11] T. Holstein, Ann. Phys. 8 (1959) 725. [12] J.M. Ziman, Electrons and Phonons, Clarendon, Oxford, 1982. [13] Y.C. Cheng, R.J. Silbey, D.A. da Silva Filho, J.P. Calbert, J. Cornil, J.L. Bre´das, J. Chem. Phys. 118 (8) (2003) 3764. [14] R.S. Mulliken, C.A. Reike, D. Orloff, H. Orloff, J. Chem. Phys. 17 (1949) 1248. [15] A. Hansson, S. Stafstro¨m, Phys. Rev. B 67 (2003) 075406. ˚ . Johansson, S. Stafstro¨m, Phys. Rev. Lett. 86 (2001) 3602. [16] A [17] S. Block, H.W. Streitwolf, J. Phys.; Condens. Matter 8 (1996) 889. [18] H. Nakamura, Nonadiabatic Transition: Concepts, Basic Theories and Applications, World Scientific, Singapore, 2002. [19] R.B. Campbell, J. Monteath Robertson, J. Trotter, Acta Cryst. 14 (1961) 705. [20] A.A. Johansson, S. Stafstro¨m, Phys. Rev. B 69 (2004) 235205. [21] K. Hannewald, P.A. Bobbert, Appl. Phys. Lett. 85 (2004) 1535. [22] A. Troisi, G. Orlandi, Phys. Rev. Lett. 96 (2006) 086601.
Polaron dynamics in highly ordered molecular crystals M. Hultell *, S. Stafstro¨m Department of Physics, Chemistry and Biology, Linko¨ping University, S-58183 Linko¨ping, Sweden Received 6 February 2006; in final form 24 May 2006 Available online 21 July 2006
Abstract From a numerical solution of the time-dependent Schro¨dinger equation and the lattice equation of motion we obtain a microscopic picture of polaron transport in highly ordered molecular crystals in the presence of an external electric field. We have chosen the pentacene single crystal as a model system, but study the transport as a function of the intermolecular interaction strength, J. We observe a smooth transition from a nonadiabatic to an adiabatic polaronic drift process over the regime 20 < J < 120 meV. For intermolecular interaction strengths above 120 meV the polaron is no longer stable and the transport becomes band like. 2006 Elsevier B.V. All rights reserved.
1. Introduction Charge transport in synthesized organic systems has been a subject of intensive research over the last three decades. In the past, conjugated polymers and charge transfer complexes were the main focus of these studies. Today, a lot of attention is also drawn to highly ordered molecular crystals, of particular interest for their high mobilities and their potential use in molecular electronic devices [1]. Recent observations of extremely high mobilities in rubrene [2] and several reports of high mobility in pentacene [3–5] showed that it is indeed possible to compete with inorganic semiconductors in this aspect. However, the high mobility state can easily by destroyed in the presence of disorder, which reduces the electronic coupling strength between neighboring molecules, induces localization of the charge carriers and/or creates traps for the charge carriers. These effects lead to hopping transport and low mobility [6]. Thus, a high degree of structural order as well as low impurity concentration is absolutely essential in order to achieve these high mobilities.
*
Corresponding author. Fax: +46 13 137568. E-mail addresses: [email protected] (M. Hultell), [email protected] (S. Stafstro¨m). 0009-2614/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.07.042
In addition to extrinsic effects there are also intrinsic properties that are important in the context of transport. The electron–phonon coupling present in most carbon based systems can give rise to self-localization of the charge carrier, i.e., polaron formation [7,8]. The constituent molecules undergo relaxation of atomic positions in order to lower the lattice energy by an amount El, defined here as the difference in lattice energy between the charged and the neutral geometrical ground states of a single molecule. In molecular crystals, this effect is competing with the electronic intermolecular interaction J, which favors charge delocalization. Depending on the relative strengths of these two quantities, we can describe intrinsic transport in three different categories [9]. In the limit of weak intermolecular interactions (J El) the small Holstein polaron localized to a single molecule is stable [7]. Since the electronic state is localized, transport occurs via hopping from one molecule to the next. This is a nonadiabatic process for which the carrier motion is slow compared to the molecular relaxation times and the mobility depends quadratically on the intermolecular interaction strength [7,10]. For intermediate intermolecular interactions (J El) the polaron is delocalized over several molecular units and the electron probability density can sample an even larger region in space. Under the influence of an external electric field the charge density associated with the polaron moves
M. Hultell, S. Stafstro¨m / Chemical Physics Letters 428 (2006) 446–450
along the direction of the electric field and changes the molecular geometry as the carrier drifts along the system [11]. Thus, in this regime, in the absence of disorder, the activation energy barrier is small and the mobility is weakly dependent on the intermolecular interaction strength. Finally, for J El the lattice relaxation can be neglected and the transport can be described by textbook band theory [12,13]. To the best of our knowledge, no single theory [7–11] of today can be applied to the full range of relations between J and El. In particular both the small polaron model by Holstein [11] and the Marcus theory [10] are based on perturbation theory and hence are valid only in the limit of weak intermolecular interaction. In this Letter, we introduce a new approach, capable of covering the intermediate range of intermolecular interactions, within which a quantum mechanical description of the charge carrier dynamics is readily derived from the simultaneous numerical solution to the time-dependent Schro¨dinger equation _ b el jWðtÞi; i hjWðtÞi ¼H
ð1Þ
ð2Þ
Here, a viscous damping term has been introduced in the lattice equation of motion in order to simulate thermal contact between the system and the surrounding medium. Since both J and El will be uniquely defined by the Hamilb , the relationship between J and El, and possibly tonian H also the external field strength, will determine if these solutions result in band transport or in the formation of a polaron, and if the polaron has an adiabatic or nonadiabatic dynamical behavior. 2. Methodology b of the system we invoke In deriving the Hamiltonian H the r–p electron separability and treat only the p-electrons associated with the carbon 2pz orbitals quantum mechanically, assuming a tight-binding model X bp ¼ H tnn0 ðtÞð^cyn^cn0 þ ^cyn0 ^cn Þ; ð3Þ nn0
with hopping integrals ( t0 aðrnn0 ðtÞ r0 Þ intra-molecular; tnn0 ðtÞ ¼ ð0Þ inter-molecular; kS nn0
low for the contribution of an external electric field E(t), such that the total system Hamiltonian reads X b latt b ¼H b p þ jej rn EðtÞ ^cy ^cn 1 þ H H n
X
n
b latt ; ^cyn ^hnn0 ^cn0 þ H
ð5Þ
nn0 2 where EðtÞ ¼ E0 exp ½ðt tc Þ =t2w ^e for t < tc, and E0^e otherwise (for the definition of ^e, see below). Note that tc = 75 fs and tw = 25 fs for the simulations in Section 3. The exact equation of motion for the monomer displacement may now be derived by differentiating the total energy of the system with respect to the atomic coordinates rn, unraveling the interdependence with Eq. (1) through the density matrix elements qnn0 ðtÞ. If we then make the ansatz P qnn0 ðtÞ ¼ p wnp ðtÞfp wn0 p ðtÞ, where p is the molecular orbital (MO) index, wnp(t) the time-dependent MO, and fp 2 [0,1,2] the time-independent occupation number of the pth MO, {wnp(t)} will be solutions to the time-dependent Schro¨dinger equation X ihw_ np ðtÞ ¼ hnn0 ðtÞwn0 p ðtÞ: ð6Þ n0
and the lattice equation of motion b jWðtÞi k_rn : M€rn ¼ rrn hWðtÞj H
447
ð4Þ
a being the electron-lattice displacement coupling constant ð0Þ and S nn0 the initial state overlap integrals between 2pz Slater-type atomic orbitals on site n and n 0 on adjacent molecules [14,15]. For the small anticipated lattice distortions we adopt a classical description and expand the r bonding energy to second order around the undimerized state. Supplemented with the constraint of fixed molecular b latt . Finally we allength [16] we denote this contribution H
Expanding the time-dependent MO’s wnp in a basis of instantaneous eigenfunctions [17] X wnp ðtÞ ¼ unp0 ðtÞap0 p ðtÞ; ð7Þ p0
P defined by n0 hnn0 ðtÞun0 p ðtÞ ¼ unp ðtÞp ðtÞ, we obtain, additionally, the time-dependent occupation number of the eigenstate as X fp0 japp0 ðtÞj2 : ð8Þ np ðtÞ ¼ p0
In ‘adiabatic dynamics’ the system stays mostly on the same Born–Oppenheimer state, as the occupation numbers are constant or change slowly [18]. This is not sufficient if the states come close to each other in energy as will be the case in our approach. We therefore rather solve the coupled differential equations (2) and (6) simultaneously, allowing for a time-dependent occupation of instantaneous eigenstates. In what follows we will use the time evolution of np(t) as a signature of the adiabacity at hand and hold the intermolecular overlap integrals constant in time and uniform in space. Furthermore, we will treat k in Eq. (4) as a parameter, providing us with the means to scale the value of the total electronic interaction strength J, here defined as half the energy level splitting of the lowest unoccupied molecular orbital (LUMO) of a molecular dimer. We thus have a straightforward procedure to investigate the transport dynamics in general and the adiabacity in particular as a function of J. Using the pentacene single crystal as a model system for molecular crystals, we aim to obtain a basic understanding of the transition from nonadiabatic to adiabatic transport. This is the main objective in this work.
448
M. Hultell, S. Stafstro¨m / Chemical Physics Letters 428 (2006) 446–450
3. Results As a starting point for our dynamics simulations we first determined the ground state of a neutral pentacene molecule for a large number of Hamiltonian parameter sets, and compared our results against those obtained from first principles density functional calculations. The best overall fit (as detailed in Table 1) displays an error margin in the average molecular bond-length deviations of less than ˚ compared to the geometry obtained from the ab ini0.02 A tio calculations. The localization energy is El = 97 meV for this parameter set. One may then readily construct any arbitrarily large single crystal pentacene structure using available experimental data on the triclinic unit cell structure [19]. Since these crystals are highly anisotropic with the most important propagation routes found along the nearest-neighbor directions, we use in our simulations an array ofi 30 pentacene molecules stacked along ^e ¼ h 1 ^a; 1 ^ b; 0 , with ^ a and ^ b defined in Ref. [19]. The system 2
2
is finite and no periodic boundary conditions are applied. Having derived the initial state geometry {rn(0)} of the system, we proceed to simulate charge transport by introducing an additional electron into the LUMO level p = 331 (such that fp = 2 for p 6 330 and fp = 1 for p = 331) of the system at rest ð_rn ¼ 0Þ, together with a momentary initial potential well of 50 meV placed on the first three molecules in the array, the presence of which will localize the charge to this end of the system. At the first timestep the potential well is removed and the accompanying lattice displacements start to develop. The transient response to the localization of charge to the end of the system is a strong distortion of molecular geometries in this region. The distortions decay however quite rapidly during the onset of the electric field (t < tc = 75 fs) and the propagating charge in a system where J = 50 meV and E0 = 5.0 · 104 V/cm finds its steady state localization after 170 fs. In Fig. 1 are shown the consecutive curves for the time evolution of: (a) the total bond-length deviations per molecule relative to the initial state DR and (b) the net charge per molecule, DQ, for molecule . . .5, 6, 7, . . . along the array of 30 pentacene molecules. At each instant of time, roughly five molecules are involved in the coupled charge-lattice Table 1 Hamiltonian parameter set Parameter
Unit
Reference hopping integral: t0 = 2.66 Reference bond-length: r0 = 1.412 Reference bond angle: h0 = 2p/3 Reference torsion angel: /0 = 0 Electron–phonon coupling constant: a = 6.8 Harmonic distance spring constant: K1 = 72.0 Harmonic bond angle spring constant: K2 = 70.0 Harmonic torsion angle spring constant: K3 = 200.0 Damping constant: k = 1.0 · 105 Mass of a CH-group: M = 1349312601.0
[eV] ˚] [A [rad] [rad] ˚] [eV/A ˚ 2] [eV/A [eV/rad2] [eV/rad2] ˚ 2] [eV as/A 2 ˚2 [eV(as) /A ]
distortions, indicating a localized charge carrier. Since there is an almost instantaneous response in DR towards changes in DQ, the charge carrier is clearly of polaronic type. The total time during which the charge carrier affects a single pentacene molecule is above 100 fs. During this time the center of charge has moved more than 5 molecular distances. From this behavior we conclude that the motion of the polaron is best described as a drift process rather than hopping. In this situation the localization energy plays a less important role than in the case of polaron hopping [10,11]. Since the electron density associated with the polaron grows gradually on the molecules in the forward direction of polaron motion there is in fact no barrier for this type of transport. The situation is similar to that of polaron motion along a polymer chain [20]. The lack of an energy barrier also implies that the process is independent on the temperature of the system, as will be discussed further below. A detailed study of the time elapsed as the center of the polaron moves from one molecule to the next through the system reveals that after a short acceleration phase, the polaron continues to drift through the system at virtually constant velocity v. The behaviour of v as a function of J for E0 = 5.0 · 104 V/cm is depicted in Fig. 2. It should be stressed that computations become increasingly demanding in the lower region of intermolecular interactions where the time of acceleration grows rapidly with decrements in J. We therefore adopt a slightly different approach in the dynamic simulations for J 6 40 meV and reduce the size of the systems to twenty molecules. The value of J is initialized at 50 meV and the charge carrier is allowed to propagate through the system for 100 fs. Then J is set to the desired value of 20, 30, and 40 meV, respectively. This procedure results in a substantial reduction in CPU time for the simulations. In the limit of weak intermolecular interactions we find charge carriers to be immobile for J 6 20 meV. We also observe an upper boundary of v(J) above which the width of the polaron becomes so large that the concept of a localized charge carrier breaks down. In this regime, i.e. above J 120 meV, charge carrier propagation is best described in terms of band transport. This is in agreement with the discussion above since these values of J exceed the localization energy (El = 97 meV). Finally, in the intermediate regime we find that v(J), although strictly increasing, is a function that is concave down on the interval of 20 < J < 120 meV (see Fig. 2). This regime most probably spans the value of J relevant for pentacene crystals [13]. We therefore conclude that the mobilities between 22 cm2/V s at J = 30 meV and 130 cm2/V s at J = 100 meV (for E0 = 5.0 · 104 V/cm) set an upper limit on the polaron mobility in pentacene single crystals, since the simulations are performed for a defect free system. Since we are dealing with systems at low temperatures it could be argued that a discretization of phonons imposed by a quantum mechanical description of the lattice is crucial for an accurate description of charge transport.
M. Hultell, S. Stafstro¨m / Chemical Physics Letters 428 (2006) 446–450
DR [A]
0.06
449
a
0.04
5
6
7
8
9
10
11
12
13
5
6
7
8
9
10
11
12
13
0.02 0
b
DQ [qe]
0.3 0.2 0.1 0 100
150
200
250
300
Time [fs] Fig. 1. The time evolution of: (a) total molecular bond-length deviation (DR) from the initial state geometry and (b) net charge (DQ) per molecule for J = 50 meV and E0 = 5.0 · 104 V/cm. Only the curves of molecules 5–13 are highlighted. 0.7 0.6
v [Å/fs]
0.5 0.4 0.3 0.2 0.1 0 0
20
40
60
80
100
120
J [meV]
Fig. 2. The dependence of J on carrier velocity v for applied field strength E0 = 5.0 · 104 V/cm.
This is true in the nonadiabatic limit of weak intermolecular interactions and may explain the observed immobility of charge carriers below J 6 20 meV. However, increasing the lattice kinetic energy using nonzero initial atomic velocities ð_rn Þ so as to simulate a temperature T ¼ hM r_ 2n i=3k B ¼ 100 K only produced negligible corrections to the functional dependence of v(J) on the interval 20 < J < 120 meV1. This behaviour is consistent with the drift process with very low activation energy barriers for polaronic motion previously noted and indicate that a classical description of the lattice is indeed sufficient for the previously derived results to be reliable. In this context though it is important to point out that J is kept constant and uniform throughout our simulations. This treatment is used to steady the intrinsic nature of polaron transport in ordered molecular crystals. In reality, 1
Note that in this case the value of k is set to zero in order to preserve the temperature by avoiding heat to dissipate from the system. Comparative simulations show that this will not effect the transport dynamics.
however, molecular motion modulation of the intermolecular transfer integrals results in a power law temperature dependence of the mobility [21,22]. In the future we aim to include also time-dependent intermolecular interactions in our model to study the transport dynamics of heated systems. Despite the described situation of a polaron drift rather than polaron hopping depicted in Fig. 1, transport in the region of low J can be described as nonadiabatic, i.e., the dynamics requires multiple electronic states. We can obtain a very direct view of the nonadiabatic behavior of the polaron dynamics by studying the time evolution of the occupation number of the instantaneous eigenstates (see Eq. 8 above). Fig. 3a shows the occupation numbers for J = 50 meV and E0 = 5.0 · 104 V/cm. The occupancy after the field has reached its maximum value E0 is in the higher region of the unoccupied spectrum of the manifold of pentacene LUMO levels. Due to the weak intermolecular interaction relative to the localization energy, the electronic states are localized. The states that localize to the high (low) potential region of the system are Stark shifted upwards (downwards) in the spectrum. The electron is injected into the high potential energy region and will consequently occupy the higher lying states. Thus, the electron is initially in an excited state in the instantaneous picture (e.g. at t = 100 as in Fig. 3a). The motion towards the low potential region has a driving force in terms of the change in the occupancy of the instantaneous eigenstates towards lower energy as a function of time. This is clearly seen from Fig. 3a. The process is driven by the external electric field and is essentially independent on the lattice vibrations since transport occurs as a result of polaron drift rather than polaron hopping. The lattice kinetic energy is instead generated as a result of the geometrical response of the molecule to the presence of electric charge. Fig. 3b shows the occupancy of the instantaneous eigenstates in the case of J = 100 meV and E0 = 5.0 · 104 V/cm. Clearly, the behavior is different from that shown in Fig. 3a.
450
M. Hultell, S. Stafstro¨m / Chemical Physics Letters 428 (2006) 446–450
strength J and localization energy El. Using the pentacene single crystal as a model system, we show that in the regime of intermediate J (i.e., J El) there is a smooth transition from a nonadiabatic to an adiabatic polaronic drift process on the region of 20 < J < 120 meV. For intermolecular interaction strength above 120 meV the polaron is no longer stable and the transport becomes band like. We thus cover a regime between J El, where text book band theory applies, and J El, where the polaron localize to a single molecule and transport is due to a nonadiabatic hopping process. Acknowledgement Financial support from the Center of Organic Electronics (COE), Swedish Foundation of Strategic Research, is gratefully acknowledged. References
Fig. 3. The time evolution of the occupation number at E = 5.0 · 104 V/ cm, (a) J = 50 meV and (b) J = 100 meV.
The occupancy is now almost constant over time indicating that the transport is adiabatic [18]. Still, however, we are dealing with a localized charge carrier. The difference as compared to the low J case is that all other states are delocalized and therefore do not exhibit such strong Stark shift as in the previous case. The only state which is really sensing the external field is the polaron state which itself drifts along the system towards the low potential energy regime. Thus, the dynamics is in this case best described as a change in the spatial location of one and the same occupied state as opposed to the nonadiabatic case in which the electron performs a sequence of hopping events between states that are localized to different regions in space. 4. Conclusions In summary, we have introduced an approach to study the transport of charge carriers in perfectly ordered molecular crystals as a function of intermolecular interaction
[1] P. Batail (Ed.), Molecular Conductors of Chemical Reviews, vol. 104, ACS Publications, 2004. [2] V. Podzorov, E. Menard, A. Borissov, V. Kiryukhin, J.A. Rogers, M.E. Gershenson, Phys. Rev. Lett. 93 (8) (2004) 086602. [3] W. Warta, L.B. Stehle, N. Karl, Appl. Phys. A 36 (1985) 163. [4] H. Klauk, M. Halik, U. Zschieschang, G. Schmid, W. Radlik, W. Weber, J. Appl. Phys. 92 (2002) 5259. [5] O.D. Jurchescu, T. Baas, T.T.M. Palstra, Appl. Phys. Lett. 84 (2004) 3061. [6] A. Miller, E. Abrahams, Phys. Rev. 120 (1960) 745. [7] T. Holstein, Ann. Phys. 8 (1959) 706. [8] D. Emin, Phys. Rev. B 61 (21) (2000) 14543. [9] V. May, O. Kuhn, Charge and Energy Transfer Dynamics in Molecular Systems, Wiley-VCH, 2004, Chapter 6. [10] R.A. Marcus, Rev. Mod. Phys. 65 (1993) 599. [11] T. Holstein, Ann. Phys. 8 (1959) 725. [12] J.M. Ziman, Electrons and Phonons, Clarendon, Oxford, 1982. [13] Y.C. Cheng, R.J. Silbey, D.A. da Silva Filho, J.P. Calbert, J. Cornil, J.L. Bre´das, J. Chem. Phys. 118 (8) (2003) 3764. [14] R.S. Mulliken, C.A. Reike, D. Orloff, H. Orloff, J. Chem. Phys. 17 (1949) 1248. [15] A. Hansson, S. Stafstro¨m, Phys. Rev. B 67 (2003) 075406. ˚ . Johansson, S. Stafstro¨m, Phys. Rev. Lett. 86 (2001) 3602. [16] A [17] S. Block, H.W. Streitwolf, J. Phys.; Condens. Matter 8 (1996) 889. [18] H. Nakamura, Nonadiabatic Transition: Concepts, Basic Theories and Applications, World Scientific, Singapore, 2002. [19] R.B. Campbell, J. Monteath Robertson, J. Trotter, Acta Cryst. 14 (1961) 705. [20] A.A. Johansson, S. Stafstro¨m, Phys. Rev. B 69 (2004) 235205. [21] K. Hannewald, P.A. Bobbert, Appl. Phys. Lett. 85 (2004) 1535. [22] A. Troisi, G. Orlandi, Phys. Rev. Lett. 96 (2006) 086601.