Supporting information Microscopic simulations of charge transport in disordered organic semiconductors Victor R¨ uhle,1 Alexander Lukyanov,1 Falk May,1 Manuel Schrader,1 Thorsten Vehoff,1 James Kirkpatrick,2 Bj¨orn Baumeier,1, 3 and Denis Andrienko1, 3, 4 1
Max Planck Institute for Polymer Research, Ackermannweg 10, 55128 Mainz, Germany 2 Oxford Centre for Collaborative Applied Mathematics (OCCAM), University of Oxford, St Giles’ 24-29, OX1 3LB Oxford, UK 3 Institute for Pure and Applied Mathematics, University of California Los Angeles, 460 Portola Plaza, Los Angeles, CA 90095, USA 4 Center for Organic Photonics and Electronics and School of Chemistry and Biochemistry, Georgia Institute of Technology, Atlanta, Georgia 30332, USA∗ (Dated: August 24, 2011)
I.
(a)
BIMOLECULAR ELECTRON TRANSFER RATE
(b)
(c)
FIG. 1: (a) Potential energy surfaces of the charge transfer complex in a dimer representation. ET is from molecule i to molecule j. In the initial state, |I00 ⟩, both molecules are in their vibrational ground states. In the final state, |Fl′ m′ ⟩, the neutral molecule i is in vibrational state l′ , while the charged molecule j is in vibrational state m′ . Initial and final states are coupled to a classical harmonic outer-sphere normal mode with mass weighted average coordinate q and reorganization energy λout ij . For small couplings VI00 Fl′ m′ the ET reaction takes place on the diabatic states (solid curves). (b) PES of molecule i as a function of the averaged normal mode qi . l and l′ enumerate vibrational modes of the initial charged and the final neutral states. (c) Same as (b) for initially neutral molecule j. ∆Ui (∆Uj ) is the internal energy difference while λcn (λnc i j ) is the intramolecular reorganization energy for discharging molecule i (charging molecule j).
In the case of a bimolecular electron transfer (ET) reaction the electron moves between two independent molecules. Therefore, one needs separate sets of coordinates for the donor and acceptor. Strictly speaking, the classical Marcus rate (see sec. Introduction of the main text) assumes a common set of vibrational coordinates and, as such, can not be used for bimolecular ET. Yet if the independent vibrational modes are harmonic, are treated classically, and the charging and discharging reorganization energies of the molecule are identical, one still obtains the Marcus-type ET rate with the intramolecular reorganization energy which is the sum of the reorganization energies of the donor and the acceptor [1]. Similarly, the classical treatment of the outer-sphere mode, which is due to rearrangement of the surrounding, allows to add its reorganization energy to the intramolecular one. However, the main issue with the classical Marcus rate is that the high-frequency intramolecular vibrational modes are energetically comparable to the C-C bond stretching mode. At room temperature h ¯ ωCC ∼ 0.2 eV ≫ kB T ∼ 0.025 eV and therefore these modes should be treated quantum mechanically. In fact, for a common set of intramolecular high-frequency (quantum-mechanical) and an outer sphere low-frequency (classical) vibrational coordinates, a
∗ Electronic
address:
[email protected] 2 mixed quantum-classical multi-channel generalization of the Marcus formula is readily available [1]. Such generalization, to the best of our knowledge, has not been made for the bimolecular ET rate, which requires independent sets of coordinates for donor and acceptor. The purpose of this section is to derive a quantum-classical expression for the ET rate with two independent, high-frequency vibrational modes and a common low-frequency outer sphere mode. Following Ref. [2] we assume that all intramolecular modes of a donor i can be averaged into a mode with mass weighted coordinate qi and energy h ¯ ωin (¯hωic ) for the molecule in a neutral (charged) state. Similar assumptions are made for the acceptor j. In addition, we allow for an averaged classical outer-sphere mode with mass weighted out coordinate q and energy h ¯ ωij ≪ kB T . This mode is common to both molecules and plays the role of the ET reaction coordinate [3]. In amorphous organic semiconductors the electronic coupling is usually small compared to both the energy of the classical vibrational mode and intermolecular reorganization energies. In this case the initial, |Ilm ⟩, and final, |Fl′ m′ ⟩, states of the ET reaction are diabatic (non-interacting) dimer states which depend on the vibrational states (with quantum numbers l, m, l,′ m′ ) of both molecules. The potential energy surfaces (PES) corresponding to these states are shown in figure 1a. The PES for intramolecular degrees of freedom for molecules i and j are shown in figure 1b and c, respectively. For the contributions of the outer-sphere mode to initial and final states we introduce Hamiltonian functions HI,F (q) =
]2 1 [ out ω (q − qI,F ) , 2
(1)
where the equilibrium position in the initial (final) state qI (qF ) corresponds to the arrangement of all nuclear coordinates of molecules surrounding the ET complex when molecule i (j) is charged. The outer sphere reorganization 2 1 out energy, defined as λout |qI − qF |] , is shown in figure 1a. It can be computed from the initial and final electric ij = 2 [ω displacement fields of the charge-transfer complex (see sec. Outer-sphere reorganization energy of the main text). The complete Hamiltonian of the ET complex can now be written as
Hij =
∞ ( ∑
∞ ) ( ) ∑ ij ji HI (q) + Elm |Ilm ⟩⟨Ilm | + HF (q) + Em |Fl′ m′ ⟩⟨Fl′ m′ | + ′ l′ l′ ,m′ =0
l,m=0 ij Elm
=
UicC
+
UjnN
ji Em UjcC + UinN ′ l′ =
) ( )] [ ( 1 1 c n +h ¯ ωi l + + + + ωj m + , 2 2 [ ( ) ( )] 1 1 + Ejel + Ejext + h ¯ ωjc m′ + + ωin l′ + . 2 2 Eiel
∑
VIlm Fl′ m′ |Ilm ⟩⟨Fl′ m′ | + h.c ,
l,m,l′ ,m′
Eiext
(2) Here, a manifold of initial states, |Ilm ⟩, with quantum numbers l (m) for intramolecular vibrations in molecule i (j) ij and energy Elm , is coupled to a classical phonon bath HI (q). Transitions to the manifold of final states |Fl′ m′ ⟩ where ij ji , and final, Em the charge has hopped from i to j are possible due to a coupling VIlm Fl′ m′ . The initial, Elm ′ l′ , energies contain internal energies UinN and UicC (UjnN and UjcC ) of molecule i (molecule j) in the neutral and charged ground states as introduced in sec. Intramolecular reorganization energy of the main text, the contributions of the external electric field, Eiext and Ejext , and electrostatic interactions, Eiel and Ejel , as introduced in sec. Site energy difference of the main text, and respective oscillator energies. Within the Born-Oppenheimer approximation, a separation in terms of electronic and nuclear degrees of freedom gives |Ilm ⟩ = |ϕci ⟩|χcil ⟩|ϕnj ⟩|χnjm ⟩ , |Fl′ m′ ⟩ = |ϕni ⟩|χnil′ ⟩|ϕcj ⟩|χcjm′ ⟩ ,
(3)
where ϕni (ϕci ) corresponds to the electronic part of the wave function, while χnil (χcil ) represents an l-th phonon mode of the neutral (charged) molecule i. The coupling element VIlm Fl′ m′ in eq. (2) can then be factorized in an electronic and nuclear parts VIlm Fl′ m′ = Jij ⟨χcil |χnil′ ⟩⟨χnjm |χcjm′ ⟩ .
(4)
The calculation of the electronic part Jij is explained in more detail in sec. Transfer integrals of the main text. Franck-Condon overlap integrals ⟨χcil |χnil′ ⟩ (⟨χnjm |χcjm′ ⟩) describe couplings of vibrational modes l, l′ (m, m′ ) of the charged and neutral configurations of molecule i (j). Exemplary modes are shown in figure 1b,c.
3 Since kB T ≪ ¯ hωic , ¯ hωjn one can restrict the initial state to the vibrational ground-states l = m = 0 while allowing tunneling to all vibrationally excited states l′ for molecule i and m′ for molecule j. In other words, a single initial state |I00 ⟩ couples to a manifold of final states |Fl′ ,m′ ⟩. This assumes that ET is sufficiently slow compared to the relaxation of the intramolecular degrees of freedom, so that there is enough time for a complex to relax to its vibrational ground state between two consecutive ETs. The energy difference driving the reaction to channel l′ m′ therefore is ij ji ∆Elij′ m′ = E00 − Em h(ωin l′ + ωjc m′ ) , ′ l′ = ∆Eij − ¯ ext el int is explained in sec. Site energy difference of the main text. where ∆Eij = ∆Eij + ∆Eij + ∆Eij out out Assuming that |VI00 Fl′ m′ | ≪ λij , ¯ hω and using Fermi’s golden rule with VI00 Fl′ m′ as a perturbation to the initial diabatic state, we obtain a multi-channel rate equation
ωij =
∫ ∞ ∑ 2π |VI00 Fl′ m′ |2 dqfI (q)δ(∆Elij′ m′ + HI (q) − HF (q)) . h ¯ ′ ′
(5)
l ,m =0
where the thermal averaging over the classical outer-sphere mode is performed by introducing a canonical distribution ∫ function fI (q) = Z −1 exp(−HI (q)/kB T ), with Z = dq exp(−HI (q)/kB T ). Energy conservation pins the transition to the crossing point of the diabatic PES (see figure 1a) resulting in 2π |Jij |2 √ ωij = ¯h 4πλout ij kB T
∞ ∑
|⟨χci0 |χnil′ ⟩|2 |⟨χnj0 |χcjm′ ⟩|2
l′ ,m′ =0
{ [ ]2 } ∆Eij − ¯h(l′ ωin + m′ ωjc ) − λout ij . exp − 4λout ij kB T
(6)
Eq. (6) is the quantum-classical expression for the bimolecular ET rate with two independent, high-frequency vibrational modes and one classical common outer-sphere mode. It is the main result of this section. If the curvatures of intramolecular PES of charged and neutral states of a molecule are different, that is ωic ̸= ωin , 1 1 n n c 2 nc c n c 2 the corresponding reorganization energies, λcn i = 2 [ωi (qi − qi )] and λi = 2 [ωi (qi − qi )] , will also differ. In this case the Franck-Condon (FC) factors for discharging of molecule i read [4] √
|⟨χci0 |χnil′ ⟩|2 =
) 2 ( )k/2 l′ ( ′ ) ( c ∑ k! l 2ωi si exp (−|si |) Hl′ −k √ cn , ωic + ωin (k/2)! k 2Si
ωic ωin 2 ′ c 2l l′ ! (ωi + ωin )
(7)
k=0 k even
√ 2 λnc λcn where Hn (x) is a Hermite polynomial, si = h¯ (ωci+ωni ) , and Sicn = λcn hωic . The FC factors for charging of molecule i /¯ i i cn c nc n j can be obtained by substituting (si , Si , ωi ) with (−sj , Sj , ωj ). In order to evaluate the FC factors, the internal reorganization energy λcn i can be computed from the intramolecular PES, as shown in figure 1b,c and explained in sec. Intramolecular reorganization energy of the main text. To conclude the section, we compare the bimolecular quantum-classical rate, eq. (6), the classical bimolecular Marcus rate, eq. (1) of the main text, and the quantum-classical Jortner rate with a common set of vibrational coordinates [1] ∞ ∑ 2π |Jij |2 1 √ ωij = h ¯ N! 4πλout ij kB T N =0
(
λint ij hω int ¯
)N
(
λint ij exp − int ¯hω
)
{ [ ]2 } ∆Eij − ¯hN ω int − λout ij exp − . 4λout ij kB T
(8)
cn If ωic = ωin = ωi (λnc i = λi = λi ), the Franck-Condon factor simplifies to
1 |⟨χi0 |χ ⟩| = ′ l! il′
2
(
λi ¯hωi
)l ′
(
λi exp − ¯hωi
) .
(9)
If this simplification is applicable for both donor and acceptor molecules, eq. (6) becomes identical to the quantumclassical rate eq. (8) with λint ij = λi + λj .
4
√ −2 FIG. 2: (a) Scaled hopping rates, ω ¯ ij = ωij Jij (2π)−1 ¯ h 4πkB T , calculated using the classical Marcus, Jortner quantum mechanical eq. (8) and bimolecular multichannel eq. (6) rate expressions. Outer sphere reorganization energy λout ij = 0.05 eV, cn nc nc λcn i = λj = 0.14 eV, and λi = λj = 0.09 eV (all added in the classical Marcus rate while the latter two are added for the Jortner rate). Intramolecular hωiint = ¯ hωjint = 0.2 eV for the Jortner rate while ¯ hωin = 0.2 eV p p cn vibrations have averaged frequency ¯ out / λ for the bimolecular rate. (b) The same but for λ = 0.1 eV. (c) Histogram of rates at a field of and h ¯ ωjc = ¯ hωin λnc ij j j 108 Vm−1 for the Marcus and Jortner rates with distance dependent λout < 0.08 eV for the neighborlist pairs (see sec. Outerij sphere reorganization energy of the main text) and constant λint = 0.23 eV. A small difference can be seen in the tail of small rates.
To compare the quantum-mechanical and classical rates, intramolecular hole reorganization√energies of Alq3 , λcn i = n c n nc /λcn . Due to the 0.14 eV and λnc = 0.09 eV were used. We also assumed that ω = 0.2 eV and ω = ω λ i i i i i i out out uncertainty in determining λout ij , two cases are considered, λij = 0.05 eV (figure 2a) and λij = 0.10 eV (figure 2b). Note that the estimate made in the main text predicts λout < 0.08 eV. In both cases we used fixed, molecularij separation independent, λout . ij Figure 2 shows that the main difference between the quantum-classical and classical rates is the tail of smaller rates for large negative ∆E (endothermal hopping) and higher rates for large positive ∆E (exothermal hopping). Figure 2c also shows the corresponding distributions of rates for all pairs from the neighbor list for 512 molecules of amorphous Alq3 . Here we used the distance-dependent λout ij from the neighbor list as computed from dielectric displacement fields (see sec. Outer-sphere reorganization energy of the main text) with the Pekar factor of cp = 0.01. One can see that the distributions are practically on top of each other (except for very small rates) and hence will lead to similar charge dynamics. In general, our observation is that for a situation with (i) intramolecular reorganization energy similar to the outer out sphere one (λint ij ∼ λij ), (ii) driving force ∆Eij small compared to the intramolecular reorganization energy, and (iii) λ ∼ ¯hω, the classical (eq. (1) of the main text) and semi-classical (eq. (8) and eq. (6)) expressions lead to quantitatively similar rates. However, for systems with large ∆Eij , such as donor-acceptor mixtures, eq. (6) or eq. (8) should be used. In this case a rather accurate estimate of the outer sphere reorganization energy is required [5]. II.
FORCE FIELD
B3LYP functional and the GAUSSIAN package [6] were used in all density functional calculations. The geometry optimization of the meridional [7] isomer, whose structure is shown in figure 3, was performed using the 6-311g(d) basis set and partial charges were obtained using the CHELPG method [8]. The Lennard-Jones parameters were taken from the OPLS [9] force field. Bonded interactions were parametrized by matching the first-principle and the force field based potential energy surface (PES) scans [10, 11], as described in sec. Material morphology of the main text. Since the meridional isomer of Alq3 is asymmetric, e.g. Oa-Al-Ob and Ob-Al-Oc angles as well as Oa-Al-Ob-Cb1 and Ob-Al-Oc-Cc1 dihedrals are different (atom names are explained in figure 3), in total 16 scans (6-311g+(d,p) basis set) were performed. The PES are shown in figure 4 for several representative degrees of freedom. For angle and dihedral potentials, a quadratic dependence, V (θ) = 21 kθ (θ − θ0 )2 , was used for fitting. To model the ionic Al-N bond, a harmonic spring potential was used. All three bonded potentials (Al-Na, Al-Nb, Al-Nc) were parameterized based on the Al-Na scan (see figure 4). The resulting constants are summarized in table I.
5
Na
Oc N
Al
Al
Oa Ob
O
Ca2 Ca1
Ca3
3
FIG. 3: Chemical structure and meridional isomer of tris-8(hydroxyquinoline aluminium) (Alq3 ). Atom labels correspond to the atom names used in the force-field.
Ca1-Oa-Al
Oa-Al-Ob
Al-Oa-Ca1-Ca2
6
energy, kJ/mol
4
B3LYP FF
10
3
4
2 5
2 1
0 110
115
120
0 90
125
95
100
angle, deg
angle, deg
Oa-Ca1-Ca2-Ca3
Oa-Al-Ob-Cb1
105
0
170
180
190
angle, deg Al-Na 100
4
energy, kJ/mol
4 75
3 3 2 1 0 170
50
2
25
1 180
190
angle, deg
0 -100
-90
-80
-70
0 0.15
angle, deg
0.2
0.25
bond, nm
FIG. 4: Potential energy scans (full lines, B3LYP functional, 6-311+g(d,p) basis set) and the corresponding fits based on the OPLS force field (dashed lines) for several representative degrees of freedom. Atom labeling is shown in figure 3.
The force-field files (in GROMACS format) are part of the supporting information.
III.
CONFORMATIONAL DISORDER
As already mentioned in sec. Intramolecular reorganization energy of the main text, the three ligands (denoted by letters a, b, and c in figure 3) of Alq3 can easily change their mutual orientations. Molecular conformations are then ‘frozen’ due to non-bonded interactions in an amorphous glass. Particularly ‘soft’ are six dihedral angles: Al-Oa-Ca1-Ca2, Al-Ob-Cb1-Cb2, Al-Oc-Cc1-Cc2, Oa-Al-Ob-Cb1, OaAl-Oc-Cc1, Oc-Al-Oa-Ca1 which determine mutual orientations of ligands (atom names are shown in figure 3, while
6 Degree of freedom Ca1-Oa-Al Cb1-Ob-Al Cc1-Oc-Al Oa-Al-Ob Ob-Al-Oc Oc-Al-Oa Al-Oa-Ca1-Ca2 Al-Ob-Cb1-Cb2 Al-Oc-Cc1-Cc2 Oa-Ca1-Ca2-Ca3 Ob-Cb1-Cb2-Cb3 Oc-Cc1-Cc2-Cc3 Oa-Al-Ob-Cb1 Oa-Al-Oc-Cc1 Oc-Al-Oa-Ca1 bond Al-Na
θ0 , deg 120 120 122 95.5 170 94 180 172 178 180 180 180 -80.42 -174 89.5
kθ , kJmol−1 rad−2 600 500 400 320 200 420 55 50 50 200 200 200 50 30 20
0.208 nm 62459, kJmol−1 nm−2
TABLE I: Parameters for bonded interactions obtained from PES scans. Atom labeling is shown in figure 3.
int FIG. 5: Distributions of (a) internal reorganization energy λint ij and (b) internal energetic disorder ∆Eij for a system of 512 molecules. Geometry optimizations (B3LYP functional, 6-311G(d,p) basis set) in charged and neutral states are performed for every molecule with six ‘soft’ dihedral angles constrained to their average values.
table I summarizes the corresponding energies). int In order to evaluate the variation of λint ij and ∆Eij due to conformational changes, the angles of these six dihedrals were averaged, for every molecule, over 100 snapshots of a 2 ns trajectory. A simulation box of total 512 molecules was used. The internal energies UinC , UinN , UicN , and UicC were calculated for every molecule i after optimizing its molecular geometry (B3LYP functional, 6-311G(d,p) basis set) in charged and neutral states with the six dihedral angles constrained to their average values.
7 A.
Internal reorganization energy
Distributions of λint ij for electrons and holes are shown in figure 5a. For holes (electrons), the maximum is at 0.21 eV (0.27 eV). Gaussian fit provides variance of 0.03 eV (0.02 eV). Computing λint ij from the PES of two unconstrained molecules leads to a similar value of 0.23 eV (0.28 eV). Since Alq3 has high energetic disorder arising from its large dipole moment, this small variance in reorganization energy does not affect charge carrier mobility or Poole-Frenkel behavior.
B.
Internal energetic disorder
int Distributions of ∆Eij for electrons and holes are shown in figure 5b. One can see that significant conformational int changes lead to an approximately Gaussian distribution of ∆Eij with a small variance of 0.011 eV (0.018 eV) for holes (electrons). For Alq3 , the internal energy disorder is small compared to the electrostatic (including polarization) energetic disorder and hence does not affect the charge carrier mobility.
IV.
MASTER EQUATION FOR MULTIPLE CHARGE CARRIERS
For multiple charge carriers, a rigorous derivation of master equation in terms of occupation probabilities pi is given in Ref. [12]. Here we recapitulate the derivation where Coulomb interactions are neglected and charge-charge interactions are taken into account by imposing the single site occupancy. We also use a mean-field approximation, i. e. ignore correlations between site occupation probabilities. Rewriting the sum over states (α, β) in terms of the sum over sites (i, j) we obtain ∑ ∂Pα ωji αj (1 − αi )Pα(i↔j) − ωij αi (1 − αj )Pα , = ∂t i,j
(10)
i̸=j
where α(i ↔ j) is a state obtained from a state α by interchanging the occupations of the sites i and j. The first term on the right hand side describes the transitions from all other sites to the state α (gain) while the second term corresponds to the transitions from the state α to all other states (loss). Within a mean field approximation, i.e. neglecting all correlations between occupation probabilities of the sites, we can write Pα = P (α1 )P (α2 ) . . . P (αN ) ,
(11)
where, by definition, P (αi = 1) = pi is a site occupation probability and P (αi = 0) = 1 − pi . By substituting eq. (11) into eq. (10) and summing over all αk except k = i we obtain ∑ ∂pi = pj (1 − pi )ωji − pi (1 − pj )ωij , ∂t j
(12)
which is a mean-field version of the master equation written ∑ in terms of site occupation probabilities. The solution of this equation, together with the normalization condition, i pi = m, where m is the number of charges in the system, can be used to obtain occupation probabilities, currents, and charge carrier mobilities in systems with more than one charge carrier per simulation box.
V.
KINETIC MONTE CARLO
As already mentioned in the main text, a combination of the variable step size method (VSSM) [13] and the first reaction method (FRM) is implemented in the code. The main benefit of using FRM is that the event groups associated with every charge are independent if rates do not change. Hence, only the group of events of a selected charge must be updated after this charge moves. A flowchart of the implemented algorithm is shown in figure 6. Its right side shows the the complete workflow, while the left side explains how VSSM and FRM are combined. The processes are first grouped according to the charges,
8 charges
compute total rates and waiting times for all charges
ωi =
P
j
ωi→j ,
τi = t − ωi−1 ln r1
total rates
ω1 , ω2 , . . . , ωm choose shortest waiting time propagate time
τk = min{τi } t = τk
first reaction method − selection of a carrier
τi = t − ωi−1 ln r1 k = mini τi t = τk
choose hopping destination
P max{l | l wkl < wk r2 }
enabled event
disabled event
can hop? yes
charge 1
charge 2
rates
rates
ω12 , ω13 , . . . , ω1n
ω21 , ω23 , . . . , ω2n
no
execute hop update list of hops for carrier k
remove event update waiting time
···
charge m
···
rates
ωm1 , ωm2 , . . . , ωmn
···
···
τk = t − ωk−1 ln r1 variable step size method − selection of a hopping site removes disabled events and updates total rates if needed
···
···
site n
···
site 2
τi = t − ωi−1 ln r1
site 1
ωi→j ,
site n
j
site 2
P
site n
ωi =
site 3
update total rates and waiting times for required charges
site 2
(e.g. due to Coulomb interactions)
site 1
P max{l | l wkl < wk r2 }
recalculate rates if required
FIG. 6: Workflow of the implemented kinetic Monte Carlo method. The right scheme illustrates how the events are grouped and how the variable step size method (VSSM) and the first reaction method (FRM) are combined hierarchically. The main benefit of FRM is that the event groups are independent (if rates do not change). Hence, only the group of events associated with a selected charge must be updated after this charges moves.
i.e. only processes which involve an occupied site are considered. For each charge i, a total escape rate is calculated as ∑ ωi = ωij , (13) j
where the sum is over all neighbors of site i in the neighbor list. At this stage the disabled processes (for example, if the neighboring site is occupied) do not have to be excluded. After this, a particular charge carrier is selected using the first reaction method. To do this, waiting times for each charge carrier are calculated using a random number r1 ∈ (0, 1] τi = t − ωi−1 ln r1 .
(14)
After this, the charge k with the smallest waiting time τk is selected and the system time is advanced to tnew = τk . Then, a particular process l is picked according to VSSM with a probability ωkl /ωk , which is equivalent to choosing the biggest l for which ∑ ωk−1 ωkl ≤ r2 . (15) l
where r2 ∈ (0, 1] is the second random number. If the process is enabled, the charge is moved to the site l. If the process is disabled, the charge remains on the current site and the disabled process is removed from the list of processes. After this, a new waiting time for the charge on the site l is calculated according to τl = tnew − ωl−1 ln r1
(16)
If the rates of the other charges did not change (i.e. no charge-charge interactions), waiting times of the remaining charges are not updated. Finally, the charge k with the smallest waiting time τk is selected again and the cycle repeats.
9 VI.
DIFFUSION IN PERIODIC BOUNDARY CONDITIONS
In a macroscopic limit the probability p(r, t) to find a charge at a position r at time t is a solution of the Smoluchowski equation ∂p(r, t) ˆ −βU (r) ∇eβU (r) p(r, t) , = ∇ · De ∂t
(17)
ˆ is the diffusion tensor, and U (r) is an external potential. where β = 1/kB T , D In our case eq. (17) is supplemented with periodic boundary conditions p(r, t) = p(r + m1 a1 + m2 a2 + m3 a3 , t), where a1 , a2 , a3 are primitive lattice vectors, and m1 , m2 , and m3 are integer numbers. In the absence of external fields eq. (17) reduces to a diffusion equation ∂p(r, t) ∑ ∂ 2 p(r, t) = Dαβ ∂t ∂rα ∂rβ
(18)
(19)
α,β
where α, β = x, y, z. Separating the variables and expanding p(r, t) in Fourier series we obtain ∑ p(r, t) = exp (−λk t) ck exp (ik · r)
(20)
n1 ,n2 ,n3
where k are vectors of the reciprocal lattice k = n1 g1 + n2 g2 + n3 g3 , a3 × a1 a1 × a2 a2 × a3 , g2 = 2π , g3 = 2π , g1 = 2π a1 · (a2 × a3 ) a2 · (a3 × a1 ) a3 · (a1 × a2 ) n1 , n2 , n3 are integer numbers, and λk =
∑
Dαβ kα kβ .
(21)
(22)
α,β
If at t = 0 the charge is injected at a site with the coordinate rm then 1 ∑ p(r, t = 0) = δ(r − rm ) = exp [ik · (r − rm )] . V n ,n ,n 1
2
(23)
3
Comparison to eq. (20) yields ck =
1 exp (−ik · rm ) V
(24)
and p(m) (r, t) =
∑ n1 ,n2 ,n3
exp (−λk t)
exp [ik · (r − rm )] . V
The Fourier coefficients of the occupation probability then read ∫ M 1 ∑ exp (ik · rm ) p(m) (r, t) exp (−ik · r) dV ⟨pk ⟩ = M m=1
(25)
(26)
where we additionally averaged over M charge injections. An estimate of this integral for an irregular lattice can be obtained by using the effective volume per site j, Vj , obtained from a Voronoi tessellation of space. For reasonably uniform lattices (uniform site densities) this volume is almost independent of the site and a constant Vj = V /N can be assumed ⟨pk ⟩ =
M 1 ∑ ik·rm V ∑ (m) e p (rj , t)e−ik·rj M m=1 N j
(27)
∑ To determine the diffusion tensor, the time dependence of pk can be fitted to exp(− α,β Dαβ kα kβ t). The modes with the smallest wavevectors k, such as g1 , g2 , and g3 will provide the best estimate of Dαβ .
10 References [1] May, V.; K¨ uhn, O. Charge and Energy Transfer Dynamics in Molecular Systems, 2nd ed.; Wiley-VCH Verlag GmbH & Co. KGaA, 2003. [2] Br´edas, J.; Beljonne, D.; Coropceanu, V.; Cornil, J. Chem. Rev. 2004, 104, 4971–5004. [3] In a more general case the curvature of PES of this mode might change depending on whether molecule i or j is charged, out out out i.e. ωij ̸= ωji and thus λout ij ̸= λji . If the difference is large one has to combine our derivation with the argument of Ref. [14]. [4] Chang, J. J. Mol. Spectrosc. 2005, 232, 102–104. [5] McMahon, D. P.; Troisi, A. J. Phys. Chem. Lett. 2010, 1, 941–946. [6] Frisch, M. J. et al. Gaussian 03, Revision C.02 ; 2004; Gaussian, Inc., Wallingford, CT, 2004. [7] Another isomer which is commonly observed in the octahedral complexes of the type MN3 O3 , where M is a trivalent metal, namely the facial isomer, has not been experimentally identified for amorphous Alq3 [15]. [8] Breneman, C. B.; Wiberg, K. E. J. Comp. Chem. 1990, 11, 361–373. [9] Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657–1666. [10] R¨ uhle, V.; Kirkpatrick, J.; Kremer, K.; Andrienko, D. Phys. Stat. Solidi B 2008, 245, 844. [11] Marcon, V.; Vehoff, T.; Kirkpatrick, J.; Jeong, C.; Yoon, D. Y.; Kremer, K.; Andrienko, D. J. Chem. Phys. 2008, 129, 094505. [12] Cottaar, J.; Bobbert, P. A. Phys. Rev. B 2006, 74, 115204. [13] Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. J. Comput. Phys. 1975, 17, 10–18. [14] Casado-Pascual, J.; Goychuk, I.; Morillo, M.; H¨ anggi, P. Chem. Phys. Lett. 2002, 360, 333–339. [15] C¨ olle, M.; Br¨ utting, W. Phys. Status Solidi A 2004, 201, 1095–1115.
