저작자표시 - 비영리 - 변경금지 2.0 대한민국 이용자는아래의조건을따르는경우에한하여자유롭게 이저작물을복제, 배포, 전송, 전시, 공연및방송할수있습니다. 다음과같은조건을따라야합니다 : 저작자표시. 귀하는원저작자를표시하여야합니다. 비영리. 귀하는이저작물을영리목적으로이용할

Transcription

1 저작자표시 - 비영리 - 변경금지 2.0 대한민국 이용자는아래의조건을따르는경우에한하여자유롭게 이저작물을복제, 배포, 전송, 전시, 공연및방송할수있습니다. 다음과같은조건을따라야합니다 : 저작자표시. 귀하는원저작자를표시하여야합니다. 비영리. 귀하는이저작물을영리목적으로이용할수없습니다. 변경금지. 귀하는이저작물을개작, 변형또는가공할수없습니다. 귀하는, 이저작물의재이용이나배포의경우, 이저작물에적용된이용허락조건을명확하게나타내어야합니다. 저작권자로부터별도의허가를받으면이러한조건들은적용되지않습니다. 저작권법에따른이용자의권리는위의내용에의하여영향을받지않습니다. 이것은이용허락규약 (Legal Code) 을이해하기쉽게요약한것입니다. Disclaimer

2 Charge Transport in Organic Photovoltaic Cells

3

4 이학박사 학위논문 Charge Transport in Organic Photovoltaic Cells 유기태양전지에서 전하 전송 연구 2014년 2월 서울대학교 대학원 화학부 물리화학전공 마누엘 슈라더

5

6 Johannes-Gutenberg-Universität Mainz Fachbereich Physik Charge Transport in Organic Photovoltaic Cells Dissertation zur Erlangung des akademischen Grades Doktor der Naturwissenschaften (Doctor rerum naturalium) verfasst und vorgelegt von Manuel chrader geb. in Wiesbaden Max-Planck-Institut für Polymerforschung Mainz, Juni 2013

7

8 본 논문은 유기 태양 전지의 전하 전송에 대한 심도 있는 이해를 추구하였다. 컴퓨터시뮬레이션을통하여유기 물질에서 전하 전달체의 역학을, 미시적 수 준의양자역학적 과정에서부터 전하 전달체의 이동도를 정량화할 수 있는 거 시적 수준까지재구성하였다. 이다중 스케일적접근 방법에 의하여, 유기 물질 의 화학 구조와 거시적 이동도의 관계 (구조-물성 관계)를 확립하였는데, 이 관 계는 태양광 효율의 개선을 지원하게 된다. 시뮬레이션 모델에는 다음 세 가지 주요구성 요소가포함되어 있다. 첫째는 형태로서, 해당 물질 내에서 분자배 열 모형을 원자 단위로 쪼개어 구성하였다. 둘째는 전하 전송의 호핑 모델로서, 전하이동을 개별 분자사이에서의 연속적인전하 전달반응으로 설명하였다. 마지막은 전하 전달의 비단열 모델로서, 전이율을 다음 세가지파라미터로 설 명하였다: 재구성 에너지, 사이트에너지, 전달 인테그럴. 전하 전송 시뮬레이션은 다이시아노비닐 치환 올리고싸이오펜의 물질적 인 부분과단결정 및 박막과 비정질/스멕틱 메조상의 형태에초점을 맞추었다. 이에따른 일반적결과는, 어셉터-도너-어셉터 순서와 유연한 올리고머주사슬 로구성된 분자구조가 분자의쌍극자모멘텀에 변화를줌으로써사이트에너 지를 변화시킨다는 것이다. 이에너지 측면에서의 무질서는 보통 결정에서 높 으며메조상에서는 더욱 높은 것으로드러났다. 단결정의 경우, 파이스태킹을 갖춘 결정 구조와 그에따른 대규모 전자 전달 인테그럴이 상대적으로 낮은 이 동도로 이어졌다. 이 반직관적인 행동은 에너지 결함이 발생하기 쉬운 전송 경 로의 형성에 기인한것이다. 박막의 경우, 위 추론이 다시 확인됨으로써 실험적 이동도에대한 미시적 이해로이어질 수 있다. 사실, 시뮬레이션결과는측정된 이동도와 태양광 효율 모두와 연관된다. 비정질/스멕틱 계의 경우, 에너지무질 서는올리고머의 길이에따라 증가하는데, 보다스멕틱 질서가커질수록 이동 도가 감소하는 이상 현상을 보인다. 그이유는 스멕틱 층이에너지무질서의 공 간적 상관관계와 충돌하기 때문으로 설명된다. 주요어: 태양전지, 유기, 전하 전송, 이동도, 시뮬레이션, 올리고싸이오펜 학번:

9

10 Zusammenfassung Die vorliegende Dissertation dient dazu, das Verständnis des Ladungstransportes in organischen olarzellen zu vertiefen. Mit Hilfe von Computersimulationen wird die Bewegung von Ladungsträgern in organischen Materialien rekonstruiert, und zwar ausgehend von den quantenmechanischen Prozessen auf mikroskopischer Ebene bis hin zur makroskopischen kala, wo Ladungsträgermobilitäten quantifizierbar werden. Auf Grundlage dieses skalenübergreifenden Ansatzes werden Beziehungen zwischen der chemischen truktur organischer Moleküle und der makroskopischen Mobilität hergestellt (truktur-eigenschafts-beziehungen), die zu der Optimierung photovoltaischer Wirkungsgrade beitragen. Das imulationsmodell beinhaltet folgende drei chlüsselkomponenten. Erstens eine Morphologie, d. h. ein atomistisch aufgelöstes Modell der molekularen Anordnung in dem untersuchten Material. Zweitens ein Hüpfmodell des Ladungstransportes, das Ladungswanderung als eine Abfolge von Ladungstransferreaktionen zwischen einzelnen Molekülen beschreibt. Drittens ein nichtadiabatisches Modell des Ladungstransfers, das Übergangsraten durch drei Parameter ausdrückt: Reorganisationsenergien, Lageenergien und Transferintegrale. Die Ladungstransport-imulationen richten sich auf die Materialklasse der dicyanovinyl-substituierten Oligothiophene und umfassen Morphologien von Einkristallen, Dünnschichten sowie amorphen/smektischen Mesophasen. Ein allgemeiner Befund ist, dass die molekulare Architektur, bestehend aus einer Akzeptor-Donor-Akzeptor- equenz und einem flexiblen Oligomergerüst, eine erhebliche Variation molekularer Dipolmomente und damit der Lageenergien bewirkt. Diese energetische Unordnung ist ungewöhnlich hoch in den Kristallen und umso höher in den Mesophasen. Für die Einkristalle wird beobachtet, dass Kristallstrukturen mit ausgeprägter π-tapelung und entsprechend großer Transferintegrale zu verhältnismäßig niedrigen Mobilitäten führen. Dieses Verhalten wird zurückgeführt auf die Ausbildung bevorzugter Transportrichtungen, die anfällig für energetische törungen sind. Für die Dünnschichten bestätigt sich diese Argumentation und liefert ein mikroskopisches Verständnis für experimentelle Mobilitäten. In der Tat korrelieren die imulationsergebnisse sowohl mit gemessenen Mobilitäten als auch mit photovoltaischen Wirkungsgraden. Für die amorphen/smektischen ysteme steigt die energetische Unordnung mit der Oligomerlänge, sie führt aber auch zu einer unerwarteten Mobilitätsabnahme in dem stärker geordneten smektischen Zustand. Als Ursache dafür erweist sich, dass die smektische chichtung der räumlichen Korrelation der energetischen Unordnung entgegensteht.

11

12 Abstract This thesis serves to deepen the understanding of charge transport in organic photovoltaic cells. Using computer simulations, the dynamics of charge carriers in organic materials is reconstructed, starting from the quantum mechanical processes on the microscopic level up to the macroscopic scale, where charge carrier mobilities can be quantified. Based on this multiscale approach, relationships between the chemical structure of organic molecules and the macroscopic mobility are established (structure-property relationships), which assist the improvement of photovoltaic efficiencies. The simulation model includes the following three key components. First, a morphology, i.e., an atomistically resolved model of the molecular arrangement within the material of interest. econd, a hopping model of charge transport, describing charge migration as a succession of charge transfer reactions between individual molecules. Third, a nonadiabatic model of charge transfer, expressing transition rates by three parameters: reorganization energies, site energies, and transfer integrals. The charge transport simulations focus on the material class of dicyanovinyl-substituted oligothiophenes and cover morphologies of single crystals, thin films, and amorphous/smectic mesophases. A general result is that the molecular architecture, consisting of an acceptor-donor-acceptor sequence and a flexible oligomer backbone, gives rise to substantial variations of molecular dipole moments and hence of the site energies. This energetic disorder is unusually high in the crystals and even higher in the mesophases. For the single crystals, it is observed that crystal structures with a pronounced π-stacking and correspondingly large transfer integrals lead to relatively low mobilities. This counterintuitive behavior is traced back to the formation of preferred transport directions which are prone to energetic defects. For the thin films, this reasoning can be confirmed and provides a microscopic understanding for experimental mobilities. In fact, the simulation results correlate with both measured mobilities and photovoltaic efficiencies. For the amorphous/smectic systems, the energetic disorder increases with the oligomer length, but also leads to an unexpected mobility reduction in the more ordered smectic state. The reason for this is elucidated by showing that the smectic layering conflicts with the spatial correlations of the energetic disorder. Keywords: solar cell, organic, charge transport, mobility, simulation, oligothiophene tudent number:

13

14 Contents Introduction 15 Related Publications 19 Notation 21 I. Background Organic Photovoltaic Cells Electronic tructure of Organic olids Conversion of olar Radiation into Electric Power Power Conversion Efficiency Foundations of Computational Chemistry eparation of Nuclear and Electronic Motions Electronic tructure Theory Hartree-Fock Theory Kohn-ham Density Functional Theory Basis et Approximation II. Methodology Charge Transfer in Molecular ystems Regimes of Charge Transfer Adiabatic Charge Transfer Nonadiabatic Charge Transfer Charge Transfer Rates Bimolecular High-Temperature Nonadiabatic Charge Transfer Further Limits of Nonadiabatic Charge Transfer

15 14 Contents 4. Charge Transport in Organic olids Regimes of Charge Transport Ordered Organic olids Disordered Organic olids Charge Transport imulations in Disordered Organic olids Force Field Morphology Hopping ites Transfer Integrals ite Energies Reorganization Energies Charge Dynamics Macroscopic Observables III. Results Charge Transport imulations in Organic Crystals Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals Morphological Disorder Charge Transfer Parameters Charge Carrier Mobility Dicyanovinyl-ubstituted Quaterthiophenes: Thin Films Crystal tructure Analysis Charge Carrier Mobility Charge Transport imulations in Organic Mesophases Dicyanovinyl-ubstituted Oligothiophenes: Amorphous/mectic Morphological Disorder Charge Transfer Parameters Charge Carrier Mobility Electric Current Pathways Conclusion and Outlook 151 Bibliography 157

16 Introduction The sun provides more energy to the earth everyhourthanmankindconsumesin an entire year. In fact, the energy resource of terrestrial solar radiation far exceeds that of all other renewable and fossil energy sources combined. 1 Harnessing the immense solar energy resource not only has the potential to accommodate the increasing global energy demand, but also holds promise to reshape the energy sector for environmental sustainability. However, a widespread adoption of photovoltaic electricity generation is only achievable through competitive pricing on the energy market. In fact, conventional inorganic photovoltaic cells, although technologically advanced, are still limited to niche applications due to high costs. The emerging technology of organic photovoltaic cells,incontrast,couldquicklyfindaubiquitousdeploymentsince organic materials offer strong potential for cost reduction. In addition to an inexpensive production, organic solar cells can inherit the advantageous physical properties of organic materials, such as light weight and mechanical flexibility. Although organic photovoltaic technology is still far from the level of maturity required to deliver these promises, the field has recently experienced such a rapid progress that it is currently transitioning from a phase of technology development to industrial production. This dynamic development is the fruit of concerted efforts in several areas, such as synthetic chemistry, producing increasingly fine-tuned organic compounds, and material processing, constantly adapting to the demands of the field. Now that organic photovoltaics is close to first commercialization, the scientific community is more than ever demanded to address the still major challenges ahead. The most critical issues of organic solar cells, as compared to their inorganic counterparts, are their shorter life spans and lower power conversion efficiencies. One of the greatest difficulties in improving such device properties is the widely lacking comprehension of how these properties are linked to the constituent organic compounds. As a result, the chemical synthesis of new or modified compounds is mostly guided by intuitive rather than rational design rules. With the aim of a rational compound design, models relating the chemical structures to macroscopic properties, so-called structure-property

17 16 Introduction relationships, become highly desired. Establishing such relationships is a central concern of this work. This necessitates both a microscopic description of organic photovoltaic cells as well as methods for linking the macroscopic properties to this description. A more detailed discussion on the challenges for improving photovoltaic device properties is provided after introducing the required background on organic photovoltaic cells in Chapter 1. Awindowintothemicroscopicworldofchemicalmatterisopenedbythefieldof computational chemistry, whichprovidesmethodsofcomputersimulation,or,ina sense, a virtual laboratory. Most fundamentally, computational chemistry considers matter as a many-particle system of two different constituents: atomic nuclei and electrons, interacting through the electromagnetic force. Modeling this chemical reality by computer simulations receives its justification and merit from the full understanding of the underlying physical principles. In fact, already in 1929, Dirac realized that the underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble. 2 Indeed, with the Dirac equation, and certainly today with quantum electrodynamics, well-elaborated theories, accounting for both quantum mechanics and relativity, have been developed. For many chemical systems, these levels of elaboration are not even necessary and one can restrict the description to classical electrodynamics and non-relativistic quantum mechanics, that is, to the chrödinger equation. The second part of the quotation might, due to the ever-increasing processing power of computers, be seen in a different light today. Although certain approximations are indeed required in order to transform an exact quantum mechanical equation of motion from its abstract form into actually tractable expressions (these fundamental approximations include in particular the Born-Oppenheimer approximation, which facilitates a decoupling of the nuclear and electronic motions), the enormous amount of calculation required for solving the resulting expressions has become an increasingly feasible task. As of today, computational chemistry techniques which are solely based on first principles of physics and fundamental approximations, so-called ab initio methods, are applicable to complex systems of microscopic size, such as molecular systems. The background of these foundations of computational chemistry, covering in particular the separation of nuclear and electronic motions, their decoupling by means of the Born-Oppenheimer approximation, as well as ab initio electronic structure theory is discussed in Chapter 2.

18 Introduction 17 One of the most important macroscopic properties of organic photovoltaic cells is their ability to produce an electrical current. The prerequisite for such a current flow is the migration of charge carriers through the organic material. This process is denoted as charge transport and is characterized in terms of the charge carrier mobility. This quantity, which is associated with the average velocity of the charge carriers, can be experimentally measured for a given sample of the organic material. However, when it comes to optimizing the material for an improved mobility, one faces the problem of missing structure-property relationships linking the mobility to the constituent molecules. With the aim of closing this knowledge gap, this work applies computer simulations to reconstruct the macroscopic process of charge transport based on its microscopic origins. These charge transport simulations are based on a model including the following three key components. First, a sufficiently large, but at the same time atomistically resolved model of the organic material, a so-called morphology. To generate such a large-scale material morphology, ab initio methods need to be supplemented by computational chemistry techniques operating on a higher level of approximation. This is achieved by molecular dynamics simulations, where the motion of atoms is governed by classical Newtonian mechanics, calibrated according to ab initio methods. The second ingredient is a model of charge transport which describes charge carrier migration within the morphology as a sequence of charge hops between individual molecules. These microscopic processes of charge movement are referred to as charge transfer reactions and the quantity characterizing their efficiency is the charge transfer rate. uch a rate is influenced by several factors: the electronic structure of the two individual molecules, their relative positions and orientations, but also their environment of surrounding molecules. The third component is an appropriate model of charge transfer, which translates these dependences into a set of tangible parameters, which are accessible by methods of computational chemistry. In some cases, the applicability of ab initio methods may be limited by the large number of molecular pairs for which charge transfer parameters need to be evaluated. The parametrization can then be assisted by semiempirical methods, which are still based on the quantum mechanical level of description, but incorporate certain empirical data to accelerate the computation. Altogether, the simulation of charge transport invokes a hierarchy of methods to scan all the required length and time scales with amanageablecomputationaleffort. Whilethetheoryofchargetransferinmolecular systems is treated in Chapter 3, the complete methodology of charge transport in organic solids is presented in Chapter 4.

19 18 Introduction C60 dcvnt NC NC n Chemical structures of buckminster fullerene (C 60)anddicyanovinyl-substituted oligothiophenes (dcvnt). At present, the most successful materials for building organic photovoltaic cells are compositions of buckminster fullerene (C 60 ) and the novel class of dicyanovinyl-substituted oligothiophenes (dcvnt) or its derivatives. In fact, based on these compounds, an ongoing series of world record power conversion efficiencies has been achieved between 2009 and 2013 by Heliatek GmbH. 3 In addition to these proprietary cells, a wide range of related devices has been published by the collaborating groups of Bäuerle at the Institute of Organic Chemistry II and Advanced Materials in Ulm, Germany, and Leo at the Institute for Applied Photo Physics in Dresden, Germany. In this work, charge transport is studied for a variety of dcvnt material morphologies associated with these devices. A simulation study on single crystals and a further one on thin films are presented in Chapter 5. For single crystals, charge transport is compared for asetoffoursystems:theterthiopheneandquaterthiophene,dcv3tanddcv4t,aswell as two methylated derivatives, dcv3t-m and dcv4t-m. For thin films, charge transport is examined in systems of the bare and methylated quaterthiophenes, dcv4t and dcv4t-m. A simulation study on amorphous and smectic systems of the compound series of thiophene to sexithiophene, dcv1t to dcv6t, is presented in Chapter 6. Parts of the methodology and the results reported in this work (ections ) are the subject of prior publications, listed on Page 19. These studies are presented here in significantly more detail. The background and the methodology (ections ) have been developed based on the textbooks and review articles provided at the beginning of the respective discussions. All chapters of this work employ a notation for symbols summarized on Page Lewis, N.. Toward Cost-Effective olar Energy Use. cience 315 (2007), Dirac, P. A.M. Quantum Mechanics of Many-Electron ystems. Proc. R. oc. London, er. A 123 (1929), Le éguillon, T., and Pfeiffer, M. Efficiency Development. Heliatek,

20 Related Publications In the course of these doctoral studies, the following journal articles were published. The publications on charge transport [1 4] are related to the methodology for charge transport simulations (ection 4.2) and to the results on single crystals (ection 5.1), thin films (ection 5.2), and amorphous/smectic systems (ection 6.1) of dicyanovinylsubstituted oligothiophenes. The independent study on proton transport [5] is not addressed in this thesis. [1] chrader, M., Fitzner, R., Hein, M., Elschner, C., Baumeier, B., Leo, K., Riede, M., Bäuerle, P., and Andrienko, D. Comparative tudy of Microscopic Charge Dynamics in Crystalline Acceptor-ubstituted Oligothiophenes. J. Am. Chem. oc. 134 (2012), [2] chrader, M., Körner,C.,Elschner,C.,andAndrienko,D.ChargeTransport in Amorphous and mectic Mesophases of Dicyanovinyl-ubstituted Oligothiophenes. J. Mater. Chem. 22 (2012), [3] Elschner, C., chrader, M., Fitzner, R., Levin, A. A., Bäuerle, P., Andrienko, D., Leo, K., and Riede, M. Molecular Ordering and Charge Transport in a Dicyanovinyl-ubstituted Quaterthiophene Thin Film. RC Adv. 3 (2013), [4] Rühle, V., Lukyanov, A., May, F., chrader, M., Vehoff,T.,Kirkpatrick,J.,Baumeier, B., and Andrienko, D. Microscopic imulations of Charge Transport in Disordered Organic emiconductors. J. Chem. Theory Comput. 7 (2011), [5] Wehmeyer, C., chrader, M., Andrienko, D., and ebastiani, D. Water-Free Proton Conduction in Hexakis-(p-phosphonatophenyl)benzene Nano-Channels. J. Phys. Chem. C 117 (2013),

21 20 Related Publications

22 Notation Indices n, m electrons a, b nuclei / atoms i, j molecules / sites α, β electronic states +, / i, f adiabatic / diabatic electronic states η, ϑ nuclear states n, m molecular orbitals σ, τ atomic orbitals µ, ν general vector and matrix elements Entity and Pair Properties r n / r electronic coordinates / multi-index p n σ n R a / R R ab P a V a M a z a p a q c a / q n a α c a / α n a Q / P / ω vib electronic momenta electronic spin nuclear coordinates / multi-index nuclear separations nuclear momenta nuclear velocities nuclear masses atomic numbers atomic dipole moments atomic partial charges of charged / neutral site atomic polarizabilities of charged / neutral site reaction coordinate / momentum / eigenfrequency

23 22 Notation r i r ij p i α n i / α c i U n i / U c i U n i / U c i Wi n E int i E elstat i / W c i / E int ij E i / E ij λ ij / λ i J ij ω ij c ij / E elstat ij site coordinates site separations site occupation probabilities site polarizability tensor in neutral / charged state internal site energies in neutral / charged state internal site energies in neutral / charged state, opposite geometry electrostatic site energies in neutral / charged state internal site energies / differences electrostatic site energies / differences site energies / differences reorganization energies / site contributions transfer integrals charge transfer rates edge currents ystem Properties and Observables ρ electron density n number of electrons a number of nuclei / atoms d mass density T temperature F external electric field t time of nuclear motion τ time of charge carrier motion λ / Λ diagonal / off-diagonal dynamic disorder σ / Σ diagonal / off-diagonal static (energetic / electronic) disorder σ / σ eff energetic disorder of neighbor list / reduced neighbor list C E Q µν / Q D µν / D site energy correlation function nematic order tensor charge carrier diffusion tensor µ µν / µ charge carrier mobility tensor

24 Notation 23 Operators Ĥ T nuc V nuc nuc T el / ˆt el (r n ) V el el / ˆv el el (r n, r m ) V nuc el / ˆv nuc el (r n ) Ĥ el (R) Ĥnuc(R) α Θ nuc(r) αβ Ĥ nuc(q) ± / Ĥnuc(Q) i,f Θ + nuc(q) / J if (Q) / J ij Ĥ 1 / Ĥ 2 Ĥ 3 / Ĥ 4 Ĥ hf / Ĥ ks / Ĥ ˆv h (r 1 ) / ˆv x (r 1 ) / ˆv xc (r 1 ) Hamiltonian operator nuclear kinetic energy nuclear-nuclear interaction electronic kinetic energy / for electron n electronic-electronic interaction / for electrons n, m nuclear-electronic interaction / for electron n electronic Hamiltonian operator nuclear Hamiltonian operator nonadiabatic coupling (nonadiabaticity operator) adiabatic / diabatic nuclear Hamiltonian operator nonadiabatic / nondiabatic coupling (transfer integral) equilibrium / phononic Hamiltonian local / nonlocal electron-phonon coupling Hartree-Fock / Kohn-ham / one-particle operator Hartree / exchange / exchange-correlation operator Wave Functions / Eigenstates Ψ(R, r, t) χ αη (R, t) ψ α (r, R) ψ ± (Q) / ψ i,f / ψ i ϕ n i (r 1) / ϕ n i / ϕf i φ τ i (r 1) / φ τ i total wave function nuclear wave functions electronic wave functions adiabatic / diabatic electronic states molecular orbitals / frontier orbitals atomic orbitals Eigenvalues and Potential Energy urfaces E α el (R) U α (R) E el ± (Q) / Ei,f el U ± (Q) / U i,f (Q) ε n i electronic eigenvalues potential energy surface adiabatic / diabatic electronic eigenvalues adiabatic / diabatic potential energy surface molecular orbital energies

25 24 Notation Functionals, Matrices, and Tensors E x [ρ] / E c [ρ] / E xc [ρ] exchange / correlation / exchange-correlation functional Exc b3lyp [ρ] b3lyp exchange-correlation functional T ab / T µ ab / T µν ab multipole interaction tensor / first / second derivative H στ / H one-electron Hamiltonian matrix Hzindo στ / H zindo diag(ε n ) / E C τn / C στ / zindo Hamiltonian matrix molecular orbital energy matrix atomic orbital matrix atomic orbital overlap matrix Photovoltaic Cell Properties and olar Parameters η pce η eqe η ff j / j sc / j mp V / V oc / V mp P solar Φ solar (E) power conversion efficiency external quantum efficiency fill factor current density / at short circuit / for maximum power voltage / at open circuit / for maximum power solar power density solar spectral photon flux density Physical Constants m el e ε 0 c ħ k B electron mass elementary charge vacuum permittivity speed of light Planck constant Boltzmann constant

26 Part I. Background

27

28 Chapter 1. Organic Photovoltaic Cells In contrast to their silicon-based inorganic counterparts, organic photovoltaic cells are manufactured from organic, i.e., carbon-based molecules. Depending on the molecular weight, there is a common classification into organic solar cells produced from polymers and from small molecules. This distinction refers to the processing techniques used for preparing the desired layers of organic molecular solids: while polymers are dissolved in solutions, which are solidified by solvent removal techniques, small molecules are mostly processed by vacuum evaporation or sublimation and subsequent material deposition. However, both types of cells share the same working principle for the photovoltaic power conversion. A third type of functionally different organic solar cells, which is not related to this work, is the class of dye-sensitized solar cells. The following discussion opens with a qualitative insight into the electronic structure of organic molecular solids (ection 1.1). ince organic solids possess relatively weak cohesive intermolecular interactions, their electronic structure can be regarded as a perturbed one of its constituent molecules. Molecules of particular interest are those which comprise π-conjugated systems, since they can enable the desired semiconducting properties of the organic solid. Then, the focus is directed to organic photovoltaic cells and their working principle for the conversion of solar radiation into electrical power (ection 1.2). The power conversion is based on four optical and electronic processes: optical absorption yielding an exciton (a bound electron-hole pair), exciton diffusion, exciton dissociation into free charge carriers, and charge transport towards the electrodes. Finally, the most important metric of a photovoltaic cell the power conversion efficiency is introduced (ection 1.3). After briefly reflecting on the theoretical upper limits for the efficiency, the currently achieved values and challenges for further improvements are discussed. Among the main challenges are the improvement of the light harvesting, the active layer morphologies, but also the fundamental understanding of how the efficiency is linked to the properties of the constituent molecules.

29 28 Chapter 1. Organic Photovoltaic Cells 1.1. Electronic tructure of Organic olids In general, organic solids [6 10] are solid-state materials which are composed of molecules falling within the scope of organic chemistry. These organic molecules are predominantly composed of carbon atoms and exist, due to the versatile bonding capabilities of carbon, in a myriad of architectures. This diversity of compounds is reflected in a wide spectrum of observed solid state order, ranging from the perfect crystalline to the amorphous phase and covering many intermediate forms, such as polycrystalline, semicrystalline, or mesomorphic phases. While most organic solids are insulators, the field of organic photovoltaics is primarily concerned with the subclass of materials acting as (semi)conductors. Thesematerials,capableofcarryinganelectriccurrent, are generally composed of molecules which have electrons delocalized over larger, socalled conjugated systems of the molecular skeleton. Important building blocks for such conjugated systems are aromatic hydrocarbons, such as the polyacenes, i.e., benzene, naphthalene, anthracene, etc., or heterocyclic compounds, such as thiophene, furane, pyrrole, etc., which are depicted in Figure 1.1. To understand the origin of electronic delocalization in conjugated molecules, one can start from the familiar viewpoint of independent electrons, described by individual wave functions, i.e., molecular orbitals. In addition, these molecular orbitals shall be composed as linear combinations of atomic orbitals (mo-lcao). In fact, these concepts constitute electronic structure theories (ection 2.2), which enable one to quantitatively derive the right linear combinations for composing the molecular orbitals. One can then verify that there are indeed delocalized orbitals. For a qualitative understanding, however, the notion of valence bond theory may be illustrative. There, pairs of overlapping atomic valence orbitals give rise to bonding molecular orbitals, i.e., shared electron pairs lead to covalent bonds. This simple picture is accompanied by the Benzene Naphthalene Anthracene Thiophene O Furan H N Pyrrole Figure 1.1. election of basic conjugated organic molecules, acting as building blocks for small molecules, oligomers, or polymers employed in organic electronic devices.

30 1.1. Electronic tructure of Organic olids 29 (a) y x p z -orbital sp2-orbital π-bonding σ-bonding (b) 6p z 2sp2 lumo homo σ π π σ Figure 1.2. (a) Atomic valence orbitals of the carbon atoms in benzene. Overlapping sp 2 -orbitals give rise to σ-bonding, while overlapping p z -orbitals lead to π-bonding. (b) plitting of the energy levels of two atomic sp 2 -orbitals leading to bonding and antibonding molecular σ- and σ -orbitals as well as of six p z -orbitals yielding π- and π -orbitals. Adapted from Reference [11]. idea of hybridization: it allows one to transform each basis of standard (hydrogen-like) atomic orbitals, by intuitive linear combinations, to equivalent bases of so-called hybrid atomic orbitals. Considering carbon, the ground state electron configuration in terms of hydrogen-like atomic orbitals reads 1s 2 2s 2 2p x 2p y with two valence electrons. Hybrid atomic orbitals are, however, derived from the excited electron configuration 1s 2 2s 2p x 2p y 2p z with four valence electrons. This is because the energy expenditure for the excitation is more than compensated by the formation of two additional bonds. Asimplelinearcombinationofthe2s-,2p x -, and 2p y -orbitals leads to three hybrid sp 2 -orbitals, which lie in the xy-plane at angles of 120. The p z -orbital remains unchanged and is perpendicular to the xy-plane. Using the example of benzene, containing six carbon atoms in a hexagonal arrangement, these orbitals are illustrated in Figure 1.2 a. Now, pairs of overlapping atomic sp 2 -orbitals in the xy-plane give rise to molecular σ-orbitals,which are localized between the respective pairs of nuclei. Figure 1.2 b shows how the energy levels of an overlapping pair of sp 2 -orbitals are split into an energetically lower level, corresponding to a bonding σ-orbital, which is doubly occupied, and a higher level, corresponding to an antibonding σ -orbital, which is vacant. Due to the strong overlap of sp 2 -orbitals, the energy splitting and the resulting energetic advantage is large, and therefore the σ-bonding a strong effect. In total, the molecular backbone of the benzene molecule involves twelve sp 2 -orbitals forming the hexagon and a further six linking the hydrogens. The remaining six atomic p z -orbitals are also overlapping, namely above and below the xy-plane, which gives rise to three bonding molecular π-orbitals, which are doubly occupied, and three antibonding π -orbitals, which are empty. Obviously, the three π-orbitals cannot be localized between three pairs of nuclei, since all six pairs of nuclei are equivalent by the molecular symmetry. In fact, the π-orbitals are instead delocalized over the molecular skeleton. As the overlap of the p z -orbitals is weak, their energy splitting is small and

31 30 Chapter 1. Organic Photovoltaic Cells lumo P E dos E E E dos homo P Molecule Ordered olid Disordered olid Figure 1.3. Electronic structure of organic solids: one observes a general shift as well as a slight broadening of the molecular energy levels. While ordered solids show narrow energy bands, disordered solids often exhibit a Gaussian distributed density of states (dos). Adapted from Reference [8]. therefore the mechanism of π-bonding comparatively weak. In many conjugated organic molecules, the highest occupied molecular orbital (homo) is a π-orbital, while the lowest unoccupied molecular orbital (lumo) is a π -orbital. The electronic structure discussed so far refers to isolated organic molecules, as they are encountered in the gas phase. In an organic solid, formed upon condensation of the molecules, the electronic structure changes, since molecules interact with each other. The interaction between molecules, causing their cohesion, is dominated by the van der Waals interaction, providedthemoleculesareneutralandarenotforming ionic bonds. Van der Waals interactions result from fluctuations in the molecular charge distributions: such fluctuating dipole moments polarize adjacent molecules, leading to an induced dipole-dipole attraction. ince these intermolecular interactions are much weaker than the strong covalent binding forces within the molecules, the molecular properties remain largely intact in an organic solid. Thus, the electronic structure of the solid is only a moderately altered one of a free molecule. The main differences are illustrated in Figure 1.3. First, one observes a general shift of the energy levels due to the polarizable environment. The homo and lumo energies, i.e., the ionization potential E ip and electron affinity E ea (in Koopman s approximation), are displaced by the polarization energies P ea and P ip,respectively.therefore,inthe solid, the difference between the ionization potential and electron affinity is usually lowered. econd, the energy levels in the solid are slightly broadened due to the weak overlap of the molecular orbitals. In the case of ordered solids, such as crystalline phases at low temperatures, narrow energy bands can emerge. In analogy to inorganic materials, these bands are sometimes referred to as the valence and conduction bands of the organic solid and the region in between as the band gap. In the case of disordered solids, such as amorphous or mesomorphic phases, the density of states (dos) is often described by Gaussian distributions. Then, the distribution tails extend into the band gap and the band edges are no longer clearly defined.

32 1.2. Conversion of olar Radiation into Electric Power Conversion of olar Radiation into Electric Power Organic photovoltaic cells [12 18] make use of organic solids to convert solar photons into electric voltage and current. A major similarity to inorganic cells, which are mostly based on silicon, is that the photoactive organic materials are semiconductors. Therefore, the photovoltaic effect can be exploited for promoting electrons across the band gap, whichpreventstherapiddecaybacktothegroundstatebyaseriesofphonons, as would occur without the gap. A key difference in organic semiconductors, on the other hand, is that a promoted electron is not free,but instead electrostatically bound to the remaining hole. The bound electron-hole pair is denoted as exciton and therefore organic solar cells sometimes as excitonic solar cells. For the separation of excitons, the most common concept is to use a junction between two different organic semiconductors, which is referred to as a heterojunction. Thisdevicedesignwasfirst proposed in 1986 by Tang in the much-cited Reference [19]. The basic working principle of a heterojunction solar cell involves four optical and electronic processes, which are illustrated in Figure 1.4: optical absorption yielding an exciton, exciton diffusion to the heterojunction, exciton dissociation into free charge carriers, and charge transport to the electrodes. First, upon optical absorption, asolarphotonpromotesanelectronwithinoneof the two different organic semiconductors across the band gap. This is possible since, due to the π-conjugation, organic semiconductors exhibit relatively low band gaps, roughly between 1 and 4 ev, which lies within the spectrum of the solar radiation received on earth. After the photoexcitation, the system rapidly relaxes to the band edges, i.e., dissipates the energy exceeding the band gap via a series of phonons as heat, and finally forms an exciton. The exciton binding energy, that is, the electrostatic interaction energy between the electron and hole, is of the order of 0.1 to 1 ev in organic materials, which is significantly higher than thermal energy at room temperature. As aconsequence,theelectronandholearenot free. Thisstrongelectrostaticattractionis a result of the low dielectric constants, i.e., the weak electrostatic screening of organic materials. In inorganic semiconductors, in comparison, exciton binding energies are of the order of 10 3 ev and photoexcited electrons and holes are free at room temperature. Compared to inorganic materials, organic semiconductors also have significantly higher absorption coefficients. As a consequence, organic photoactive layers can be much thinner. A thickness of the order of 100 nm is usually sufficient to absorb most incident photons whose energy bridges the band gap.

33 32 Chapter 1. Organic Photovoltaic Cells Figure 1.4. Optical and electronic processes taking place in an organic photovoltaic cell: first, optical absorption yielding an exciton; second, exciton diffusion to the heterojunction; third, exciton dissociation into free charge carriers; and fourth, charge transport to the electrodes. econd,exciton diffusion to the heterojunction is required. As an exciton is a neutral quasiparticle, which is not affected by any electric fields, its migration is a purely diffusive process. The diffusion length is determined by the finite lifetime of the exciton and is of the order of 10 nm. Within this length scale, the exciton, traveling within one of the two semiconductors, must reach the interface to the other one, otherwise it is lost due to radiative recombination. In planar heterojunction architectures, where the two semiconductors are arranged in two layers on top of each other, the exciton diffusion length obviously requires thinner layers than are needed for efficient photon absorption (100 nm). It is therefore necessary to find a compromise for the layer thickness. In order to avoid such a trade-off, one can employ bulk heterojunction architectures [20], where the two semiconductors are mixed to an interpenetrating network, as sketched in Figure 1.4. This design allows the interface area to be increased, while at the same time tuning the layer thickness for optimal absorption. Third, exciton dissociation can take place once the exciton has reached the heterojunction of the two semiconductors. The rationale behind this heterojunction, as introduced by Tang, is to provide appropriate energetic steps between the ionization potentials and electron affinities, aligned such as to overcome the exciton binding energy and therefore to facilitate the separation of the electron-hole pair. Figure 1.5 a depicts the required level alignment of the two semiconductors, which are henceforth referred to as the electron donor and acceptor,respectively.theenergydifferencebetweenthe ionization potential of the donor and the electron affinity of the acceptor, E d ip E a ea, must be more than the binding energy lower than the band gap of either material, i.e., E d ip E d ea and E a ip E a ea, provided excitons are generated in both materials. The illustration shows the case where the exciton is formed within the donor: since the energetic step in the electron affinity at the donor-acceptor heterojunction exceeds the binding

34 1.2. Conversion of olar Radiation into Electric Power 33 (a) (b) E E Cathode V Cathode Anode E Donor E Anode Acceptor Figure 1.5. (a) Energy level alignment of a donor-acceptor heterojunction solar cell required to facilitate exciton dissociation into charge carriers. (b) chematic energy diagram of the solar cell under operating conditions leading to drift currents of charge carriers towards the electrodes. energy, the separation of the electron and hole is an energetically favorable process. Therefore, the electron can be transferred from the donor to the acceptor, while the hole remains on the donor. Conversely, if the exciton is formed within the acceptor, the hole can be transferred from the acceptor to the donor, while the electron remains on the acceptor. Fourth, charge transport of the free electron and hole towards their respective electrodes occurs as a result of diffusion and drift [12, 21]. While charge diffusion, similar to the migration of excitons, occurs independently of electric fields, drift currents of the charge carriers are a result of the electric potential gradient inherent in the device. As illustrated in Figure 1.5 b, this potential gradient arises once the anode and cathode are either short-circuited, as indicated by the dashed line, or connected to an external circuit with a voltage drop V. The higher the voltage drop across the external circuit, the lower the internal potential gradient and thus the drift currents. If the voltage drop nearly cancels the internal potential gradient, the migration of electrons and holes is dominated by diffusion currents. Finally, the charge carriers are collected at their respective electrodes, i.e., the electron at the cathode and the hole at the anode. The electrodes, as conductors, are solely characterized by their Fermi levels, or their work functions. In an idealized model, the work function of the cathode matches the electron affinity of the electron acceptor, while the work function of the anode fits to the ionization potential of the donor. In practice, the cathode is often manufactured from aluminum, while the common choice for the anode is indium tin oxide (ito), which is not only conductive, but also transparent for the incident light.

35 34 Chapter 1. Organic Photovoltaic Cells 1.3. Power Conversion Efficiency The power conversion efficiency of any photovoltaic cell [22 24] depends on the device characteristics when operating in an electric circuit. In principle, both inorganic and organic devices generate a photocurrent under illumination, while they exhibit rectifying properties of a diode in the dark. This similarity is because the energetic step in an organic device, due to the heterojunction of two different semiconductors, is essentially similar to the step arising in an inorganic device upon contacting a p-and a n-doped material to a pn-homojunction [14]. As a consequence, any ideal solar cell can be modeled by an equivalent circuit consisting of a currentsourcein parallelwith a diode, as illustrated in Figure 1.6 a. The current-voltage characteristic of a solar cell exposed to light thus corresponds to a shifted diode characteristic, as seen in Figure 1.6 b. If the electrodes of the solar cell are connected, that is, R = 0, no voltage between them can be established and the cell delivers the short-circuit current density j sc.(notethat the current density j is used instead of the current I, since the photocurrent is ideally proportional to the illuminated area.) Conversely, if the electrodes are isolated, that is, R =, no current can flow and the cell develops the open-circuit voltage V oc. This case corresponds to Figure 1.5 a. For any intermediate applied resistance R, the cell generates a voltage V and a current density j = j(v), according to the currentvoltage characteristic, such that R = V/I. This general case corresponds to Figure 1.5b. At any point on the current-voltage characteristic, the electric power density supplied by the solar cell is given by the product of j and V. Thepointmaximizingthisproduct j j dark (a) j j j j (b) Maximum Power V j R V dark V V light V Figure 1.6. (a) Equivalent circuit of an ideal solar cell consisting of a current source in parallel with a diode. Under illumination, the cell generates a current density j and voltage V. (b) Current-voltage characteristic j = j(v). In the light,a shifted characteristic of an ideal diode is encountered.in the dark, the ideal diode characteristic is obtained when a voltage is applied.

36 1.3. Power Conversion Efficiency 35 determines the current density and voltage for maximum power, which are denoted as j mp and V mp.undertheseoperatingconditions,theratioofthemaximumelectric power density and the incident solar power density P solar defines the Power Conversion Efficiency η pce = j mpv mp P solar = η ff j sc V oc P solar. (1.1) Here, the fill factor η ff is introduced to easily reflect the shape of the current-voltage characteristic. It is defined as the quotient of the two rectangular areas in Figure 1.6 b: η ff = j mpv mp j sc V oc. (1.2) Upper Limits for the Efficiency As pointed out by hockley and Queisser [25], atheoreticallimitingefficiencyforanidealsolarcellcanbedeterminedbythreebasic assumptions. First, the device exhibits perfect absorption, i.e., each incident photon produces an exciton, provided the photon energy bridges the optical gap E gap of the absorbing semiconductor (i.e., the electronic gap minus the exciton binding energy). econd, there is no internal device resistance and each electron-hole pair is instantaneously collected at the electrodes if they are short-circuited. With these assumptions, the ideal short-circuit current density simply equals the elementary charge times the number of absorbed photons per time, which can be written as j sc = e E gap deφ solar (E). (1.3) Here, Φ solar is the solar photon flux density in spectral distribution and the lower integration limit reflects the minimum energy of absorbed photons, as shown in Figure 1.8. The third assumption refers to the case where an external resistance is applied to the electrodes and hence charge carriers can no longer be collected instantaneously. In this case, an inevitable process, occurring in addition to absorption, is the spontaneous emission of photons as a result of radiative recombination of electron-hole pairs. By relating generation and recombination rates according to the principle of detailed balance, the current-voltage characteristic of the ideal solar cell (shown in Figure 1.6 b) can be parametrized. With the current-voltage function j(v) at hand, the ideal open-

37 36 Chapter 1. Organic Photovoltaic Cells ev η Figure 1.7. Upper limit of the power conversion efficiency η pce as a function of the optical gap E gap of the absorbing semiconductor and the voltage loss due to exciton dissociation. The values along the abscissa, i.e., for zero voltage loss ( = 0), correspond to the hockley-queisser limit for inorganic photovoltaic cells. Adapted from Reference [26]. E gap ev circuit voltage is defined by the point of vanishing current, which is V(j = 0) 1 e E gap. Although this voltage represents a strict limit for the ideal solar cell, irrespective of its inorganic or organic nature, one can argue more precisely for organic cells. In fact, exciton dissociation in organic devices involves intermediate charge-transfer or chargeseparated states and thus entails further inevitable energy losses. To account for these inherent losses, conceivable as the driving force for exciton dissociation, the considerations of hockley and Queisser can be extended by a voltage loss parameter [26]: V oc = V(j = 0) 1 e. (1.4) With the short-circuit current (1.3) and the open-circuit voltage (1.4), the fill factor (1.2) is, of course, determined and one realizes that the power conversion efficiency (1.1) becomes a function of the optical gap E gap and the voltage loss. This function, shown in Figure 1.7, indicates that the maximum efficiency for a given is achieved for some intermediate gap E gap.thisisbecausetheshort-circuitcurrentgoestozeroforlarge gaps, while the open-circuit voltage vanishes for small gaps. For zero voltage loss, i.e., along the abscissa, the hockley-queisser limit for inorganic solar cells is reproduced, which is 33.7% at an optical gap of 1.34 ev. If the voltage loss is = 0.2 ev, organic cells can theoretically achieve efficiencies slightly above 25% for optical gaps between 1.1 and 1.7 ev [26]. It should be mentioned that these upper limits apply to solar cells with a single absorbing semiconductor and can be surpassed by tandem cells.

38 1.3. Power Conversion Efficiency 37 Quantum Efficiency (η ) η Φ solar C60 dcv5t-m Photon Energy (E ev) Photon Wavelength (λ nm) s 1 m 2 nm 1 Φsolar Figure 1.8. pectral distribution of the terrestrial solar photon flux density Φ solar and external quantum efficiency η eqe of an organic photovoltaic cell prepared from a dcv5t-m donor and a C 60 acceptor [28]. Challenges for Efficiency Improvements The highest reported power conversion efficiency of organic solar cells is currently 12%, reached by tandem cells developed by Heliatek GmbH [27]. The absorbers employed are proprietary derivatives of dicyanovinyl-substituted oligothiophenes (dcvnt), acting as donor materials, while fullerene C 60 is the acceptor. For this class of compounds, depicted on Page 18, remarkable efficiencies have also been published for single-absorber devices, the highest of which is 6.9%, reported for a methylated quinquethiophene (dcv5t-m) [28]. In view of the first donor-acceptor heterojunction device by Tang with 1% efficiency, these achievements demonstrate the significant progress in the field of organic photovoltaics. Now that organic solar cells are close to commercialization, the scientific community is more than ever challenged to further advance the power conversion efficiencies in the direction of the theoretical limiting efficiencies. This goal requires improvements in all parameters entering the efficiency: the short-circuit current density, the open-circuit voltage, and the fill factor. The short-circuit current density of a real solar cell is, contrary to the ideal form (1.3), determined by the external quantum efficiency η eqe, which is defined as the ratio of collected charge carriers and incident photons of a given energy E: j sc = e deη eqe (E) Φ solar (E). (1.5) Figure 1.8 depicts, as an example, the quantum efficiency of the mentioned solar cell prepared from dcv5t-m and C 60. The path to improved quantum efficiencies is to

39 38 Chapter 1. Organic Photovoltaic Cells minimize the losses of photons upon optical absorption, of excitons during exciton diffusion and dissociation, and of charge carriers during charge transport; in short: losses in all processes involved in the conversion of radiation into electricity (ection 1.2). Optical absorption with little loss critically depends on an optimized optical gap of the absorbing semiconductor. In fact, most organic semiconductors suffer from overly large band gaps, resulting in a lack of light harvesting in the low-energy, i.e., the high-wavelength region of the solar spectrum. Fullerene C 60,forinstance,can only harvest a tiny fraction of solar photons, as can be seen in Figure 1.8. Fortunately, in the last decade, great progress has been made in synthesizing absorbers with lower optical gaps. The class of dcvntcompoundsisoneexample,wherethisisrealizedby attaching electron-withdrawing dicyanovinyl moieties to conjugated oligothiophene cores. As a result, the light harvesting reaches the visible range of the solar spectrum, as seen in Figure 1.8. However, synthetic chemistry still faces many challenges, in particular in engineering materials extending the absorption to the infrared region. Exciton diffusion and dissociation with little loss requires both efficient exciton transport and an optimized solar cell architecture. In planar heterojunction layouts, the thickness of the active layer must satisfy a compromise between optimal extinction of light and lossless exciton diffusion to the interface. In bulk heterojunction layouts, the layer thickness can be optimized for light extinction, however, the scale of the donor-acceptor phase separation in the morphology must meet a trade-off between maximum interfacial area and continuous percolation paths to the electrodes. While optimal scales of phase separation have been determined theoretically [29, 30], it remains a major challenge for the field of material processing to control the bulk heterojunction morphology during layer preparation in the desired manner. Charge transport with little loss requires, apart from proper percolation paths, high electron and hole mobilities in the acceptor and donor domains, respectively. This is particularly important since mobilities in organic semiconductors can vary by many orders of magnitude. Among the materials with the highest electron mobilities are fullerenes and its derivatives, which are therefore the most common candidates as acceptors. However, when engineering associated donors, these need to be primarily tuned for optical absorption, making it problematic to simultaneously ensure optimal hole mobilities. Although increased hole mobilities are often highly desirable, there are few systematic strategies for this purpose. A major obstacle is the lack of structureproperty relationships, linking chemical structures to mobilities. Closing this knowledge gap is an important challenge for computational chemistry.

40 1.3. Power Conversion Efficiency 39 The open-circuit voltage of a real solar cell is, contrary to that of the ideal one (1.4), usually approximated by the difference between the ionization potential of the donor and the electron affinity of the acceptor, since these band edges are thought to characterize the hole and electron energies after the exciton dissociation: V oc = 1 e ( Ed ip E a ea ). (1.6) However, several sometimes insufficiently understood effects lead to an empirical voltage loss, which is typically within 0.3 to 0.7 ev [31 33]. An important part of this loss stems from the energy differences between the band edges and the actual charge transport levels. In fact, band transport is seldom seen in organic materials. This is because effects of disorder often influence the electronic structure of organic solids (ection 1.1). Most particularly, disordered solids, such as amorphous or mesomorphic materials, are typically characterized by Gaussian distributed densities of states rather than by sharp energy bands. Then, hopping transport is observed, where charge transport levels are located in the distribution tails. But even organic crystals, if heated to room temperature, tend to exhibit hopping or so-called polaronic transport rather than band transport. As a consequence, improving the open-circuit voltage can benefit from an in-depth understanding of the nature of charge transport, as can be provided by computational chemistry. The fill factor of a real solar cell is, as opposed to the ideal cell, no longer a function of the short-circuit current and the open-circuit voltage, since the current-voltage characteristic has no predefined functional form. Instead, the shape of the current-voltage characteristic is negatively affected by mainly two parasitic effects [14]. First, leakage currents flowing in the opposite direction of the photocurrent lead to a reduction of the maximum power current j mp.thiseffectisrelatedtothequalityofthesolarcell and may be reduced by careful device preparation. econd, the finite resistance of the device leads to a voltage drop under operation, i.e., a reduction of the maximum power voltage V mp. This effect is mostly due to the finite conductivity of the semiconducting materials and can thus be reduced by improving the charge carrier mobilities. In light of these challenges for improving power conversion efficiencies, this work is devoted to the field of charge transport using methods of computational chemistry. While the characterization of morphological disorder allows the nature, or regime of charge transport, to be specified, the simulation of charge carrier dynamics delivers the missing links between chemical structures and mobilities.

41 40 Chapter 1. Organic Photovoltaic Cells

42 Chapter 2. Foundations of Computational Chemistry Computational chemistry has its foundation in viewing the structure of matter as a many-body system of point-like nuclei and electrons. These constituents carry an electrical charge and thus interact through the electromagnetic force, expressed as an instantaneous field according to classical electrodynamics. The quantum nature of the constituents is taken into account by the equation of motion of non-relativistic quantum mechanics, i.e., the chrödinger equation, complemented by spin if required. With this powerful framework, virtually all chemical and related physical system properties become, in principle, accessible by solving the chrödinger equation. In what follows, the initial step towards a solution is recalled: the separation of nuclear and electronic motions (ection 2.1). This approach is motivated by the large mass ratio of nuclei and electrons, implying that nuclear positions remain nearly fixed while the electronic motion takes place. The separate chrödinger equation for the electrons thus depends on the nuclear degrees of freedom only as a fixed parameter. If this electronic chrödinger equation is solved, i.e., the electronic states and energies are obtained, and that for a large set of fixed nuclear configurations, the separate nuclear equation of motion can be parametrized. To solve the nuclear equation of motion, one can often, but not always, neglect electronic state transitions, which is known as the adiabatic or Born-Oppenheimer approximation. In the classical limit, these approximations lead to the Newtonian equation of motion for the nuclei. ubsequently, an insight into electronic structure theory is offered (ection 2.2). This field focuses on finding solutions of the electronic chrödinger equation at fixed nuclear positions. Two formalisms are discussed: Hartree-Fock and Kohn-ham density functional theory. Both theories reduce the many-electron chrödinger equation to a set of one-electron equations governed by an effective one-electron operator.

43 42 Chapter 2. Foundations of Computational Chemistry 2.1. eparation of Nuclear and Electronic Motions A chemical system can be regarded as an assembly of electrostatically interacting atomic nuclei and electrons, governed by quantum mechanical dynamics. The full Hamiltonian operator of such a system depends on the cartesian coordinate and conjugate momentum operators of the nuclei, R a and P a,aswellasonthoseoftheelectrons, r n and p n [34 37]. It can be written as a sum Ĥ= T nuc + V nuc nuc + V nuc el + T el + V el el (2.1) of kinetic energy terms for the a nuclei with masses M a and atomic numbers z a as well as the n electrons with masses m el, T nuc = a a=1 P 2 a 2M a, Tel = n n=1 p 2 n 2m el, (2.2) further, terms accounting for the repulsive electrostatic pair interaction of the nuclei and of the electrons, respectively, a a V nuc nuc = 1 2 a=1 b a z a z b e 2 4πε 0 R a R b, V el el = 1 2 n n n=1 m n e 2 4πε 0 r n r m, (2.3) where e is the elementary charge and ε 0 the vacuum permittivity, and finally the attractive electrostatic interaction between the nuclei and electrons, a n V nuc el = a=1 n=1 z a e 2 4πε 0 R a r n. (2.4) With the Hamiltonian operator Ĥ of the system, the dynamics of nuclei and electrons are governed by the chrödinger equation Ĥ Ψ(R, r, t) =iħ Ψ(R, r, t), (2.5) t where the full wave function of the system (which is sought after) is denoted as Ψ and the sets of nuclear and electronic coordinates are combined in the multi-indices R (R 1,...,R a ) (R 1,...,R 3a ), r (r 1,...,r n ) (r 1,...,r 3n ). (2.6)

44 2.1. eparation of Nuclear and Electronic Motions 43 Adiabatic Electronic Basis ince nuclei have more than three orders of magnitude larger mass than electrons, one expects only little variations in nuclear coordinates while the electronic motion takes place. If the nuclei were entirely fixed, the dynamics of electrons would be governed by an electronic chrödinger equation Ĥ el (R) ψ α (r, R) =E α el (R) ψα (r, R), (2.7) where the electronic Hamiltonian operator (including the interaction with the nuclei) Ĥ el (R) = Ĥ T nuc V nuc nuc = V nuc el + T el + V el el, (2.8) its eigenfunctions ψ α,andtheassociatedenergyeigenvalueseel α depend on the nuclear configuration R only parametrically. Theseeigenfunctionsψ α are known as adiabatic electronic wave functions and their complete orthonormal set as adiabatic electronic basis. The naming stems from the conceptual correspondence between moving electrons at fixed nuclear positions on the one hand and electrons responding instantaneously, or adiabatically, tothemuchslowermovements ofthenucleiontheother. To study the slow nuclear motion in this adiabatic picture, the full wave function Ψ of the system (2.5) can be expanded in the adiabatic electronic basis, Ψ(R, r, t) = α χ α (R, t) ψ α (r, R), (2.9) where the expansion coefficients χ α can be identified with the nuclear wave functions of the system in the respective adiabatic electronic states ψ α. Inserting the expansion (2.9) in the full chrödinger equation (2.5), applying the product rule of differential calculus for the nuclear momentum operators P a = iħ a,andmakinguseof the electronic chrödinger equation (2.7) leads to a nuclear equation of motion Ĥnuc(R) α χ α (R, t) + Θ αβ (R) χ β (R, t) =iħ β α t χα (R, t) (2.10) for each nuclear wave function χ α.clearly,suchawavefunctioncarriesafunctional dependence on the nuclear coordinates R. The same applies to the two operators introduced in (2.10): the nuclear Hamiltonian operator, Ĥnuc(R) α = T nuc + V nuc nuc + Θ αα (R) +Eel α (R), (2.11) U α (R)

45 44 Chapter 2. Foundations of Computational Chemistry acting on χ α and depending in particular on the electronic energy eigenvalue E α el,and the so-called nonadiabaticity operator, Θ αβ (R) = ψ α (R) T nuc ψ β (R) + a a=1 1 M a ψ α (R) P a ψ β (R) P a, (2.12) which acts on the other nuclear wave functions χ β and thereby couples the nuclear equations of motion. Here integrals over electronic coordinates are abbreviated by ψ α (R) ψ β (R) dr ψ α (r, R) ψ β (r, R). (2.13) The dynamics of the nuclei, governed by the nuclear equation of motion (2.10), can be understood as follows: while the system is in a specific adiabatic electronic state ψ α, the nuclear motion is determined by the related nuclear Hamiltonian operator (2.11), hence evolves in the effective potential U α.thefunctionofthispotentialdefinesa hypersurface in the space of nuclear coordinates, which is called adiabatic potential energy surface (pes) [34]. However, while the nuclear motion proceeds, it may occur that the system transitions into another adiabatic electronic state ψ β with associated pes U β. uch transitions, for example from the electronic ground state to the first excited state, are mediated by the nonadiabaticity operator (2.12). Technically, solutions of the full chrödinger equation (2.5) follow from the expansion (2.9), once solutions of the nuclear equation of motion (2.10), that is, wave functions χ αη,withη being the nuclear quantum number, are obtained. The solving of the nuclear equation of motion requires, in turn, the prior determination of the electronic wave functions ψ α and energy eigenvalues Eel α. These solutions of the electronic chrödinger equation (2.7) are needed for each parameter R, i.e., nuclear configuration, covered by the motion of the nuclei. olving the electronic chrödinger equation for specific nuclear positions is the remit of electronic structure theory (ection 2.2). Adiabatic Approximation For many chemical systems, the ratio of electronic and nuclear energy scales is approximately given by (M a /m el ) 1/2 as follows from simple arguments [34, 35]. In cases where this estimate is reliable, the large mass difference between electrons and nuclei leads to an electronic energy spacing at least two orders of magnitude larger than the energy scale of nuclear motion. These systems remain during nuclear motion, in good approximation, in the same eigenstate of the electronic Hamiltonian; in other words, electronic transitions due to the nonadiabaticity operator are negligible. The neglect of the nonadiabaticity operator Θ αβ in the nuclear equa-

46 2.1. eparation of Nuclear and Electronic Motions 45 tion of motion (2.10) is referred to as the adiabatic approximation. If in addition the diagonal terms Θ αα in the nuclear Hamiltonian operator (2.11) are disregarded, the approximation is known as the Born-Oppenheimer approximation [38]. 1 In both cases, the nuclear equation of motion becomes decoupled, yielding for each adiabatic electronic state ψ α an individual nuclear chrödinger equation: [ T nuc + U α (R)] χ α (R, t) =iħ t χα (R, t). (2.14) Note that these approximations are not applicable if electronic energy levels become so close that nuclear vibrational energy suffices for initiating electronic transitions, which is exactly the case in charge transfer reactions in molecular systems (Chapter 3). Classical Nuclei Approximation For nuclei at room temperature, the range of quantum phase coherence, provided by the de Broglie wavelength (2πħ/M a k B T) 1/2, is at least an order of magnitude shorter than typical interatomic distances [35]. Hence, nuclei can often be treated as classical particles, which is accomplished by considering anuclearwavefunctionχ α as Hartree product of incoherent wave functions for individual nuclei, localized at the classical particle positions R a (identified with the mean of the coordinate operator). Then, owing to the Ehrenfest theorem [39], the chrödinger equation (2.14) transforms into a Newtonian equation of motion for the nuclei: M a 2 R a t 2 = a U α (R). (2.15) Numerical integration requires knowledge of the pes U α of the selected electronic state ψ α as a function of R. Inprinciple,onecanconstructaninterpolationbysolving the electronic subsystem for a grid of sampling points in the space of nuclear coordinates. However, since the number of sampling points scales exponentially with the dimension of this space, such an approach is prohibitive in practice. A viable option is to approximate the pes as truncated expansion in terms of many-body interactions [37], U α (R) a a=1 U α 1 (R a )+ 1 2 a a,b a U2 α (R a, R b )+ 1 6 U3 α (R a, R b, R c )+..., (2.16) denoted as force field (ection 4.2.1). The determination of expressions for the included many-body terms allows for classical molecular dynamics simulations (ection 4.2.2). 1 Note that in literature, the Born-Oppenheimer approximation is occasionally not clearly distinguished from the adiabatic approximation, or even the introduction of the adiabatic electronic basis. a,b,c

47 46 Chapter 2. Foundations of Computational Chemistry 2.2. Electronic tructure Theory Electronic structure theory [35, 40 43] aims at finding solutions to the electronic chrödinger equation of a many-body problem, such as a molecular system, with fixed nuclear positions. This electronic chrödinger equation, provided in (2.7), implies that for any solution ψ, the associated electronic energy is determined by E el [ψ] = ψ Ĥ el ψ dr ψ (r) Ĥ el ψ(r), (2.17) that is, the expectation value of the electronic Hamiltonian operator (2.8), which reads Ĥ el = n n=1 ˆv nuc el (r n ) V nuc el + n ˆt el (r n ) n=1 T el + 1 n n 2 n=1 m n ˆv el el (r n, r m ) V el el. (2.18) Here, the following terms, corresponding to the definitions (2.2) to (2.4), are introduced: ˆv nuc el, denoting the one-electron energies in the potential of the fixed nuclei, ˆt el,theone-electronkineticenergies,andˆv el el,thetwo-electroninteractionenergies. 2 The following discussion is limited to the task of retrieving the electronic ground state, that is, the wave function with lowest energy. For this case, a solving strategy is the minimization of the energy functional (2.17) with respect to a trial set of wave functions. Wave functions coming into consideration are normalized as well as antisymmetric under the exchange of two electrons. The latter originates from the Pauli principle applying to fermions. 3 In principle, the most general antisymmetric wave function, defining a complete trial set, can be represented as a linear combination of an infinite set of later determinants. In practice, however, tractable approximations are required. In what follows, the traditional Hartree-Fock theory is outlined (ection 2.2.1). Here, the simplest antisymmetric wave function is constructed, a single later determinant. This yields a model of independent electrons, each moving inamean field of all other electrons. While this approach includes the exchange interaction between electrons, their dynamical correlation is neglected. In the subsequently discussed Kohn-ham density functional theory (ection 2.2.2), these exchange and correlation contributions are separated in an energy functional of the electron density. Thisso-calledexchangecorrelation functional allows the dynamical correlation of electrons to be treated in an approximate way. Finally, the basis set approximation is introduced (ection 2.2.3). 2 Note that in (2.17) and (2.18) the parametric dependences on the nuclear coordinates are omitted. 3 The Pauli principle represents an independent postulate within non-relativistic quantum mechanics.

48 2.2. Electronic tructure Theory Hartree-Fock Theory Amulti-electronwavefunctionofn electrons can be constructed as a product, a socalledhartree product, of n independent wave functions for the individual electrons. These one-electron wave functions ϕ n are denoted as molecular orbitals, each describing a single electron with space and spin coordinate r n. 4 To ensure the antisymmetry of the n-electron wave function, i.e., that the exchange of two electron coordinates causes only a sign change, the Hartree product can be antisymmetrized by a later determinant. This simplest antisymmetric wave function is used in Hartree-Fock theory: ψ hf (r)= 1 n! ϕ 1 (r 1 ) ϕ 2 (r 1 ) ϕ n (r 1 ) ϕ 1 (r 2 ) ϕ 2 (r 2 ) ϕ n (r 2 ). (2.19) ϕ 1 (r n ) ϕ 2 (r n ) ϕ n (r n ) Inserting this n-electron wave function in the energy functional (2.17), with the electronic Hamiltonian operator provided by (2.18), that is, working out the expression E el [ψ hf ]= ψ hf V nuc el ψ hf + ψ hf T el ψ hf + ψ hf V el el ψ hf, (2.20) allows one tocollect similarintegrals over the n individual electronic coordinates and to abandon their numbering. The first two terms of (2.20), comprising the one-electron energy contributions ˆv nuc el and ˆt el,simplybecomesumsofthesecontributions: 5 n ψ hf V nuc el ψ hf = 1 ϕ n (r 1 ) ˆv nuc el (r 1 ) ϕ n (r 1 ), n=1 (2.21) ψ hf T n el ψ hf = 1 ϕ n (r 1 ) ˆt el (r 1 ) ϕ n (r 1 ). n=1 (2.22) The third term of (2.20), that is, the electron-electron interaction composed of the two-electron contributions ˆv el el,yieldstwoparts: n ψ hf V el el ψ hf = 1 2 dr 1 dr 2 ϕ n (r 1 )ϕ m (r 2 ) ˆv el el (r 1, r 2 ) ϕ n (r 1 )ϕ m (r 2 ) (2.23) n,m 1 2 n n,m dr 1 dr 2 ϕ n (r 1 )ϕ m (r 2 ) ˆv el el (r 1, r 2 ) ϕ m (r 1 )ϕ n (r 2 ),(2.24) 4 Here, r n =(r 1 n, r 2 n, r 3 n, σ n) captures the three spatial and one spin coordinate for electron n collectively. 5 In Hartree-Fock literature the terms (2.21) and (2.22) are typically combined into a core Hamiltonian.

49 48 Chapter 2. Foundations of Computational Chemistry where the first part, denoted as the Hartree energy,represents the Coulomb energy due to all independent electrons and the second part is referred to as the exchange energy. Now, to find the wave function ψ hf with lowest energy, the energy functional (2.20) is minimized with respect to variations of the orbitals ϕ n,whichareconstrainedto orthonormality. Employing Lagrange multipliers ε n,henceforthinterpretedasorbital energies, allows the derivation of an eigenvalue equation for the orbitals, known as the Hartree-Fock Equation Ĥ hf ϕ n (r 1 )=ε n ϕ n (r 1 ), (2.25) Ĥ hf = ˆt el (r 1 )+ˆv nuc el (r 1 )+ˆv h (r 1 )+ˆv x (r 1 ), (2.26) where Ĥ hf is called the Fock operator. Itsfourcontributionsarerelatedtothefourenergy terms (2.21) to (2.24) and the newly specified ones act on the orbitals as follows: n ˆv h (r 1 ) ϕ n (r 1 ) = [ dr 2 ϕ m (r 2 ) ˆv el el (r 1, r 2 ) ϕ m (r 2 )] ϕ n (r 1 ), (2.27) m=1 ˆv x (r 1 ) ϕ n (r 1 ) = [ n m=1 dr 2 ϕ m (r 2 ) ˆv el el (r 1, r 2 ) ϕ n (r 2 )] ϕ m (r 1 ). (2.28) The Fock operator is an effective one-electron operator that governs the dynamics of the independent electrons. Each electron experiences the fixed nuclear potential ˆv nuc el, as well as the Hartree potential ˆv h and the exchange potential ˆv x.thehartreepotential mimics theexact electron-electroninteractionbyamean field of all other independent, or noninteracting electrons, whose charge can be seen as spread over the system. The exchange potential describes a certain interaction the exchange interaction: it prevents the near vicinity of parallel-spin electrons, and thus reflects the Pauli principle enforced bytheantisymmetrization of the wave function. However, the Fock operator neglects any dynamical correlation of electrons due to their electrostatic interaction. In other words, opposite-spin electrons can come arbitrarily close to each other. Technically, the Fock operator depends on its own eigenfunctions, i.e., the orbitals which are actually determined by the Hartree-Fock equation. Hence, the Hartree- Fock eigenvalue equation needs to be solved self-consistently (ection 2.2.3). Once self-consistent solutions are obtained, the later determinant ψ hf constructed from the n orbitals with lowest orbital energies provides the Hartree-Fock ground state.

50 2.2. Electronic tructure Theory Kohn-ham Density Functional Theory The electronic Hamiltonian operator (2.18), describing a many-body system such as a molecular system, is completely defined by the specification of both the fixed nuclear coordinates, i.e., the nuclear potential V nuc el,andthenumberofelectronsn. The key insight giving rise to density functional theory is that the Hamiltonian operator is alternatively determined by the knowledge of the ground state electron density ρ,that is, the probability of finding any of the n electrons at a particular point in space: dr 1 ρ(r 1 )=n. (2.29) This is proven by the first Hohenberg-Kohn theorem [44], providing the equivalent statement that the electron density ρ (which already specifies n) uniquely determines the nuclear potential V nuc el, thus the Hamiltonian operator, and therefore all information of the system, in particular its ground state energy through a necessarily existing functional (2.30). A second theorem states that this functional obeys the variational principle, i.e., is minimized by ρ yielding the ground state energy: E el [ρ] = dr 1 ρ(r 1 ) ˆv nuc el (r 1 )+T el [ρ] +V el el [ρ]. (2.30) As opposed to the nuclear potential, one lacks explicit expressions of the kinetic energy T el and the electron-electron interaction V el el in terms of the electron density. Within the Kohn-ham approach [45], this central problem of density functional theory is addressed by introducing an auxiliary system of independent, noninteracting electrons with molecular orbitals ϕ n,representingthesame electron density: ρ(r 1 )= n n=1 ϕ n (r 1 ) ϕ n (r 1 ). (2.31) In this way, portions of the unknown energy functionals, corresponding to the noninteracting auxiliary system, can be expressed in terms of the orbitals using the results of Hartree-Fock theory for the kinetic energy (2.22) and the Hartree energy (2.23), T el [ρ] = n n=1 dr 1 ϕ n (r 1 ) t el (r 1 ) ϕ n (r 1 )+ T el [ρ], (2.32) V el el [ρ] = 1 2 dr 1 dr 2 ρ(r 1 ) ρ(r 2 ) ˆv el el (r 1, r 2 )+Ṽ el el [ρ], (2.33)

51 50 Chapter 2. Foundations of Computational Chemistry while all remaining many-body contributions, due to the exchange interaction and dynamical correlation, are incorporated in the exchange-correlation energy functional of the density,whichis apurelyformal expression at this point: E xc [ρ] = T el [ρ] +Ṽ el el [ρ]. (2.34) To find an expression of the ground state density (2.31), the energy functional (2.30) is minimized with respect to variations of the auxiliary orbitals ϕ n,constrainedtoorthonormality. The method of Lagrange multipliers ε n leads to an eigenvalue equation for the orbitals, which requires self-consistent solving (ection 2.2.3) and is called the Kohn-ham Equation Ĥ ks ϕ n (r 1 )=ε n ϕ n (r 1 ), (2.35) Ĥ ks = ˆt el (r 1 )+ˆv nuc el (r 1 )+ˆv h (r 1 )+ˆv xc (r 1 ). (2.36) Here, the effective one-electron operator Ĥ ks, denoted as the Kohn-ham operator, describes the electron-electron interaction by the Hartree potential ˆv h,representingthe electronic mean-field, plus the exchange-correlation potential ˆv xc,whichisthefunctional derivative of the exchange-correlation energy functional (2.34): ˆv h (r 1 )= dr 2 ρ(r 2 ) ˆv el el (r 1, r 2 ), ˆv xc (r 1 )= δe xc[ρ] δρ(r 1 ). (2.37) Exchange-Correlation Energy Functional In principle, self-consistent solutions of the Kohn-ham equation give access to the exact ground state of the interacting many-body system, provided an exact expression for the exchange-correlation energy functional E xc is used. Unfortunately, such an expression is unknown but anumberofapproximationsareavailable. Thesimplestoftheseisthelocal density approximation (lda), where the electron density is approximated as locally homogeneous in order to use the model of the homogeneous electron gas at the same density, for which the exchange functional Ex lda is known analytically [46] and the correlation functional Ec lda can be parametrized with high accuracy [47]. To address the issue of the still neglected nonlocal density dependences of the exchange-correlation energy, improved approximations can be performed. In a first step, semilocal dependences can

52 2.2. Electronic tructure Theory 51 be taken into account by establishing expressions for the density gradient and higher order derivatives, referred to as the generalized gradient approximation (gga); examples are the exchange functional Ex gga of Becke [48] and the correlation functional Ec gga of Lee, Yang, and Parr [49]. The next step of accounting for truly nonlocal density dependences poses the greatest challenge; however, at least for the exchange energy, the nonlocal expression given by Hartree-Fock theory Ex hf (2.24) can be incorporated in the description, resulting in hybrid Hartree-Fock Kohn-ham exchange-correlation energy functionals. Particularly successful is the Becke three-parameter hybrid functional with Lee, Yang, and Parr correlation (b3lyp) [50], which is composed as E b3lyp xc = E lda x + E lda c + α 1 (E hf x Ex lda )+α 2 (E gga x E lda x )+α 3 (Ec gga Ec lda ).(2.38) The parameters α 1 = 0.2, α 2 = 0.72 and α 3 = 0.81 are determined by fitting to experimental data for atomization energies, ionization potentials, and proton affinities of up to 56 different molecules, and total energies of 10 atoms [50]. In this work, Kohn-ham density functional theory with the b3lyp hybrid functional is used, inter alia, to obtain many-body terms, as in (2.16), for the parametrization of force fields (ection 4.2.1), to determine molecular point charge distributions and polarizabilities for the calculation of site energies (ection 4.2.5), and to assess the nuclear rearrangement of molecules upon electron attachment for the quantification of reorganization energies (ection 4.2.6) Basis et Approximation The discussed electronic structure theories reduce the electronic chrödinger equation of the n-electron system to the Hartree-Fock (2.25) and Kohn-ham equation (2.35). Both of these represent a set of one-electron eigenvalue equations for molecular orbitals ϕ n and associated orbital energies ε n,governedbyaneffectiveone-electron Hamiltonian operator Ĥ,whichdependsonitsowneigenfunctions: Ĥϕ n (r 1 )=ε n ϕ n (r 1 ). (2.39) This type of equation is routinely solved after conversion into an algebraic equation by means of a finite set of predefined basis functions φ τ. These basis functions are often chosen as atomic orbitals of the isolated atoms of the system, or are based on them, and are therefore conventionally denoted as atomic orbitals. Byexpandingthe

53 52 Chapter 2. Foundations of Computational Chemistry unknown molecular orbitals as a linear combination of the atomic orbitals (mo-lcao), only the expansion coefficients C τn remain to be determined: ϕ n (r 1 )= C τn φ τ (r 1 ). (2.40) τ Now, inserting this expansion in the eigenvalue equation (2.39), multiplying from the left with φ σ and integrating over r 1 leads to the desired algebraic equation, also known as the Roothaan-Hall matrix equation [51, 52], HC = CE, (2.41) for the expansion coefficients C =(C τn ) and molecular orbital energies E = diag(ε n ). H =(H στ ) is the Hamiltonian matrix and =( στ ) the atomic orbital overlap matrix, H στ = φ σ Ĥ φ τ dr 1 φ σ (r 1 ) Ĥφ τ (r 1 ), (2.42) στ = φ σ φ τ dr 1 φ σ (r 1 ) φ τ (r 1 ), (2.43) where the latter arises, since atomic orbitals are generally nonorthogonal, which implies that nonvanishing terms of differential overlap,that is,for σ τ,exist.the application of a similarity transformation T, whichorthogonalizesthebasissetofatomic orbitals T T = I, leading to a new coefficient matrix C = T 1 C and a transformed one-electron Hamiltonian matrix H = T HT,resultsinamatrixequation[40], H C = C E, (2.44) which yields the new coefficients and the orbital energies by diagonalization of H. Of course, this matrix still depends itself on the coefficients sought, hence the diagonalization process requires an initial guess of the coefficients. To obtain a self-consistent solution, i.e., coefficients matching the initial guess, the diagonalization process needs to be iterated until the coefficient matrix is converged to a desired level of precision. These iterative methods are implemented in standard software packages for computational chemistry, such as the gaussian package [53], which is used in this work. These methods of electronic structure theory are often referred to as ab initio or first principles methods, although density functional theory may include certain empirical parameters, depending on the employed exchange-correlation energy functional.

54 2.2. Electronic tructure Theory 53 Basis ets In practice, computational chemistry packages provide a comprehensive choice of basis sets. Particularly popular are the Pople style basis sets [54], whose basis functions are, as is customary, based on the atomic orbitals of the isolated atoms of the system. The correct functional form for describing atomic orbitals are later type orbitals, whose radial part decays exponentially with the distance from the atomic nucleus. However, for computational efficiency, it is beneficial to approximate a later type orbital by multiple Gaussian type orbitals, which decay exponentially with the squared distance. According to this approximation, the minimal Pople style basis sets, designated as sto-ng, describe the later type orbitals by n Gaussian type orbitals. These minimal basis sets are, however, often unsuitable for molecular systems, since chemical bonding in molecules breaks the spherical symmetry of the atomic orbitals. To address this issue, the valence orbitals, which facilitate the bonding effects, can be split into two or more functions with varying decay factors. The corresponding split-valence Pople style basis sets are designated as n-mlg, where n indicates the number of Gaussian type orbitals used to represent the core orbitals, while each of the subsequent variables m, l, etc. referstooneofthefunctionsintowhichthevalence orbitals are split. The basis sets can be further improved by including polarization functions, indicated by an asterisk, while two asterisks signify that they are added to light atoms as well. In this work, the 6-311g basis set is typically employed. emiempirical Methods Based on Hartree-Fock theory, there exists a variety of semiempirical methods [55 58], which incorporate further approximations and empirical information in order to reduce the computational expense of ab initio methods. On the one hand, such semiempirical methods often involve what is called the neglect of differential overlap (ndo), where the calculation of the atomic orbital overlap matrix (2.43) is initially avoided, that is, στ = δ στ. This implies that the Roothaan- Hall matrix equation (2.41) is directly in the standard eigenvalue form (2.44). On the other hand, the matrix elements (2.42) of the Fock operator (2.26) are replaced by significantly simplified expressions including empirical relationships. These simplified expressions may also include some special, again empirically parametrized, atomic orbital overlap terms στ to compensate for the major errors resulting from the ndo. An example is Zerner s intermediate neglect of differential overlap (zindo) [59, 60], which is referred to at a later point; in this work, it is employed to calculate coupling strengths between electronic states, specifically the transfer integrals (ection 4.2.4).

55 54 Chapter 2. Foundations of Computational Chemistry

56 Part II. Methodology

57

58 Chapter 3. Charge Transfer in Molecular ystems Charge transfer [34, 61 64] refers to a microscopic process of electron movement, which is elementary for many macroscopic charge transport phenomena. The microscopic systems hosting charge transfer processes are molecular systems,such as molecular fragments, entire molecules, or their aggregates. A more precise definition of such a microscopic system is the charge transfer complex. It denotes a molecular system comprising a negative or positive excess charge, i.e., an electron or hole, which is localized to a certain degree, such that (in the simplest case) two spatial localization centers can be identified: the donor and the acceptor. The process of charge transfer is then understood as the spatial displacement of the excess charge from the donor to the acceptor part of the charge transfer complex, and is written, in the case of electron transfer, as d a da. ince such a charge redistribution alters the electrostatic field within the complex, it is associated with a rearrangement of the nuclear coordinates of the molecular system. This leads to a central property of charge transfer: the two nuclear configurations corresponding to the charge localization at the donor and acceptor can be thought of as minima in the potential energy surfaces for the nuclear motion. The energetic barrier separating these two minima acts as reaction barrier for the transfer process and critically determines the characteristic reaction time, in other words the charge transfer rate. The formulation of analytical expressions for the charge transfer rate is the main objective of charge transfer theories. To this end, the key features characterizing the charge transfer complex and its environment need to be translated into tangible parameters. The discussion starts by introducing the regimes of charge transfer (ection 3.1). This overview illustrates that a charge transfer reaction should be treated on a different theoretical basis, depending on the degree of localization of the excess charge. Based on these considerations, it follows the explicit formulation of charge transfer rates in the regime of strongly localized charges (ection 3.2).

59 58 Chapter 3. Charge Transfer in Molecular ystems 3.1. Regimes of Charge Transfer The microscopic processes occurring in a charge transfer complex are governed by the chrödinger equation of its constituent electrons and atomic nuclei (ection 2.1). uch a full chrödinger equation (2.5) can be expanded by an adiabatic electronic basis of delocalized states, defined by an electronic motion with fixed nuclei, i.e., a parametric dependence on the nuclear coordinates. This leads to separate but coupled equations of motion for electrons (2.7) and nuclei (2.10), which may be decoupled by the Born-Oppenheimer (or adiabatic) approximation. Now, considering a charge transfer reaction, the validity of the Born-Oppenheimer approximation is not expectable in general since the charge displacement and nuclear rearrangement occur concertedly. Moreover, even the use of the adiabatic electronic basis can be unsuitable for representing electronic wave functions describing the shift of charge localization. Hereafter, the adiabatic electronic basis is applied to describe the transfer of rather delocalized charges, which is the regime of adiabatic charge transfer (ection3.1.1). The transfer of rather localized charges, in contrast, falls into the regime of nonadiabatic charge transfer (ection 3.1.2). In this case, it proves beneficial to choose an electronic basis of localized states, which is denoted as diabatic electronic basis Adiabatic Charge Transfer To treat the present case of two localization centers (the donor and acceptor), it suffices to introduce an adiabatic electronic basis in two-state approximation, consisting of two delocalized states ψ ± (Q). These states correspond to adiabatic wave functions ψ ± (r, R)= r ψ ± (R) defined in (2.7), which are parametrically dependent on the nuclear coordinates R. Notethatthis3a-dimensional set of nuclear coordinates is here expressed by a one-dimensional linear combination Q,the reaction coordinate. While detailed later, it can for now be imagined to pass through the two nuclear configurations, where the excess charge is localized at the donor and the acceptor. Hamiltonian in Adiabatic Representation In the thus-defined adiabatic electronic basis, the full Hamiltonian of the charge transfer complex expands as follows: 1 Ĥ= Ĥnuc(Q) ψ α α (Q) ψ α (Q) + Θ αβ (Q) ψ α (Q) ψ β (Q). (3.1) α=± α=± β 1 This results from the nuclear equation of motion (2.10), as seen by taking γ (3.1) ψ γ (Q) χ γ (Q), using the full chrödinger equation (2.5) and the expanded full wave function (2.9).

60 3.1. Regimes of Charge Transfer 59 U (Q) U + (Q) 2J E Q Q Figure 3.1. Adiabatic potential energy surfaces of a charge transfer complex as functions of the reaction coordinate Q. ThepotentialsU and U + govern the nuclear motion in the adiabatic electronic basis states ψ and ψ +. At the nuclear configuration Q, the functions approach to a minimum separation of 2J and the potential U exhibits an energetic barrier E,separatingtwominima. This implicates that the electronic coupling between the two adiabatic basis states ψ and ψ + is caused by the nonadiabaticity operator Θ ±,whichisgivenin(2.12). If the system is in one of these basis states, its nuclear motion is governed by the respective nuclear Hamiltonian operator, which, according to (2.11), reads Ĥ nuc(q) ± = T nuc + V nuc nuc + Θ ±± (Q) +E el ± (Q). (3.2) U ± (Q) A prototypical example of the corresponding adiabatic potential energy surfaces (pes) U and U + is depicted in Figure 3.1. Here, the pes of the electronic ground state U exhibits two minima, separated by a barrier E at the nuclear configuration Q. Adiabatic Charge Transfer Reaction Now, if the excess charge is rather delocalized, the charge transfer complex is in the almost pure electronic ground state ψ, while electronic coupling to the next state is only a small term. 2 Then, the nuclear dynamics is governed by the pes U and is conceivable as a vibration around the, say, left minimum until, after a characteristic time, the barrier to the right minimum can be surmounted. During this nuclear rearrangement along the reaction coordinate Q, the electronic state ψ whichisaparametricfunctionofq is continuously deformed, such that the excess charge moves from the donor to the acceptor part of the charge transfer complex. This continuous modification of the electronic probability density of a single electronic state is denoted adiabatic charge transfer. Theassociatedcharge transfer rate is proportional to the nuclear transition rate, and hence takes, according to transition state theory [65, 66], the simple form of an Arrhenius relation: ω e E /k B T. (3.3) 2 In Figure 3.1 this corresponds to an energy splitting larger than nuclear energy quanta: J ħω vib.

61 60 Chapter 3. Charge Transfer in Molecular ystems Nonadiabatic Charge Transfer If the excess charge, in contrast to the preceding case, is rather localized, the system is no longer in an almost pure delocalized state of the adiabatic electronic basis. Instead, one expects substantial mixing of these states, in other words strong electronic coupling, mediated by the nonadiabaticity operator. As the nonadiabaticity operator (2.12) is a differential operator with respect to Q, the derivation of an eligible charge transfer rate is no simple task. It is therefore advisable to introduce an electronic basis of localized states, which is referred to as diabatic electronic basis. Diabatic Electronic Basis taying in the two-state approximation, diabatic electronic states ψ i,f are inferred from the adiabatic ones ψ ± by a basis transformation, ( ψi (Q 0 ) ψ f (Q 0 ) ) = T(Q)( ψ (Q) ψ + (Q) ), (3.4) where the unitary matrix T is chosen such as to eliminate their parametric dependence on the nuclear coordinates Q, in other words to inhibit any electronic coupling due to the Q-derivative nonadiabaticity operator [67 69], which, in turn, enables electronic localization. 3 Typically, the diabatic states depend on a fixed nuclear reference configuration Q 0 (henceforth omitted), irrespective of the actual nuclear positions. Thus, they are no longer eigenstates of the electronic Hamiltonian (2.8), except at Q 0 [34], Ĥ el (Q 0 ) ψ i,f =E i,f el ψ i,f, (3.5) with eigenvalues E i,f el, while the parametric dependence on Q is separated in a term V, Ĥ el (Q 0 )=Ĥ el (Q) V(Q, Q 0 ). (3.6) Hamiltonian in Diabatic Representation In the thus-defined diabatic electronic basis, the full Hamiltonian of the charge transfer complex expands as follows: 4 Ĥ=Ĥ i nuc(q) ψ i ψ i +Ĥ f nuc(q) ψ f ψ f +J if (Q) ( ψ i ψ f + ψ f ψ i ). (3.7) 3 One can show that the absence of the nonadiabatic coupling terms is related to the diagonalization of operators associated with charge localization, such as the dipole moment operator [70 75]. 4 The expression (3.7) results from expanding Ĥ=Î Ĥ Î using the completeness relation of the diabatic electronic basis and the definition of the electronic Hamiltonian operator (2.8).

62 3.1. Regimes of Charge Transfer 61 U (Q) E U + (Q) Q 2J U i (Q) U f (Q) Q Figure 3.2. Diabatic potential energy surfaces of a charge transfer complex as functions of a reaction coordinate Q. The potentials U i and U f (solid black and red lines) govern the nuclear motion in the diabatic electronic basis states ψ i and ψ f. The functions are crossing at the nuclear configuration Q, defining the energetic barrier E.In contrast, the adiabatic potential energy surfaces (dashed gray lines) avoid crossing at Q with an energetic separation of 2J. While any terms with the nonadiabaticity operator vanish by definition, one now encounters electronic coupling between the diabatic basis states ψ i and ψ f due to the Transfer Integral J if (Q) = ψ i Ĥ el (Q) ψ f = ψ i V(Q, Q 0 ) ψ f. (3.8) If the system is in one of the diabatic basis states, the nuclear dynamics is described by the respective nuclear Hamiltonian operator, which is given by Ĥ i,f nuc(q) = T nuc + V nuc nuc + J ii,ff (Q) +E i,f el U i,f (Q). (3.9) Prototypical diabatic potential energy surfaces U i and U f are depicted in Figure 3.2. Apparently, these two functions are crossing at the nuclear configuration Q and give rise to the shape of their adiabatic counterparts, which avoid crossing at this point. Nonadiabatic Charge Transfer Reaction Now, if the excess charge is rather localized, the charge transfer complex is in an almost pure diabatic basis state, say, ψ i, while the electronic coupling to the other one, i.e., the transfer integral, is only a small term. 5 Then, nuclear vibrational motion is governed by the pes U i and leads, after a certain time, to the nuclear transition state Q.There,anelectronic transition to the state ψ f can take place, 6 such that nuclear motion proceeds in the pes U f.notethat 5 In Figure 3.2 this corresponds to an energy splitting smaller than nuclear energy quanta: J ħω vib. One can show that this splitting equals the transfer integral at the transition state: J = J if (Q ) [34]. 6 An electronic transition is restricted to this diabatic crossing point due to the simultaneous fulfillment of the Franck-Condon principle [76], i.e., the condition of constant nuclear coordinates (vertical transition in Figure 3.2), and the energy conservation (horizontal transition).

63 62 Chapter 3. Charge Transfer in Molecular ystems nuclear motion along the reaction coordinate Q causesno alteration of the basis states, which follows from their independence of Q. Instead,itisthetransitionbetweenthe basis states which describes the movement of the excess charge from the donor to the acceptor. This instantaneous modification of the electronic probability density due to the transition between electronic states is denoted nonadiabatic charge transfer. The associated charge transfer rate not only reflects the probability of reaching the nuclear transition state Q, which is the familiar Arrhenius relation (3.3), but also the probability for the occurrence of the electronic transition mediated by the transfer integral: ω J if 2 e E /k B T. (3.10) When deriving rate expressions of this type, the transfer integral as a small quantity can be conveniently treated as a perturbation to the initial state (ection 3.2). Reaction Coordinate What follows is a supplement to the one-dimensional reaction coordinate Q,whichcombinesthe3a-dimensional set of nuclear coordinates R. To describe nuclear vibrational motion in a pes U(R) around some minimum R 0,one may approximate the pes in the vicinity of R 0 by a second order Taylor expansion [34], U(R)=U h ab (R a R 0 a)(r b Rb 0 )=U a,b ξ (ω vib ξ (Q ξ Q 0 ξ ))2 = U(Q),(3.11) which allows the transformation, by diagonalization of the Hessian matrix (h ab ),to 3a 6 (mass weighted) normal mode coordinates Q ξ with eigenfrequencies ω vib ξ.note that 6 eigenvectors represent overall translational and rotational motions. Then, along with the conjugate momentum operators P ξ, the nuclear Hamiltonian operator reads Ĥ nuc (Q)= 1 2 ξ P 2 ξ + U(Q). (3.12) Achargetransferreactionisusuallytriggeredbythenuclearmotionofasmallsubset of normal mode coordinates, the active normal mode coordinates. An appropriate linear combination of these active coordinates often yields a suitable one-dimensional reaction coordinate Q with average eigenfrequency ω vib and conjugate momentum operator P, such that the nuclear Hamiltonian (3.12) takes the simple form Ĥ nuc (Q)= 1 2 P2 + U(Q), U(Q)=U (ωvib (Q Q 0 )) 2. (3.13)

64 3.2. Charge Transfer Rates Charge Transfer Rates A charge transfer reaction can take place in numerous variants, which is reflected in the existence of a wide variety of rate expressions, tailored to the characteristics of the charge transfer complex. These distinctive features are worth mentioning. A first distinction is the degree of localization of the excess charge: according to the discussion on the regimes of charge transfer (ection 3.1), this specifies if the charge transfer reaction isanadiabatic or nonadiabatic process and suggests the appropriate, i.e., adiabatic or diabatic, representation of the Hamiltonian operator. A second factor is the actual composition of the molecular system: for unimolecular systems, where the donor and acceptor are parts of the same molecule, the nuclear vibrational dynamics can be described by a common set of normal mode coordinates, which usually leads to a single reaction coordinate. For bimolecular systems, in contrast, it is generally required to introduce two separate sets of normal mode coordinates with independent reaction coordinates. A third aspect is the surrounding of the charge transfer complex: if the complex is embedded in a static environment,it sufficestoaccountfor changes inthe electrostatic interactions with the environment, related to the charge redistribution in the complex. However, if the complex is formed in a polar solvent, thereorientation of surrounding solvent molecules causes significant nuclear adjustments, which requires the introduction of an outer-sphere reaction coordinate in addition to the internal ones. A fourth important system property is the temperature: it allows the determination of whether the nuclear vibrations along the internal and outer-sphere reaction coordinates should be described by classical or quantum mechanics. In a hightemperature limit, the vibrations related to all reaction coordinates can be treated classically. At medium temperatures, thehigh-frequencyvibrationsalongtheinternalreaction coordinates should be treated quantummechanically,while the low-frequency dynamics related to the outer-sphere reaction coordinate may still be described classically. In a low-temperature limit, the nuclear vibrations may be entirely neglected. In the studies of this work, charge transfer occurs between individual molecules, i.e., in bimolecular systems, surrounded by solid material, i.e., a static environment,which is at room temperature, identified as high-temperature limit, and exhibits structural disorder, causing strong charge localization and hence nonadiabatic charge transfer. For this particular type of charge transfer, the derivation of an appropriate rate expression is detailed (ection 3.2.1). This is followed by a brief comparison to rates in the medium and low temperature limits of nonadiabatic charge transfer (ection 3.2.2).

65 64 Chapter 3. Charge Transfer in Molecular ystems Bimolecular High-Temperature Nonadiabatic Charge Transfer The theoretical elaboration of the charge transfer reaction starts from the Hamiltonian operator of the charge transfer complex in diabatic representation (3.7). Accordingly, the electronic basis consists of an initial and a final electronic state (i and f), corresponding to charge localization at the donor and acceptor, respectively. In treating thebimolecular case, the donor and acceptor are taken as two independent molecules (i and j), whose nuclear vibrational motions are described by two separate sets of normal mode coordinates, which are collected in two independent (mass weighted) reaction coordinates Q i and Q j with conjugate momentum operators P i and P j.hence, the nuclear Hamiltonian operators for the initial and the final electronic states (3.9) contain two independent contributions of kinetic energy, potential energy and of interaction energy with the environment, where the latter accounts for a static environment of the charge transfer complex by constant electrostatic interactions Wi, c,n j with the two molecules in their charged and neutral states (c and n): Ĥnuc(Q i i, Q j )= 1 2 P2 i P2 j + Ui(Q i i )+U i j(q j )+Wi c + Wj n, E i (Q i,q j ) Ĥnuc(Q f i, Q j )= 1 2 P2 i P2 j + Ui f (Q i )+U f j(q j )+Wi n + Wj c. E f (Q i,q j ) (3.14) Here, E i and E f represent two-dimensional potential energy surfaces, which, in a harmonic oscillator approximation, are expressed by quadratic functions with respect to the individual reaction coordinates of the two molecules, such as in (3.13): U i,f i, j(q i, j )=U c,n i, j (ωvib i, j (Q i, j Qc,n i, j ))2. (3.15) Figure 3.3 provides an illustration of these four quadratic functions and defines the minimum energies Ui, c,n j at the nuclear equilibrium configurations Qi, c,n j given in (3.15). The functions Ui i and Ui j, drawn in black, govern the nuclear vibrations of the two molecules (around the configurations Qi c and Q n j ) in the initial electronic state of the charge transfer complex, where the excess charge is localized at the donor molecule i. In the final electronic state of the complex, where the charge has passed to the acceptor molecule j, thenuclearvibrationsofthetwomolecules(aroundq i n and Q c j )are described by the functions Ui f and Uf j,markedinred.

66 3.2. Charge Transfer Rates 65 (a) Ui(Q i i ) U f (b) i (Q i ) U f j(q j ) U i j(q j ) U n i U c j λ i λ j E int i Q i U c i U n i E int j Q c i Q n i Q n j Q c j Q j U c j U n j Figure 3.3. (a) Diabatic potential energy surfaces for the donor molecule i, (b) for the acceptor molecule j. For each molecule, the potential energy surface of the initial electronic state is drawn in black, while that of the final electronic state is drawn in red. Adapted from Reference [4]. The high-temperature limit refers to the limiting case in which the energy spacing of the harmonic oscillators is small compared to the thermal energy (ħω vib i, j k B T) and the nuclear vibrations can be treated classically. Then, an average treatment of the kinetic energy is possible by the reaction coordinate distribution functions f i, j (Q i, j )= 1 Z i, j e U i i, j (Q i, j )/k B T, Z i, j = dq i, j e U i i, j (Q i, j )/k B T. (3.16) Characteristic fornonadiabatic charge transfer is a weak electronic coupling between the initial and final diabatic electronic states: this allows the transfer integral J if (3.8) to be treated as a perturbation to the initial state, switched on at time zero. Hence, the time-dependent electronic chrödinger equation, describing the localization shift of the electronic probability density from the donor to the acceptor, can be expressed by means of time-dependent perturbation theory. In first-order approximation, Fermi s golden rule [76 78] provides the transition rate from the initial electronic state with energy E i to the final state with energy E f, 7 where summations over quantum level manifolds can be replaced by integrations over the averaged vibrational energies: ω if = 2π ħ Jif 2 dq i dq j f i (Q i ) f j (Q j ) δ (E i (Q i, Q j ) E f (Q i, Q j )). (3.17) 7 To apply Fermi s golden rule, the transfer integral J if is considered constant with respect to the reaction coordinate, which is referred to as the Condon approximation [79].

67 66 Chapter 3. Charge Transfer in Molecular ystems To carry out the integration (3.17), one inserts the expressions (3.14), (3.15), and (3.16) and demands that the argument of the delta function vanishes [80 82]. This leads to achargetransferratethatisreferredtoasthe Marcus Rate ω ij = 2π ħ Jij 2 4πλij k B T exp ( E ij λ ij ) 2 4λ ij k B T. (3.18) Notice that the rate ω ij and the transfer integral J ij are henceforth labeled by the index pair of the molecules at which the excess charge is localized (i and j)intherespective diabatic electronic state (i and f). The same applies to the energy difference between the nuclear equilibrium configurations of these electronic states, that is, the driving force for the charge transfer reaction, which can be expressed as the ite Energy Difference E ij = (E int i + E elstat int i ) (E j + E elstat j ), (3.19) E i E j where energetic contributions are grouped with respect to molecules i and j, defining the site energies E i and E j.thesiteenergyofaparticularmoleculecomprisesaninternal and an electrostatic contribution, each of which characterizing the charged state (c) of the molecule relative to its neutral state (n), as indicated in Figure 3.3: E int i, j = Uc i, j Un i, j, E elstat i, j = W c i, j Wn i, j. (3.20) In addition, the Marcus rate accounts for the energetic expense associated with the nuclear rearrangement from the equilibrium configurations to the nuclear transition states, located at the diabatic crossing points. This is covered by the Reorganization Energy λ ij = (Ui n Ui n ) c + (U j U c j ). (3.21) λ i λ j

68 3.2. Charge Transfer Rates 67 The reorganization energy λ ij of the charge transfer complex can be considered as a sum of reorganization energies λ i and λ j for the donor and acceptor molecules. These molecular reorganization energies denote the energy of a molecule in the final electronic state, but at the nuclear equilibrium configuration of the initial electronic state, as is indicated in Figure 3.3 by gray dots. In other words, the reorganization energy λ i represents the energy of molecule i in its neutral state with nuclear configuration of its charged state (n ), and λ j is the energy of molecule j in its charged state with nuclear configuration of its neutral state (c ). These energetic values contain the information about the curvature and mutual offset of the potential energy surfaces and hence about the height of the energetic barrier at their diabatic crossing points. 8 To summarize, the Marcus rate (3.18) is a simple yet powerful analytical rate expression for charge transfer that depends on three parameters of the charge transfer complex. First, the transfer integral J ij (3.8), specifying the electronic coupling between the two diabatic electronic states, in which the excess charge is localized at the molecules i and j,respectively.econd,thesiteenergydifference E ij (3.19), denoting the difference between the energies of these two electronic states with molecules being in their nuclear equilibrium configurations. And third, the reorganization energy λ ij (3.21), accounting for the energetic barrier between these nuclear equilibrium configurations, which needs to be overcome by nuclear vibrational motion. Note that the latter is also referred to as electron-phonon coupling. Depending on the relation of the site energy difference and the reorganization energy, one distinguishes between two regimes of the Marcus rate. If the site energy difference is smaller than the reorganization energy ( E ij < λ ij ), the rate is an increasing function with the site energy difference, which is denoted as the normal regime. However, if the site energy difference exceeds the reorganization energy ( E ij > λ ij ), the rate turns into a decreasing function with increasingly favorable driving force, which is called the inverted regime [83]. This is because high site energy differences lead to the formation of final states with high vibrational energy, which is experimentally established [84]. In Figure 3.3, the inverted regime corresponds to situations where the lower potentials cross the energetically higher ones on the left, and no longer on the right side with respect to their minima. 8 Note that the values have no direct physical meaning; rather, one can imagine a process in which a molecule in its nuclear equilibrium configuration of the initial electronic state (black dot) experiences, by violation of energy conservation, a vertical transition into the final electronic state (gray dot) and then reorganizes its nuclear coordinates to the new equilibrium configuration (red dot) under dissipation of the reorganization energy.

69 68 Chapter 3. Charge Transfer in Molecular ystems Further Limits of Nonadiabatic Charge Transfer The previously derived Marcus rate of nonadiabatic charge transfer covers bimolecular reactions under the constraints of a static environment and high temperature. To highlight the resulting scope of the Marcus rate, a comparison to closely related expressions is instructive. To this end, two prominent examples are presented: the more general Marcus-Levich-Jortner rate and the simpler Miller-Abrahams rate. AgeneralizationoftheMarcusrateshould,ontheonehand,accountforthenuclear relaxation of the environment associated with the charge transfer reaction. This is achieved by introducing outer-sphere normal mode coordinates, which are common to both molecules and are therefore treated by a single outer-sphere reaction coordinate, coupled to the internal ones. The outer-sphere reaction coordinate captures comparatively slow reorientations of molecules surrounding the charge transfer complex, hence its average eigenfrequency ω vib out is much lower than the eigenfrequencies of the internal reaction coordinates, ω vib i and ω vib j.ontheotherhand,ageneralizationmay embrace medium temperatures, whereinternalnuclearvibrationshavelarge energy quanta compared to thermal energy (ħω vib i, j k B T) and are therefore treated quantum mechanically. The nuclear dynamics related to the low-frequency outer-sphere reaction coordinate may, however, still be described classically (ħω vib out k B T). This situation implies that in the initial electronic state, the molecules i and j are restricted to their vibrational ground states with quantum numbers 0, while the final electronic state can be formed in a manifold of vibrational states with quantum numbers η and ϑ, respectively. Consequently, Fermi s golden rule comprises two summations over quantized manifolds of vibrational states and one integral over a reaction coordinate distribution function, such as (3.16), describing the continuous outer-sphere energy. Carrying out the integration over the outer-sphere reaction coordinate [4, 85, 86], which ensures energy conservation, yields a generalized (to the bimolecular case) Marcus-Levich-Jortner Rate ω ij = 2π ħ J 2 ij χ 0 i 4πλ out ij k B T χη i 2 χ 0 j χϑ j 2 η,ϑ=0 exp ( E ij ħ(ηω vib i + ϑω vib j ) λ out 4λ out ij k B T ij ) 2. (3.22)

70 3.2. Charge Transfer Rates 69 Here, the transition between the initial and final electronic states is associated with a nuclear rearrangement of the environment, captured by an outer-sphere reorganization energy λ out ij,andwithnuclearcouplingterms formoleculei (and similarly for j), relating the nuclear ground state wave function χ 0 i in the initial electronic state with the nuclear wave functions χ η i in the final electronic state, known as Franck-Condon factors: χ 0 i χη i 2 = 1 η! ( λ i η ħω vib ) i exp ( λ i ħω vib i ). (3.23) The Marcus-Levich-Jortner rate (3.22) should be preferred to the Marcus rate (3.18) if the classical treatment of the internal nuclear vibrations is no longer appropriate. In fact, internal normal mode frequencies are often comparable to the frequency of carbon-carbon bond stretching with vibrational energy quanta ħω vib cc 0.2 ev, which is nearly an order of magnitude larger than thermal energy at room temperature [63]. However, it has been shown for these cases that the Marcus rate still yields quantitatively comparable results, as long as the site energy difference is not large compared to the internal reorganization energy, i.e., charge transfer does not occur far outside of the normal regime [4]. If the nuclear relaxation of the environment remains to be taken into account, one can extend the Marcus rate by adding the outer-sphere reorganization energy λ out ij to the internal reorganization energy λ ij given in (3.21). However, in solid materials, contrary to polar solvents, the outer-sphere reorganization energy is typically only a small contribution [87 89]. For these reasons, the studies of this work can rely on the Marcus rate. AsimplerexpressionthantheMarcusratearisesinthelimitoflow temperatures with weak electron-phonon coupling. This comprises the neglect of nuclear reorganization towards a nuclear transition state, where the charge transfer reaction actually takes place, as well as the sole consideration of vibrational ground states, which are bridged exclusively by tunneling. Note that energy conservation of the reaction is considered as satisfied by the absorption or emission of a single lattice phonon. In this limiting case, time-dependent perturbation theory leads to a rate [90] known as the Miller-Abrahams Rate ω ij = ω 0 exp ( 2γ ij r ij ) exp ( E ij k B T ), for E ij > 0, 1, for E ij 0. (3.24)

71 70 Chapter 3. Charge Transfer in Molecular ystems Here, the electronic coupling is expressed as J 2 ij = exp( 2γ ijr ij ),whereγ ij denotes the inverse localization radius of the excess charge and r ij is the distance between its localization centers, i.e., between the donor and acceptor molecules. The prefactor ω 0, referred to as the hopping attempt frequency, can be set as the ground state vibrational frequency or be considered as a material constant. For fixed electronic coupling, the Miller-Abrahams rate decreases exponentially for a transfer reaction upwards in the energy ( E ij > 0), while a reaction downwards in energy ( E ij 0) is independent of the site energy difference. Contrary to the Marcus rate, the Miller-Abrahams rate lacks aparameterforthereorganizationenergy.hence,itprecludesaquantitativecalculation of charge transfer rates at ambient temperatures and is therefore unsuitable for this work. Nonetheless, it is often used in conjunction with phenomenological charge transport models, such as the Gaussian disorder models (ection 4.1.2).

72 Chapter 4. Charge Transport in Organic olids The electrical conductivity of organic solids depends on their ability to support charge carrier migration over macroscopic distances, that is, charge transport [6, 7, 91 96]. The efficiency of charge transport is specified by the charge carrier mobility, defined as the drift velocity attained by a charge carrier per unit electric field applied. The charge carrier mobility in organic materials depends significantly on the electronic structure of the constituent molecules, but no less critically on their packing scheme on a large scale. In addition, for a given material sample, the mobility is influenced by further parameters, such as the applied electric field and the temperature. While the molecular packing of the material, being somewhere in between perfectly crystalline and completely amorphous, is related to so-called static disorder, the temperature of the system gives rise to dynamic disorder. Theabsoluteandrelativestrengthsofthese types of disorder are often employed to distinguish between several regimes of charge transport. Within these regimes, charge carrier migration differs in its very nature, which is a reason that, as of today, there is no unified theory of charge transport for organic solids, but instead a variety of different theories. Many of these charge transport theories are based on charge transfer as an elementary process. While some theories lead to generic, material-independent expressions of the mobility as a function of, for example, the electric field and temperature, others further aim at predicting mobility magnitudes for concrete material morphologies. The following discussion starts with an overview of charge transport regimes and the related transport theories (ection 4.1). For their classification, a general charge transport Hamiltonian is introduced which formalizes both static and dynamic disorder. Then, particular attention is directed to a predictive charge transport model which includes both types of disorder and is the core methodology of this work (ection 4.2). To parametrize the charge transport Hamiltonian, this model relies on the extensive use of computational chemistry methods, all of which are presented in detail.

73 72 Chapter 4. Charge Transport in Organic olids 4.1. Regimes of Charge Transport From a microscopic point of view, charge transport of a single excess charge in a largescale organic solid can be described by a generalization of the two-state charge transfer Hamiltonian, such as (3.7), to a large basis of electronic states. In the diabatic representation, these are localized states ψ i with non-vanishing amplitude at spatial localization centers, or sites i, and the charge transport Hamiltonian can be written as Ĥ = Ĥ 1 + Ĥ 2 + Ĥ 3 + Ĥ 4, where (4.1) Ĥ 1 = i Ĥ 2 = ξ E i ψ i ψ i + J ij ψ i ψ j, (4.2) i, j 1 2 (P2 ξ +(ωvib ξ Q ξ ) 2 ), (4.3) Ĥ 3 = λ 1/2 iξ ωvib ξ Q ξ ψ i ψ i, i,ξ (4.4) Ĥ 4 = Λ 1/2 ijξ ωvib ξ Q ξ ψ i ψ j. i, j,ξ (4.5) Ĥ 1 is the equilibrium energy of the system, where all nuclear positions are fixed at their energetic minimum. This energy is determined by the site energies E i and transfer integrals J ij. Ĥ 2 represents the nuclear vibrational energy in terms of phonons with frequencies ω vib ξ, normal mode coordinates Q ξ,andconjugatemomentumoperatorsp ξ. 1 This term collectively describes low-frequency intermolecular vibrations, i.e., mutual displacements of entire molecules, and high-frequency intramolecular vibrations, and that by means of reference oscillators which are centered at the origin in the space of normal mode coordinates. Ĥ 3,knownaslocal electron-phonon coupling, accounts (in linear approximation) for modulations of the site energies E i due to intramolecular vibrations and reflects the offset of intramolecular oscillators from the origin. 2 The coupling constants λ iξ are thus associated with the reorganization energies (3.21), generally by a relationship λ i ξ λ iξ,wherethesumvanishesundertheassumption that each site is only coupled to a single phonon branch. Ĥ 4,referredtoasnonlocal electron-phonon coupling, describes modulations of the transfer integrals J ij due to lowfrequency intermolecular vibrations with the coupling constants Λ ijξ. 1 Alternatively, a conjugate pair of dimensionless operators C ξ and C ξ for the creation and annihilation of phonons with energy ħω vib ξ may be introduced, as usual, by C ξ = Q ξ (ω vib ξ /2ħ) 1/2 + ip ξ (1/2ħω vib ξ ) 1/2. 2 Ĥ 1 + Ĥ2 + Ĥ3 corresponds to the two-state Hamiltonian (3.7), where the nuclear Hamiltonians (3.14) are linearly expanded in the offset from the nuclear equilibrium configurations Q c,n i, j.

74 4.1. Regimes of Charge Transport 73 (a) Band Transport λ J (b) Polaron Transport λ J Figure 4.1. In ordered organic solids, the periodic crystal lattice (gray) allows for delocalization of the charge carrier (red) over multiple sites. (a) For weak dynamic disorder, i.e., λ J, the carrier rarely interacts with the lattice and migrates by band transport. (b) For strong dynamic disorder, i.e., λ J, the carrier deforms the lattice, leading to polaronic transport. The charge transport Hamiltonian Ĥ (4.1) allows the distinction between static and dynamic disorder. tatic disorder, which refers to fluctuations evolving slowly compared to the motion of a charge carrier, is described by the equilibrium Hamiltonian Ĥ 1. The variations of site energies E i and transfer integrals J ij are denoted as diagonal and off-diagonal static disorder, or as energetic disorder σ and electronic disorder Σ: σ = E 2 E 2, Σ = J 2 J 2. (4.6) Dynamic disorder, which refers to rapid fluctuations relative to charge dynamics, is captured by the local and nonlocal electron-phonon couplings Ĥ 3 and Ĥ 4.Therelated modulations of site energies and transfer integrals are denoted as diagonal and offdiagonal dynamic disorder and are quantified by the coupling constants λ and Λ. Note that the discrimination between static and dynamic disorder is often not obvious, since organic solids may comprise phononic modes which evolve slower than charge dynamics and should be treated as static disorder. This is particularly important for organic crystals, where static disorder due to frozen molecular motion is absent. Charge transport in organic solids may not only differ in the magnitude of the charge carrier mobility µ, but also in the functional dependence of µ on the temperature T or external field F, andthereforebyitsverynature.thenatureoftransport depends on the interplay of the average electronic coupling J,which enables charge migration and delocalization, and the influences of disorder, that is, σ, Σ, λ,and Λ, which may counteract these effects. The overview of charge transport regimes is divided into statically ordered (ection 4.1.1) and disordered organic solids (ection 4.1.2) Ordered Organic olids Ordered organic solids are devoid of any static disorder, that is, σ J and Σ J, including phonon modes evolving slower than charge carrier dynamics. This applies to certain molecular crystals which have no structural imperfections whatsoever.

75 74 Chapter 4. Charge Transport in Organic olids Weak Dynamic Disorder If dynamic disorder is weak, λ J and Λ J, the charge carrier leaves the periodic crystal lattice almost unaffected, as is sketched in Figure 4.1 a. This gives rise to the formation of an electronic band structure of allowed and forbidden energy regions. The charge carrier is then completely delocalized, described by Bloch states extending over the entire crystal, and charge transport is understood as band transport [97 99]. Transitions between different Bloch states are initiated by rare events of phonon scattering, i.e., perturbative lattice vibrations. As the scattering probability increases with the temperature, the charge carrier mobility is a decreasing function with temperature, µ T s,wheres > 1 depends on the type of scattering processes. The concepts of band transport, originally developed for inorganic materials, apply to organic solids only in a few cases of ultrapure single crystalline phases of rigid small molecules at low temperatures, where charge carrier mobilities of 10 1 cm 2 /Vs and higher are observed. trong Dynamic Disorder If dynamic disorder is strong, λ J and Λ J, where the latter is usually neglected, a charge carrier can induce local lattice distortions, which act back on the carrier and impede extended delocalization, as indicated in Figure 4.1 b. uch a self-trapped charge carrier together with its surrounding lattice deformation is termed a polaron [100, 101]. Charge carrier migration is then conceivable as polaron transport,whichisextensivelydiscussedbyalargenumber ofpolaron theories [ ]. At low temperatures, polaron theories typically result in charge delocalization over multiple sites, adiabatic transport, and a band-like temperature dependence, i.e., µ T s.inthiscase,experimentallow-temperaturemobilitiesoforganic crystals, typically of the order of 10 1 to 10 1 cm 2 /Vs, can often be reproduced. At ambient temperatures, it is predicted that polarons can localize on individual sites and transport occurs by nonadiabatic hopping processes, which are thermally activated, thus leading to an increasing mobility with temperature: µ exp( λ /k B T). In this case, quantitative results often disagree with experiments and even the qualitative model is inconsistent with certain room-temperature experiments, where small charge carrier mean-free paths point to strong charge localization, but a band-like temperature dependence is observed at the same time [93, 110]. Indications suggest that these inconsistencies can be resolved by accounting for the usually neglected nonlocal electron-phonon coupling Λ. A corresponding charge transport model, known as diffusion limited by dynamic disorder,hasonlybeendevelopedinrecentyears,firstin one-dimensional [111], then in two-dimensional space [112], allowing the reproduction of room-temperature mobilities in rubrene crystals of the order of 10 1 cm 2 /Vs.

76 4.1. Regimes of Charge Transport Disordered Organic olids Disordered organic solids possess significant static disorder, i.e., σ J and Σ J, which can apply to different phase states of matter. In organic crystals, static disorder can arise if nuclear dynamics is rich in time scales and low-frequency vibrations fall below the time scale of charge carrier dynamics [ ]. Further sources of static disorder are structural imperfections, such as lattice dislocations, residual impurities of foreign molecules, and, in polycrystalline materials, grain boundaries between adjoining crystallites. In mesomorphic or amorphous organic solids, where molecules are arranged irregularly, static disorder is, of course, particularly strong. Unlike in perfectly ordered semiconductors, where energetic bands are sharply delimited from the forbidden gap, static disorder leads to a broadening of the band edges, i.e., the formation of band tails extending in the forbidden gap. The tail states located in the forbidden gap correspond to strong charge carrier localization on the individual sites i,knownasandersonlocalization[116].theformationofthesetrapstatesresults inhopping transport occurring as a sequence of nonadiabatic charge transfer reactions between diabatic electronic basis states ψ i. To describe the evolution of the charge carrier in time τ, one adopts a quantum statistical picture, where each diabatic state is associated with an occupation probability p i, i.e., statistical weight, and the statistical mixture of these pure states is characterized by a density operator ρ. It is recalled that the time-dependent chrödinger equation for the pure states leads to the von Neumann equation for the time evolution of ρ: iħ ρ τ = [Ĥ, ρ ], ρ = p i ψ i ψ i. (4.7) i Under the Markov approximation of incoherent charge transfer processes, this can be cast into an equation of motion for the occupation probabilities p i [34, 117], known as Master Equation p i τ = j p j ω ji p i ω ij, i, (4.8) j where nonadiabatic charge transfer rates ω ij,associatedwithtransitionsbetweenelectronic states ψ i and ψ j, take the form of Fermi s golden rule. In order to apply explicitly derived rate expressions (ection 3.2), it is required to specify the strength of dynamic disorder in the statically disordered organic solid.

77 76 Chapter 4. Charge Transport in Organic olids (a) Hopping Transport Miller-Abrahams Regime λ J (b) Hopping Transport Marcus Regime λ J Figure 4.2. In disordered organic solids, the uneven energetic landscape (gray) leads to localization of the charge carrier (red) resulting in hopping transport. (a) For weak dynamic disorder, i.e., λ J, charge transfer is governed by tunneling, according to the Miller-Abrahams rate. (b) For strong dynamic disorder, i.e., λ J, charge transfer is activated by vibrational site energy modulations, according to the Marcus rate. Weak Dynamic Disorder For weak dynamic disorder, λ J and Λ J, and low temperatures, the master equation contains the Miller-Abrahams rate (3.24), leading to hopping transport by tunneling, as is illustrated in Figure 4.2 a. Obviously, when aiming at wide temperature ranges, this regime can only constitute phenomenological charge transport theories. uch theories are in particular the Gaussian disorder models [118], which use simple and intuitive relations to determine the parameters in the transport Hamiltonian (4.1) entering the rate expressions. First, site energies E i are drawn from a Gaussian distribution of postulated width σ. econd, a distribution of transfer integrals J ij with standard deviation Σ is determined by the isotropic expression Jij 2 = exp( 2γ ijr ij ),wheretheinversechargelocalizationradiusγ ij = γ i + γ j is composed of site-specific contributions, drawn from empirical Gaussian distributions. The inter-site distances r ij are typically obtained by arranging the hopping sites on a regular lattice. Numerical calculations solving the master equation have led to an empirical mobility function of the temperature T and electric field F [119, 120], 3 µ(t, F) =µ 0 exp { ( 2σ 3k B T ) 2 + β F (( σ k B T ) 2 Σ 2 )}, (4.9) including the static disorder σ and Σ as well as coefficients µ 0 and β as parameters. While the mobility increases with temperature, following a super-arrhenius relation, the field dependence ln µ F 1/2 is referred to as the Poole-Frenkel relationship. Within certain ranges, these relations are observed experimentally for numerous disordered organic materials [120, 124], and respective field ranges often become reproducible by imposing spatial correlations of site energies [ ]. The analytic function obtained by Gaussian disorder models (4.9) can be fitted to temperature and field-dependent mobility measurements of many materials, e.g., to characterize their static disorder, however it provides no material-specific structure-property relationship. 3 A relation for the charge carrier density is established by extended Gaussian disorder models [ ].

78 4.1. Regimes of Charge Transport 77 trong Dynamic Disorder For strong dynamic disorder, λ J and Λ J, where the off-diagonal part is typically negligible compared to its static pendant Σ, charge transfer is activated by intramolecular vibrations, as is sketched in Figure 4.2 b. Hence, at high temperatures, where vibrational motion behaves classically, the master equation (4.8) is equipped with the Marcus rate (3.18), while at medium temperatures, where quantum mechanical effects are important, the more general Marcus-Levich- Jortner rate (3.22) is applied. ince this regime, often with the Marcus rate, has its physical justification at ambient temperatures, one can formulate charge transport theories that transcend the phenomenological model of Gaussian disorder with randomly drawn transport parameters. In this spirit, microscopic charge transport models start from atomistically resolved molecular arrangements, i.e., material morphologies for specific chemical compounds [4, 130, 131]. With the knowledge of the chemical structures and atomistic morphologies, and the aid of computational chemistry methods, it becomes possible to explicitly determine the parameters of the charge transport Hamiltonian (4.1), i.e., the site energies E i,transferintegralsj ij, and reorganization energies λ ij, that enter the transfer rates. Thus, for a specific pair of sites, the transfer rate no longer depends only on the inter-site distance, but is also sensitive to the actual arrangement of the underlying molecular system. olving the master equation now allows one to predict the material-specific charge transport properties, in other words to link the chemical structure and the material morphology to the macroscopic charge carrier mobility. A realization of this microscopic modeling paradigm, establishing structure-property relationships, is the core methodology of this work (ection 4.2). Microscopic charge transport studies have already been performed in recent years for a variety of disordered organic solids. In amorphous phases of conjugated polymers [132] and small molecules [ ], experimentally measured charge carrier mobilities in the range of 10 9 to 10 2 cm 2 /Vs are reproduced with often considerable accuracy. Inmesomorphic phases of columnarly arranged discotic molecules [ ] and carbazole macrocycles [142], both theoretical and experimental mobilities lie in similar ranges of 10 2 to 10 0 cm 2 /Vs. In semicrystalline phases of conjugated polymers [143] and crystalline phases of several small molecules [ ], theoretical mobilities in the range of 10 1 to 10 1 cm 2 /Vs show reliable directional dependences. ince charge transport models typically neglect grain boundaries and impurities, they tend to systematically overestimate experimental results for crystals. In cases where experimental mobility measurements are performed on highly ordered and pure material samples, such as rubrene, mobilities of the order of 10 1 cm 2 /Vs agree very well.

79 78 Chapter 4. Charge Transport in Organic olids 4.2. Charge Transport imulations in Disordered Organic olids Charge transport models for disordered organic solids, which are aimed at establishing structure-property relationships, engage atomistic material morphologies for specific chemical compounds in order to parametrize the charge transport Hamiltonian (4.1). Once a compound of interest is selected, this concept requires computational chemistry techniques, such as ab initio methods (ection 2.2) but also, in order to overcome prohibitive computing times, semiempirical and classical approximations. In what follows, the workflow of charge transport simulations, depicted in Figure 4.3, and the employed computational methods are presented. Any charge transport study starts from the chemical structure of a particular organic semiconductor, e.g., dcv4t. On the experimental side, chemical synthesis, material processing, and structure analysis may lead to the crystal structure, provided that crystallization occurs. For amorphous substances, experiments may identify other expedient material characteristics, such as the density or glass transition temperature. On the computational side, density functional theory (dft) serves as the basis for developing an atomistic force field for the chemical structure (ection 4.2.1). Using both the force field and an experimentally inspired starting configuration, molecular dynamics simulations yield an atomistically resolved morphology of the material (ection 4.2.2). The morphology allows for the identification of charge localization centers, that is, the hopping sites, and the charge transfer complexes, which define a neighbor list of hopping site pairs (ection 4.2.3). Then follows, for each pair of hopping sites, the explicit evaluation of the charge transfer parameters. The transfer integrals J ij can be determined using the semiempirical zindo method (ection 4.2.4). The site energy differences E ij = E i E j result from the site energies E i and E j, which are derived from distributed point charges and dipole moments interacting via the Thole model (ection 4.2.5). The reorganization energies λ ij are computed by dft calculations on the charge transfer complexes (ection 4.2.6). These parameters lead to the charge transfer rates and thus to the parametrization of the master equation. olving the master equation can be performed by the kinetic Monte Carlo method and results in the time-dependent site occupation probabilities, i.e., the charge dynamics (ection 4.2.7). This information can finally be related to macroscopic observables,suchasthecharge carrier mobility (ection 4.2.8). This workflow is the subject of prior publication [4] and is presented below in more detail using partially similar terms. As part of this and many other works, the methods were implemented in the votca-ctp software package, while dft and molecular dynamics methods stem from the gaussian [53] and gromacs packages [148].

80 4.2. Charge Transport imulations in Disordered Organic olids 79 Force Field Chemical tructure Crystal tructure θ1 θ2 θ3 dft NC dcv4t NC yn., Exp. Molecular Dynamics tarting Configuration Morphology Hopping ites (and Neighbor List of Hopping ite Pairs) zindo Thole Model dft Transfer Integrals J ij ite Energies E i Reorganization Energies λ ij i j i i j Charge Transfer Rates ω ij = ω ij (J ij, E ij, λ ij ) Master Equation Kinetic Monte Carlo Charge Dynamics Macroscopic Observables (e.g., Charge Carrier Mobility) Figure 4.3. Workflow for microscopic simulations of charge transport in disordered organic solids.

81 80 Chapter 4. Charge Transport in Organic olids (a) R ab (b) (c) (d) (e) R ab θ abc θ abcd θ abcd Figure 4.4. Atomic interactions contributing to a force field. (a) Nonbonded interaction between two atoms. (b) Bond stretching between two atoms. (c) Angle bending between three atoms. (d) Dihedral torsion between four atoms. (e) Out-of-plane bending between four atoms Force Field The molecular modeling of a material morphology at an atomistic scale relies on simulating the dynamics of the atoms, or equivalently of the nuclei, which are adiabatically followed by the electrons, according to the Born-Oppenheimer (or adiabatic) approximation (2.14). To describe this nuclear or atomic motion, one needs knowledge about the potential energy surface U(R), which captures, for each atomic configuration R, the ground state energy of the electronic system. However, using ab initio methods to globally determine the pes is an infeasible endeavor for practical systems with highdimensional atomic configuration space. A simplified yet very effective alternative is to employ a force field: itreferstoananalyticalfunction,decomposedintoempirical atomic few-body contributions that are fitted to ab initio or experimental results [54]. The underlying notion is to fit the function only to the minima of the pes that correspond to atomic configurations of covalent bonding to the desired molecules, while effects of bond breaking and formation are neglected. Hence, the bonding pattern is definedapriori, which allows the distinction between few-body contributions for covalently bonded and not directly bonded atoms. Most force fields employ a two-body term for nonbonded and two-, three-, and four-body terms for bonded contributions: U ff (R) = U nb (R ab )+ U bond (R ab )+ U ang (θ abc )+ U dih (θ abcd ).(4.10) a,b nonbonded a,b bonds a,b,c angles a,b,c,d dihedrals Figure 4.4 illustrates the meaning of these few-body terms. U nb is the potential energy associated with interactions between two nonbonded atoms, which comprises the very short-range Pauli repulsion emerging from overlapping electron orbitals, the shortrange van der Waals attraction due to induced dipole-dipole interactions, and the longrange electrostatic interaction between atomic multipole moments arising from intramolecular charge redistributions. U bond describes the energy for stretching a bond from its equilibrium length, U ang for bending a bond angle formed of three serially

82 4.2. Charge Transport imulations in Disordered Organic olids 81 n n cn hc hc hc hc cn cc θ1 cb cb θ2 cb cb s ca ca s hn c ca ca ca ca s ca ca s c θ3 hn cb cb cb cb cc cn hc hc hc hc cn n n Figure 4.5. Atom types in the force field for dicyanovinyl-substituted oligothiophenes. Also indicated are the dihedral angles θ 1 (cc c ca s), θ 2 (s ca ca s), and θ 3 (cn cc c ca). Adapted from Reference [1]. linked atoms, and U dih for rotating a dihedral angle defined by four serially connected atoms. An energy penalty for out-of-plane bending defined by four radially bonded atoms, as seen in Figure 4.4 e, may be incorporated by improper dihedral angles. Many popular force fields, such as amber [149], charmm [150], gromos [151], or opls-aa [152], employ energy potentials which are essentially of the prototypical form U nb (R ab )= A ab R 12 B ab ab R 6 + q a q b, (4.11) 4πε ab 0 R ab U bond (R ab )= 1 2 C ab (R ab R0 ab )2, (4.12) U ang (θ abc )= 1 2 C abc (θ abc θ0 abc )2, (4.13) U dih (θ abcd )= 5 k=0 C abcd,k ( cos θ abcd ) k, (4.14) where the Pauli and van der Waals interactions are modeled by the Lennard-Jones potential [153], the electrostatic interaction is described by atomic partial charges, i.e., distributed monopole moments, the bond and angle terms are represented by harmonic potentials, and the dihedral torsion is expressed by a periodic function, here the Ryckaert-Belleman potential [154]. While force fields may share the same functional form, they differ in the parametrization of the potentials. uch a parametrization is based on atom types, defined not only by the atomic number, i.e., the chemical element, but also by the bonding situation, hybridization state, or local environment. This enables the identification of recurring units within different molecules, thus gives force fields a versatile applicability to various chemical compounds, but also limits their accuracy. In particular for organic small molecules, a compound-specific refinement is essential. In this work, a customized force field for dicyanovinyl-substituted oligothiophenes (dcvnt) is developed [1], which is based on the opls-aa force field and the atom types defined in Figure 4.5. Regarding the nonbonded interaction energy (4.11), the Lennard-

83 NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC 82 Chapter 4. Charge Transport in Organic olids Charge (qa/e) n cn cc c hn s ca cb hc Figure 4.6. Atomic partial charges q a for rotamers of dicyanovinyl-substituted quaterthiophene (dcv4t), obtained by the chelpg method and dft calculations using the b3lyp functional and 6-311g basis set. Jones parameters A ab and B ab are adopted from the opls-aa force field, where two atoms a and b of different types are treated by the geometric combination rules A ab =(A aa A bb ) 1/2 and B ab =(B aa B bb ) 1/2.Theatomicpartialchargesq a,whichare representing the molecular charge distribution, are specifically determined for the dcvnt molecules. Tothisend,themoleculargeometriesarefirstoptimizedusing dft calculations with the b3lyp functional and the 6-311g basis set resulting in the atomic coordinates of the stable molecules, the corresponding electron densities, and the associated electrostatic potentials. Then, point charges at the atomic coordinates are fitted to reproduce the molecular electrostatic potentials on cubic grids according to the chelpg method [155, 156]. Finally, the partial charges are averaged for atoms of the same type and adjusted to ensure charge neutrality of the dcv and thiophene units [1], resulting in the values listed in Table 4.1. To verify that their validity is not significantly affected by conformational variations, the procedure is repeated for rotamers of the dcv4t molecule, defined by torsions around the bonds linking dcv and thiophene units. Figure 4.6 shows that the dcv4t rotamers have very similar partial charges. Concerning the bonded interaction energies (4.12) to (4.14), all equilibrium bond lengths R 0 ab and bond angles θ0 abc are obtained from the optimized molecular n cn cc c hn s ca cb hc q a Table 4.1. Atomic partial charges q a for dicyanovinyl-substituted oligothiophenes (dcvnt), obtained by the chelpg method and dft calculations using the b3lyp functional and 6-311g basis set. The values are given in units of e.

84 4.2. Charge Transport imulations in Disordered Organic olids 83 Potential Energy / ev U dft dih U ff dih U ff dih Dihedral Angle (θ 1, θ 2 / ) Figure 4.7. Potential energy for the dihedrals θ 1 and θ 2, from density functional theory, Udih dft, and the initial and final force fields, Udih ff and Udih,respectively[1]. ff θ 1 θ 2 θ 1 θ 2 θ 3 C C C C C C Table 4.2. Ryckaert-Belleman parameters C k for the dcvnt dihedral angles θ 1, θ 2,andθ 3, as defined in Figure 4.5. The values are given in units of ev. geometries, while the force constants for bond stretching C ab and angle bending C abc are adopted from the opls-aa force field. This is justified since these degrees of freedom are fairly rigid, so that only small deviations from the equilibrium values are expected. The dihedral angles, in contrast, have often comparatively low torsional energy barriers and therefore a major impact on molecular conformations. For this reason, the Ryckaert-Belleman parameters C abcd,k are specifically determined for the dcvnt molecules. It is essential that such a parametrization takes into account the existing nonbonded interactions as they are intimately coupled to a torsional potential. Therefore, constrained geometry optimizations for a set of fixed dihedral values need to be performed by both density functional theory and also the tentative force field, which yields functions Udih dft and Uff dih,respectively. Then,theRyckaert-Bellemanpotential U dih,givenin(4.14),canbefittedtothedifferenceudih dft Uff dih,suchthatthefinal force field results in a function Udih ff that resembles the dft curve. This procedure is performed for the three dihedral angles indicated in Figure 4.5, θ 1 (cc c ca s), θ 2 (s ca ca s), and θ 3 (cn cc c ca), using dft calculations with the b3lyp functional and the 6-311g basis set. The obtained Ryckaert-Belleman parameters are listed in Table 4.2, while the aforementioned functions are depicted in Figure 4.7 for the two dihedrals θ 1 and θ 2,whichhavecomparativelylowenergies.Thepotentialenergy of the dcv-thiophene dihedral θ 1 has a global minimum at 0, which corresponds to the cis configuration, while the thiophene-thiophene dihedral θ 2 prefers the trans configuration [157, 158]. Both potentials have local minima in the opposite state, separated by energetic barriers which can be overcome at elevated temperatures.

85 84 Chapter 4. Charge Transport in Organic olids Morphology In a microscopic charge transport model, the transport Hamiltonian is parametrized based on an atomistic material morphology of the organic solid of interest. Particularly important is that the morphology provides a realistic picture of the disorder inherent to the material. Hence, to capture the related characteristic length scales and achieve adequate statistics, system sizes much larger than a unit cell are required. For the generation of such mesoscopic morphologies, ab initio methods are generally too complex; instead, methods of choice are Monte Carlo [159, 160] or, as used in this work, molecular dynamics simulations [161, 162]. Molecular Dynamics imulation The classical molecular dynamics simulation refers to the study of atomic movements based on a potential energy surface in force field representation (4.10) and an equation of motion (2.14) in classical approximation, which is the Newtonian equation of motion (2.15). This differential equation in time t is numerically solved, subject to appropriate boundary conditions, typically of periodic nature, and initial conditions for the atomic coordinates R a and velocities V a = R a / t. The numerical integration is performed iteratively and over discrete time steps: at the futuretime t + t the coordinates are expressed by the Taylor expansion R a (t + t) R a (t) +( R a t R a t ) t, (4.15) t2 where the term in the bracket can be identified with a Taylor expansion of the velocities at time t t yielding the relation (4.17). Performing a similar expansion of the velocities at time t 1 2 t and using the Newtonian law (2.15) to express the accelerations by forces, i.e., gradients of the force field potential, leads to an integration scheme V a (t t) = V a(t 1 2 t) au ff (R) M a t, (4.16) R a (t + t) = R a (t) + V a (t t) t, (4.17) known as theleap-frog algorithm [163]. In this scheme, velocities and coordinates are alternately evaluated at time points that are mutually shifted by half a time step. In one step, previous velocities and current forces are used to calculate new velocities, and in the next step, these new velocities are used to determine new coordinates. If velocities at the integer time points are required, a third evaluation step can be included [164].

86 4.2. Charge Transport imulations in Disordered Organic olids 85 An atomistic trajectory generated in this way describes a point sequence in phase space, and as such samples microstates of a statistical ensemble. The natural ensemble which is maintained when solving the energy-conserving Newtonian equation is the microcanonical ensemble, which corresponds to systems with constant number of particles, volume, and energy, i.e., isolated systems. In experiments, however, systems are usually in contact with an environment, which, as the time evolves, leads to fluctuations in some of these extensive variables and to conservation of their respective intensive variables, i.e., chemical potential, pressure, or temperature. In other words, the environment can act as chemostat, barostat, or thermostat. In fact, most experiments are performed under conditions of constant particle number, temperature, and either volume or pressure. To conduct corresponding simulations, the leap-frog integration scheme needs to be adapted to simulate a thermostat, which appropriately rescales the velocities, i.e., the kinetic energies, thus leading to a canonical ensemble, or an additional barostat, which also rescales the coordinates, i.e., the box volume, leading to an isobaric-isothermal ensemble. In this work, a stochastic velocity-rescaling thermostat [165] and a Berendsen barostat [166] are employed. The implementation of the adjusted leap-frog algorithm is provided by the gromacs package [148]. Any molecular dynamics integration scheme is afflicted with an error resulting from the discretization of time. In order to keep this error small, the time step t must be significantly shorter than the vibration period of the fastest processes in the system. ince carbon-hydrogen bond stretching vibrations have a period of roughly 10 fs, and carbon-carbon bonds still have about 20 fs, the time step is typically set in the order of 1 fs. This time interval limits the real time a simulation can achieve, given a certain computing power. When performing a number of 10 6 time steps, which is perfectly feasible for mesoscopic systems these days, a real time of 1 ns is reached. 4 tarting Configuration In order to integrate the equation of motion, initial conditions for the atomic coordinates and velocities are required. In principle, one might consider a random point in phase space, which complies with the fixed ensemble variables, and rely on the simulation to drive the conjugate thermodynamic variables to their (fluctuating) equilibrium values and the thermodynamic potential associated with the ensemble to its global minimum. However, when simulating molecular solids at fixed temperatures, one is in a regime, where the thermodynamic potential, i.e., the 4 If such time scales are insufficient to equilibrate the atomistic morphologies of interest, coarse-graining techniques can be employed, provided they allow for back-mapping to the atomistic resolution [167].

87 86 Chapter 4. Charge Transport in Organic olids Figure 4.8. Crystalline morphology of dcv4t molecules at room temperature, obtained by a molecular dynamics simulation. free energy surface, is landscape-dominated by the potential energy surface U ff and has highly complex energy barriers on many scales. This implies that an equilibration to the global minimum is impeded by both short real times in the nanosecond range and also insufficient kinetic energy for barrier crossing. In other words, ergodicity is violated and the region of the configuration space, which is sampled during the simulation, is critically determined by the initial conditions. While the starting velocities are typically straightforwardly assigned by a Maxwell-Boltzmann distribution, the starting configuration must be carefully prepared in view of the desired morphology. For molecular crystalline morphologies, which typically represent the global energy minimum, one depends on experimental x-ray scattering results providing lattice spacings, or ideally the exact atomic-scale crystal structure. 5 A starting configuration is then prepared as a supercell of perfectly arranged molecules on a lattice. For mesomorphic or amorphous systems, which are kinetically arrested in a metastable nonequilibrium state, it is customary to perform a preceding simulation at elevated temperatures, where a liquid crystalline or liquid phase is adopted. Representative starting configurations can then be obtained by annealing the system to lower temperatures. With a starting configuration at hand, the molecular dynamics simulation at the desired temperature is performed, comprising an equilibration and a subsequent production run. As an example, Figure 6.2 depicts a snapshot of a dcv4t organic crystal. The trajectories generated in a production run can be employed to investigate the time scales of static and dynamic disorder. 5 Alternatively, crystal structure prediction may be facilitated by advanced free energy sampling techniques, such as metadynamics [168], or Monte-Carlo-based umbrella sampling [169, 170].

88 4.2. Charge Transport imulations in Disordered Organic olids Hopping ites The morphology allows the construction of the charge transport Hamiltonian (4.1), Ĥ 1 + Ĥ 2 + Ĥ 3,andtheparametrizationofthemasterequationforchargecarrierdynamics (4.8). The first step is the identification of hopping sites for the charge carrier, in other words subsystems of the morphology where the diabatic electronic states are spatially localized. Typically, the morphology is partitioned into individual conjugated molecules or, more generally, into conjugated segments on which charge localization is expected by physical intuition. 6 The conjugated segments can often be further divided into relatively rigid, planar π-conjugated systems, which are referred toas rigid fragments. The bonded degrees of freedom linking these fragments often evolve on time scales much slower than charge carrier dynamics. In some cases, e.g., mesomorphic or amorphous systems, these slow motions are frozen due to nonbonded interactions with the surrounding molecules. Within a rigid fragment, in contrast, motions of bonded interactions, such as bond length vibrations, are often much faster than the dynamics of charge carriers. Therefore, rigid fragments can help to separate the effects of static and dynamic disorder. To this end, the rigid fragments in the morphology are replaced by corresponding geometry-optimized equilibrium copies, obtained from ab initio calculations; their alignment is achieved by matching the centers of mass and gyration tensors. This enables the separate construction of the equilibrium Hamiltonian Ĥ 1,i.e.,thetransfer integrals and the site energies, on the one hand and the dynamic contribution Ĥ 2 +Ĥ 3, i.e., the reorganization energies, on the other hand. To illustrate the concept of conjugated segments and rigid fragments, three representative molecular architectures are presented. Figure 4.9 a depicts the first: a benzothienobenzothiophene molecule consisting of three benzene and two thiophene units fused to a rigid, planar π-conjugated structure. This molecule is both a single conjugated segment and a rigid fragment. Figure 4.9 b shows a compound relevant for this work: a dicyanovinyl-substituted quaterthiophene (dcv4t), i.e., a short π-conjugated oligomer comprising four thiophene and two dcv moieties. ince a charge carrier is delocalized over the whole oligomer, the molecule is again a single conjugated segment. However, while the individual thiophene and dcv units are relatively rigid, the dihedral angles in between can be reoriented at ambient temperatures, as results from the previously developed force field (ection 4.2.1). Hence, each of the six units is a 6 In principle, more fundamental deductive approaches for obtaining diabatic electronic states may assist the partitioning of the morphology [171].

89 88 Chapter 4. Charge Transport in Organic olids (a) (b) (c) Rigid Fragment (rf) rf rf rf rf rf rf = Conjugated egment Conjugated egment rf rf rf rf rf rf rf rf rf rf rf rf Conjugated egment Conjugated egment Conjugated egment Figure 4.9. The concept of conjugated segments and rigid fragments. Conjugated segments are indicated by red and gray, rigid fragments by blue and cyan colors. (a) A benzothienobenzothiophene molecule is both a single conjugated segment and a rigid fragment. (b) A dicyanovinyl-substituted quaterthiophene molecule is a conjugated segment comprising six rigid fragments. (c) A polythiophene molecule can consist of multiple conjugated segments, while each repeat unit is a rigid fragment. separate rigid fragment. Figure 4.9 c illustrates a more general example, a long conjugated polymer, such as polythiophene. In this case, one molecule consists of multiple conjugated segments, since the π-conjugation along the polymer backbone can be broken due to large out-of-plane twists between adjacent repeat units. The partitioning of a molecule on individual conjugated segments can be performed by empirical criteria [ ], such as the dihedral angle [132]. Neighbor List of Hopping ites Having determined the hopping sites, represented by the molecules or their conjugated segments, the next step is the generation of a list of selected pairs of hopping sites which are sufficiently close to form a charge transfer complex. To this end, a distance criterion is employed which is based on the centers of mass of the rigid fragments: apairofmolecules,orconjugatedsegments,is included in the neighbor list if the distance between the centers of mass of any pair of mutual rigid fragments is below a certain cutoff. In this way, pairs are selected based on their minimum distance, rather than their center-of-mass distance, which is useful for molecules with anisotropic shapes. For each pair of the neighbor list, the charge transfer parameters and rates are evaluated, as discussed below.

90 4.2. Charge Transport imulations in Disordered Organic olids Transfer Integrals Apairofmoleculesi and j forming a charge transfer complex provides two localization centers for the charge carrier, which are associated with the diabatic electronic states ψ i and ψ j. In this diabatic representation (ection 3.1.2), the transfer integral (3.8) describes the electronic inter-state coupling and is mediated by the electronic Hamiltonian operator Ĥ el of the charge transfer complex: J ij = ψ i Ĥ el ψ j. (4.18) To determine the transfer integral, electronic structure methods are employed to acquire information about the three objects: the Hamiltonian of the molecular pair and the two diabatic states localized on the individual molecules. If the diabatic states are expressed by one-electron wave functions ϕ m i (r 1 )= r 1 ϕ m i and ϕ m j (r 1)= r 1 ϕ m j, i.e., molecular orbitals m, the frozen core approximation can be employed to argue that the charge transfer is accomplished exclusively by the frontier molecular orbitals f: ϕ f i and ϕf j.uchafrontiermolecularorbitalreferstothehighestoccupiedmolecular orbital (homo) in the case of hole transfer, and to the lowest unoccupied molecular orbital (lumo) in the case of electron transfer. Under this assumption, the transfer integral (4.18) takes the simplified form J ij = ϕ f i Ĥ ϕf j, (4.19) where Ĥ is an effective one-particle Hamiltonian operator of the pair of molecules, such as a Fock or Kohn-ham operator, introduced before (ection 2.2) as Ĥ hf (2.26) and Ĥ ks (2.36), respectively. Note that the transfer integral is very sensitive to the arrangement of the pair of molecules within the morphology, that is, to their distance and mutual orientation. As a consequence, the electronic disorder Σ (4.6), resulting from the distribution of transfer integrals for theneighborlist,istypicallyseveralorders of magnitude large. Dipro Method The two diabatic electronic states required for calculating a specific transfer integral can be obtained by performing ab initio calculations subject to electron density constraints [171]. A common approach is to strictly constrain the density to the individual molecules by simply isolating the molecules and carrying out two separate ab initio calculations. Then, the adiabatic electronic ground states

91 90 Chapter 4. Charge Transport in Organic olids of the isolated systems are identified with the diabatic states of the charge transfer complex and yield in particular the desired frontier molecular orbitals ϕ f i and ϕf j. To determine the Hamiltonian operator Ĥ, one can perform a third,unconstrained ab initio calculation on the molecular pair. With the resulting molecular orbitals ϕ n and corresponding orbital energies ε n, the one-particle Hamiltonian operator can be written in its spectral representation: Ĥ = ε n ϕ n ϕ n. (4.20) n With these results, the definition of the transfer integral (4.19) immediately becomes J ij = ϕ f i Ĥ ϕf j = n ε n ϕ f i ϕ n ϕ n ϕ f j. (4.21) Hence, the transfer integral follows from projecting the frontier molecular orbitals ϕ f i and ϕf j of the monomers onto the molecular orbitals ϕn of the dimer. This procedure, denoted as the dimer projection method (dipro), is typically realized within the framework of Hartree-Fock [176] or density functional theory [177]. To calculate the transfer integrals for the entire neighbor list, an extensive series of ab initio calculations on all individual and all pairs of molecules is required [ ]. uch a brute-force ab initio scheme is computationally demanding and may not be practicable. An alternative approach, offering a compromise between quantitative accuracy and computational effort, is the zindo method. Zindo Method An approximate method for computing a transfer integral [176] can be formulated within Zerner s intermediate neglect of differential overlap (zindo), a semiempirical electronic structure theory (ection 2.2.3). This approximate method eliminates the need for a self-consistent calculation on the molecular pair, such that only self-consistent calculations on the individual molecules i and j are required. Performing these two calculations leads, similarly to the dipro method, to the molecular orbitals of the individual molecules, ϕ m i = Ci τm φ τ i, ϕ m j = τ τ C τm j φ τ j, (4.22) which includes in particular the required frontier molecular orbitals ϕ f i and ϕf j. Here, the molecular orbitals are expanded in respective basis sets of atomic orbitals

92 4.2. Charge Transport imulations in Disordered Organic olids 91 { φ τ i } and { φτ j }, according to the previous definition (2.40). Now, before constructing the Hamiltonian operator of the pair, two preliminary considerations are necessary. First, the basis set for the pair of molecules is given by the joint basis set { φ τ i } { φτ j }, where the atomic orbitals from the individual molecules are numbered sequentially. Hence, in this joint basis set, the molecular orbitals of the individual molecules (4.22) are represented by zero-extended column vectors C m i = ( C 1m i,...,c (n 2 )m i, 0,..., 0 ) T, C m j = ( 0,..., 0, C (n 2 +1)m j,...,c nm j ) T. (4.23) The second aspect concerns the unknown molecular orbitals of the pair of molecules: each pair of molecular orbitals of this combined system, for instance ϕ n and ϕ n+1, is assumed to be formed from a bonding-antibonding combination of a pair of molecular orbitals of the individual molecules [176]: ϕ n = ϕ m i + ϕ m j, ϕn+1 = ϕ m i ϕ m j. (4.24) This is reasonable, since the pair of molecules is expected to be characterized by two charge localization centers on the individual molecules. The implication is that the molecular orbitals of the pair are represented by column vectors C n C (n+1) = ( C 1m i,...,c (n 2 )m i = ( C 1m i,...,c (n 2 )m i, C (n 2 +1)m j, C (n 2 +1)m j,..., C nm j ) T,,..., C nm j ) T. (4.25) Now, inserting the frontier molecular orbitals of the individual molecules ϕ f i and ϕf j, expressed as linear combination of atomic orbitals (4.22), into the definition of the transfer integral (4.19) leads to J ij = ϕ f i Ĥ ϕf j = n 2 σ=1 n Ci σf τ= n 2 +1 C τf j φ σ i Ĥ φτ j, (4.26) where the one-particle Hamiltonian operator Ĥ of the pair of molecules appears in its matrix representation in the joint basis set of atomic orbitals: H στ = φ σ i Ĥ φτ j. Obviously, due to definition (4.23), reflected in the summation indices, only the off-

93 92 Chapter 4. Charge Transport in Organic olids Figure Highest occupied molecular orbital of a dcv4t molecule, depicted as an isosurface of the probability density with red and gray colors representing the positive and negative phases of the wave function, respectively. The orbital extends phasealternating over the entire molecule. diagonal elements σ τ need to be evaluated for calculating the transfer integral. According to the semiempirical Fock matrix provided by the zindo technique [59, 60], off-diagonal matrix elements associated with atomic orbitals φ σ i and φτ j,locatedon different atoms a and b, whichisnecessarilythecasehere,sincetheatomsbelongto the different molecules a i and b j,taketheform H στ zindo = 1 2 (β a + β b ) στ Γ ab C σn C τn. (4.27) The values β a and β b are tabulated bonding parameters for the atoms a and b,which depend on their atomic number, the atomic orbital overlap matrix στ contains specifically weighted overlap terms that are related to the ordinary terms στ = φ σ i φτ j, and Γ ab is the Mataga-Nishimoto potential [183]. The key insight regarding the zindo Hamiltonian matrix (4.27) is that the sum over the expansion coefficients of the molecular orbitals of the molecular pair, i.e., the right part, vanishes due to the specification of the molecular orbitals according to (4.25). This implies that the zindo Hamiltonian matrix is no longer dependent on its own eigenvectors, which enables its computation without any iterative, self-consistent procedure, as is usually required. The remaining task for calculating the required matrix elements is the comparatively facile determination of the weighted atomic orbital overlap terms, according to the parametrization by the zindo technique. To calculate the transfer integrals for the entire neighbor list, the zindo method significantly increases computational efficiency since self-consistent calculations are only required for the individual molecules, while the calculations for the pairs, which are computationally most demanding, are avoided. However, carrying out self-consistent n

94 4.2. Charge Transport imulations in Disordered Organic olids 93 Transfer Integral (Jij/eV) (a) x dipro/homo zindo/homo dipro/lumo zindo/lumo Distance (x/nm) Axial Displacement (z/nm) (b) z Figure quared transfer integrals for a pair of dcv4t molecules in cofacial alignment, calculated using the zindo and dipro methods for the highest occupied molecular orbital (homo) and the lowest unoccupied molecular orbital (lumo). (a) Transfer integral as a function of the intermolecular distance. (b) Transfer integral at a distance of 0.35 nm as function of an axial displacement. calculations for all individual molecules may still pose a serious bottleneck in charge transport simulations. Therefore, a further simplification is employed: calculations are only performed once for each chemically different type of molecule (or conjugated segment). Then, their division into rigid fragments is utilized to subject the resulting molecular orbitals to rotational coordinate transformations leading to the actual molecular conformations within the morphology. A comparison of the zindo and the dipro methods is drawn for a charge transfer complex consisting of a pair of dcv4t molecules in cofacial alignment and with varying mutual molecular positions. The ab initio calculations, performed for the dipro method, are based on the b3lyp hybrid functional and the 6-311g basis set. As an example, Figure 4.10 visualizes the homo of a dcv4t molecule. Figure 4.11 a depicts the squared transfer integral Jij 2 for the homo and the lumo as a function of the intermolecular distance ranging from 0.3 to 0.4 nm. It is seen that both the zindo and the dipro methods lead to an exponential decrease of the transfer integral with increasing intermolecular distance. While this behavior is expected and qualitatively obtained by both methods, the zindo method quantitatively underestimates the squared transfer integral by up to a factor of two, which is a well-known observation [132, 177]. Figure 4.11 b illustrates the squared transfer integral at an intermolecular distance of 0.35 nm as a function of an axial displacement ranging from 0 to 1 nm. Again, both methods lead to qualitatively similar behavior, in this case oscillations reflecting the alternating in and out of phase overlap of the frontier molecular orbitals.

95 94 Chapter 4. Charge Transport in Organic olids ite Energies The site energy E i is the energy which is absorbed or released due to the localization of the charge carrier on the molecule i. According to charge transfer models (ection 3.2), the site energy can be decomposed into internal and electrostatic parts (3.19), which refer to the molecule itself and to interactions with the environment, respectively: E i = E int i + E elstat i. (4.28) The distribution of site energies for all molecules in the morphology is usually a normal distribution of Gaussian shape whose width defines the energetic disorder σ (4.6). Internal ite Energy The internal site energy is the total energy of the molecule i in its charged state (c) with respect to its neutral state (n), as defined in (3.20): E int i = U c i U n i. (4.29) As such, the internal site energy corresponds to the molecular ionization potential in the case of hole transfer, and to the electron affinity in the case of electron transfer. The contributions U c i and U n i can be estimated by ab initio calculations with geometry optimization. Note that for one-component systems, internal site energies often need not be evaluated since they cancel in the site energy differences, provided conformational variations of the molecules are negligible, and elsewise they are often substantially smaller than the electrostatic site energy [4, 132]. In donor-acceptor systems, internal site energies should be treated with particular attention. Electrostatic ite Energy The electrostatic site energy is the difference in electrostatic interaction energies of the charged and neutral molecule i with its neutral environment, as defined in (3.20): E elstat i = W c i W n i. (4.30) Due to the long-range nature of electrostatic interactions, the energies W c i and W n i have to be determined based on a sufficiently large volume surrounding the molecule i. In this work, a spherical surrounding volume of radius 3.5 nm is employed. Within this volume, electrostatic interactions are, in analogy to the force field (ection 4.2.1), calculated between atomic partial charges q c,n a on atoms a i of the charged or neutral

96 4.2. Charge Transport imulations in Disordered Organic olids 95 central molecule i and atoms a i of the neutral surrounding molecules. In addition to this Coulomb energy, it is essential to account for the induction energy resulting from molecular polarization. To this end, the molecules are additionally equipped with atomic polarizabilities α c,n a,whichallowfortheinductionofatomic dipole moments p a. Note that a simple screened Coulomb potential containing the relative permittivity of the macroscopic material is not adequate on a microscopic scale since energetic contributions from the first coordination shell are then underestimated. uppose the atomic polarizabilities and induced dipole moments are all known, the electrostatic interaction energies follow from the charge-charge and charge-dipole interactions, W c,n i = j a j k j b k (T ab q a q b + T µ ab q a pµ b + T µ ab pµ a q b ), (4.31) with q a = q c,n a q n a and α a = α c,n a, for a i, α n a, for a i, (4.32) where vector components are marked by superscripts, i.e., p µ a is the µ-component of p a, for which summations will be implicitly assumed (Einstein s summation convention). The summations over molecules j and k and their atoms a j and b k are explicitly given and refer to atomic partial charges q a and polarizabilities α a of the charged or neutral molecule i and of the neutral surrounding molecules (4.32). The interaction tensors for two multipoles on atoms a and b with separation R ab are expressed by the Thole Model T ab = 1 4πε 0 R ab, T µ ab = κ 3 R µ ab 4πε 0 R 3 ab 3R µ abr ν ab, T µν ab = κ 5 4πε 0 R 5 ab δ µν κ 3 4πε 0 R 3, (4.33) ab where damping coefficients κ 3 and κ 5 smoothen the potential of the dipole moments in order to prevent polarization catastrophes at short interatomic distances [184, 185]. In a modified version, the damping coefficients take the following form [186, 187]: κ 3 = 1 e wr3 ab, κ 5 = 1 (1 + wr 3 ab)e wr3 ab, R ab = R ab, (4.34) (α a α b ) 1/6 where R ab is an effective interatomic distance that depends on the atomic polarizabilities and w is a damping constant, which is set to 0.39.

97 96 Chapter 4. Charge Transport in Organic olids dcv1t dcv2t dcv3t dcv4t dcv5t dcv6t dcv3t-m dcv4t-m αh n αn n αc n αs n α c /α n tr(α n ) tr(α c ) tr(αdft) n tr(αdft) c Table 4.3. Atomic polarizabilities α n for hydrogen, nitrogen, carbon, and sulfur atoms of neutral dcvnt and dcvnt-m molecules. The atomic polarizabilities α c for the charged molecules are scaled as specified. Also listed are the traces of the molecular polarizability tensors resulting from the atomic polarizabilities and corresponding traces obtained from dft calculations. All values are given in units of Å 3. Atomic Polarizabilities To calculate the electrostatic interaction energies (4.31), the first step is to determine the (isotropic) atomic polarizabilities α c,n a.theirroleis to model the (anisotropic) molecular polarizability tensor α c,n of charged or neutral molecules, which is expressed via the relay matrix A 1 by α c,n =(α µν ) with α µν = (A 1 ) µν and A µν ab ab ab = (αc,n a ) 1, for a = b, T µν ab, for a b, (4.35) as is deduced within the Applequist model of distributed polarizabilities [ ]. Based on this relationship, the atomic polarizabilities for single charged and neutral molecules are adjusted such that the molecular polarizability tensor is in accordance with its counterpart obtained from dft calculations using the b3lyp functional and the 6-311g basis set. For a series of dcvntanddcvnt-m molecules in charged and neutral states, resulting atomic polarizabilities are listed in Table 4.3. Also provided are traces of the associated molecular polarizability tensors (4.35) and of their dft equivalents. Note that the atomic polarizabilities are fitted to the chemical elements of the atoms, in this case hydrogen, nitrogen, carbon, and sulfur, and under constraining the ratios for different chemical elements to the ratios of model atomic polarizabilities from the (revised) Thole model [185]. These model atomic polarizabilities have proven successful in reproducing molecular polarizability tensors simultaneously for several representative test molecules.

98 4.2. Charge Transport imulations in Disordered Organic olids 97 Induced Atomic Dipole Moments The second step before evaluating the electrostatic interactions of a charged or neutral molecule i with its environment (4.31) is the determination of the induced atomic dipole moments p a.theprerequisitefor this are the atomic partial charges q a and polarizabilities α a of the charged or neutral molecule i and of the neutral surrounding molecules labeled according to (4.32). Then, the permanent electrostatic field arising from the atomic partial charges leads, owing to the atomic polarizabilities, to the induction of atomic dipole moments p a on the atoms a. These induced dipole moments are, however, created not only by the permanent field due to the atomic partial charges but also from the polarization field arising from the other induced atomic dipole moments [193]. Hence, the induced atomic dipole moments are determined by p µ a = α a ( T µ ab q b + T µν ab pν b ), (4.36) b j b a where the first sum is the permanent field at atom a in molecule j originating from the atomic partial charges q b on the other molecules, each creating a field µ T abq b = T abq µ b, and the second sum is the polarization field at atom a due to the other, including intramolecular, induced dipole moments p b,eachcontributingafield µ Tν abp ν b = T ab µν p ν b. To solve (4.36), a self-consistent polarization field needs to be determined, which is achieved by iterative methods [194]. At the beginning, the induced dipole moments p a are initialized by evaluating the permanent electric field due to the atomic partial charges, i.e., the first sum, at the positions of the atoms a. Then, the induced dipole moments are refined iteratively, according to the second sum: p a(0)=α µ a T µ ab q b, b j pµ a(n +1) =Ωα a T µν ab pν b (n) +(1 Ω) pµ a(n). (4.37) b a The iteration is stopped if the accuracy of the induced dipole moments is adequate: a p µ a(n +1) p µ a(n) < 10 4 D. Note that the convergence of the iterative procedure is enforced and accelerated by the technique of successive overrelaxation using a relaxation factor Ω [187, 195]. Figure 4.12 a illustrates, as an example, converged induced atomic dipole moments for a system of three dcv4t molecules, where the central molecule is positively charged. Also shown are induced molecular dipole moments, which are composed additively of the atomic dipole moments of the respective molecule. Figure 4.12 b depicts, as another example, the induced molecular dipole moments on

99 98 Chapter 4. Charge Transport in Organic olids (a) (b) x -y z Figure Induced dipole moments resulting from a positively charged molecule, illustrated in black. (a) Induced atomic dipole moments, indicated by small arrows, and resulting molecular dipole moments, by large arrows, on two surrounding neutral molecules, shown in white. (b) Induced molecular dipole moments in a crystalline dcv4t system. dcv4t molecules in a volume of radius 3.5 nm surrounding a positively charged molecule. For dcvntsystems,arelaxationfactorω = 0.3 has proven to be a reliable compromise between speed and stability of the iterative convergence procedure. ite Energy Difference According to theories of charge transfer (ection 3.2), achargetransferreactionbetweenapairofmoleculesi and j is driven by the difference in their site energies (3.19), E ij = E i E j = E int ij + E elstat ij + E ext ij. (4.38) Here, the site energy difference is decomposed into internal and electrostatic contributions, and an additional term accounting for an externally applied electric field F, E int ij = E int i E int j, E elstat ij = E elstat i E elstat j, E ext ij = q F r ij, (4.39) where q =±e is the charge and r ij = r j r i a vector connecting the molecules i and j. For typical distances between small molecules, of the order of 1 nm, and moderate fields, of up to about 10 6 V/cm, this term is always smaller than 0.1 ev. The distribution of site energy differences for all molecular pairs of the neighbor list provides an alternative definition of energetic disorder: σ = E 2 E 2. (4.40)

100 4.2. Charge Transport imulations in Disordered Organic olids 99 ince site energies are often spatially correlated, this standard deviation of site energy differences (4.40) is usually smaller than the deviation of site energies (4.6), i.e., σ σ. In fact, charge transport is only affected by neighbor list pairs of molecules, thus σ is generally a more appropriate quantity to characterize the energetic disorder. In cases where only a subset of the neighbor list contributes to charge transport, for example when the charge carrier visits only correlated low-energy regions of a disordered site energy landscape, a reduced neighbor list should be employed to calculate an effective energetic disorder σ eff σ (Chapter 6) Reorganization Energies In an energy-conserving charge transfer complex formed of two molecules i and j, charge transfer occurs at an instant of coinciding molecular energies, which is induced by the site energy modulations due to nuclear vibrations. In the high-temperature limit (ection 3.2.1), the required change of nuclear coordinates is described by the reorganization energy (3.21), λ ij = λ i + λ j = U n i U n i + U c j U c j, (4.41) where Ui n is the energy of the neutral molecule i in the nuclear geometry of its charged state, and similarly U c j is the energy of the charged molecule j in the nuclear configuration of its neutral state. Ui n and U c j are the usual energies at the nuclear equilibrium geometries, which also contribute to the internal site energy (4.29). To compute a reorganization energy (4.41), four ab initio calculations are required: two with geometry optimization, yielding the energies at the nuclear equilibrium configurations, and two with geometry constraints to the respective opposite nuclear configuration. The determination of the reorganization energies for all molecular pairs of the neighbor list is a computationally intensive task. However, for one-component systems with negligible conformational variations, the reorganization energies are identical and only a single evaluation is required. This is still a good approximation if molecular conformations vary, but energetic disorder is substantially larger than the variance of the reorganization energies. Nonetheless, from a general point of view, where potential energy surfaces are different for each molecule, the reorganization energies are properties of the charge transfer complexes: λ ij λ ii.moreover,themolecularreorganization energies λ i and λ j can vary for charging and discharging a molecule, hence the reorganization energies of the reverse processes may be different: λ ij λ ji.

101 100 Chapter 4. Charge Transport in Organic olids Charge Dynamics Having determined the list of hopping sites and for each pair i and j the charge transfer parameters as well as the rates ω ij,thenextstepistostudythedynamicsofthecharge carrier using the master equation. As already introduced, the master equation (4.8) is a set of coupled differential equations in time τ for the occupation probabilities p i that are associated with the localized electronic states ψ i. A stochastic solving approach is provided by the kinetic Monte Carlo method. 7 In this procedure, the time evolution of the charge carrier is explicitly simulated leading to a realization of the master equation, a so-called Markov chain. A Markov chain refers to a sequence of states ψ k,startingfromarandominitialstate, ψ k1 τ 1 ψ k2 τ 2 ψ k3 τ 3, (4.42) which comprises time intervals of residence τ,suchthatattimeτ aspecificstate ψ i occurs with the probability p i that satisfies the master equation. Thus, generating an ensemble of Markov chains allows the deduction of the occupation probabilities at time τ by p i = n i /n,wheren is the total number of Markov chains and n i is the number of chains in the state ψ i at time τ.thestationary solution of the master equation is obtained for large times τ,whereoccupationprobabilitiesbecometime-invariant, p i / τ = 0. In the case of ergodic systems, where Markov chains at large times τ become independent of the randomly chosen initial state, the stationary occupation probabilities can be alternatively derived from a single, sufficiently long Markov chain, which is evaluated over its time.thestationaryprobabilitiesthenfollowfromp i = l i /l, where l is the temporal length of the Markov chain and l i is the accumulated time of residence in the state ψ i. Note that in the case of isolated systems, i.e., systems without sources, sinks, or circular currents, the stationary solution corresponds to the thermodynamic equilibrium, where rates are known to obey the condition of detailed balance,thatis,p j ω ji = p i ω ij. In practice, this can be used to test whether the system is ergodic or not by correlating ln(p i ) and the site energy E i.indeed,theratiosoftheforwardandbackwardrates are determined solely by the site energy differences since ω ji /ω ij = exp( E ij /k B T), as follows from the Marcus rate (3.18) and the assumption that the reorganization energies are equal for both directions, λ ij = λ ji. 7 Alternatively, the master equation can be solved numerically [196]. Numerical algorithms may, however, become unstable for strongly varying rates, as occurs in systems with high energetic disorder.

102 4.2. Charge Transport imulations in Disordered Organic olids 101 (a) 1. Initialization (b) ω ik1 et an initially occupied site i ω ik2 et the starting time τ et the stopping time l 2. Residence Time Calculate the total escape rate to neighboring sites: ω i = Perform the transition: i j 4. Continuation ω ij j Generate a residence time interval: τ = ln r, r (0, 1] Evolve the time: τ τ + τ 3. Transition Generate target site j: biggest j for which r, ω i j k γ ω ikγ ω i r (0, 1] r ω i (c) r ω i ω ik3 ω ij ω ik1 ω ik1 +ω ik2 ω ik1 +ω ik2 +ω ik3 ω ik1 + +ω ij Continue if the stopping time is not yet reached: τ < l top otherwise ω i Figure (a) Flowchart of the kinetic Monte Carlo scheme according to the variable step size method. (b, c) Graphical representation of the target site selection (step 3). Adapted from Reference [197]. Kinetic Monte Carlo imulation The kinetic Monte Carlo method [197, 198] is tasked to simulate the charge carrier dynamics and generate desired Markov chains. Asuitablesimulationalgorithmisthevariablestepsizemethod[ ],alsoknown as the n-fold way, which is depicted as a flowchart in Figure 4.13 a. In the first step, the initiallyoccupied site i is selected, the starting time τ is set, without loss of generality τ = 0, and the stopping time l is specified, defining the length of the Markov chain. The second step is the determination of the residence time before the charge carrier is propagated. According to the Markov approximation underlying the master equation, the charge carrier is memoryless and has in each short time interval the same prob-

103 102 Chapter 4. Charge Transport in Organic olids ability of escaping from site i, which gives rise to exponential decay statistics[197]. Hence, the probability that site i is still occupied at a later time τ + τ is given by exp ( ω i τ), whereω i is the total escape rate from site i to all accessible sites from the neighbor list. The time interval of residence τ can thus be obtained by solving exp ( ω i τ)=r with r being a random number in the unit interval. The third step is the determination of the target site to which the carrier is propagated. Based on the set of all accessible sites k γ,thetargetsitej has to be selected with a probability proportional to ω ij /ω i.figure4.13billustratesthisselectionschemebyastackofobjects, each representing, by its object height, a possible escape rate ω ikγ from site i to an accessible neighboring site k γ,suchthatthetotalstackheightisω i [197]. Generating a random number r in the unit interval yields a random position r ω i along the stack, which points with the desired probability to one of the objects, say ω ij,andtherefore indicates the selected target site j. Figure4.13cillustratesthepracticalimplementation of this scheme, employing an array of successively extending partial sums, where each array element contains the accumulated height of all preceding objects up to the current object. Consequently, the target site j is obtained by the array element with the biggest index j whose partial sum is not larger than r ω i.afterpropagatingthe the charge carrier from site i to j,thesituationcorrespondstothestartingpointofthe algorithm. Thus, in the fourth step, the procedure is continued, provided the stopping time l is not yet reached Macroscopic Observables Knowing the occupation probabilities p i of the ensemble of states ψ i at the time τ, and thus the density operator ρ (4.7), one is in the position to determine ensemble averages of observable quantities a with associated operators â: a =tr (â ρ ) = i p i ψ i â ψ i = p i a i, (4.43) i where a i denotes the expectation value of a in the pure state ψ i.iftheobservablea does not explicitly depend on time, the master equation (4.8) can be invoked to express the average a by the following relation: a τ = i, j (p j ω ji p i ω ij ) a i = p i ω ij (a j a i ). (4.44) i, j

104 4.2. Charge Transport imulations in Disordered Organic olids 103 Electric Current If the position of the charge r is an observable, the time derivative of its average is the charge velocity v = r τ = p i ω ij (r j r i )= 1 2 (p i ω ij p j ω ji ) r ij, (4.45) i, j i, j where r ij = r j r i corresponds to a vector pointing from site i to j and the symmetrization of the summation on the right hand side follows from r ij = r ji.multiplication with the charge q =±e yields the total current in the system, j = q v = 1 2 c ij r ij, c ij = q (p i ω ij p j ω ji ), (4.46) i, j and the contribution from a specific pair of sites i and j is denoted as edge current c ij. Charge Carrier Mobility The zero-field mobility tensor µ of the charge carrier is associated with the diffusion tensor D by the Einstein relation, µ = 1 k B T D, (4.47) and can thus be obtained by studying particle diffusion in the absence of external fields. Using the squared particle displacement r 2 as an observable, the diffusion tensor takes the following form, with superscripts indicating the Cartesian components: 6D µν = rµ r ν τ = p j ω ji ( r µ i rν i rµ j rν j ) = p j ω ji (r µ i rν i rµ j rν j ). (4.48) i, j i, j Alternatively, the diffusion tensor can be directly determined from the ensemble of Markov chains with temporal length l,each yielding a charge displacementvector r by the vector connecting the initially and finally occupied sites, such that 6D µν l = r µ r ν. (4.49) This method has the advantage that it can be immediately extended to systems with periodic boundary conditions. In this case, the charge displacement vector is obtained by unwrapping the diffusion trajectory defined by the Markov chain, as is exemplified in Figure 4.14 for a crystalline system of 2880 dcv4t molecules.

105 104 Chapter 4. Charge Transport in Organic olids x y r z Figure Unwrapped charge diffusion trajectory in a crystalline system of 2880 dcv4t molecules with periodic boundaries, indicated as a cuboid. The charge displacement vector r connects the initially and finally occupied sites, marked in violet and red, respectively. The Markov chain has a temporal length l = 10 3 s and the axes are given in units of nm. The charge carrier mobility tensor µ under an externally applied electric field F can be calculated by the following relation using the average charge velocity (4.45): v =µf. (4.50) Alternatively, the field-dependent mobility tensor can be directly determined from a Markov chain. To this end, the charge velocity is calculated from the charge displacement vector between the initially and finally occupied sites divided by the time length l. Projecting this velocity on the direction of the field yields the field-dependent charge carrier mobility in this particular direction. To ensure adequate statistics, the mobility can be averaged over an ensemble of Markov chains. Again, this method is particularly suitable if periodic boundary conditions are employed, where the charge displacement vector follows from the unwrapped charge carrier trajectory.

106 Part III. Results

107

108 Chapter 5. Charge Transport imulations in Organic Crystals Within the field of organic photovoltaic cells based on small molecules (Chapter 1), a particular successful class of donor compounds are dicyanovinyl-substituted oligothiophenes (dcvnt). When preparing pure dcvntmaterialsamples,as usedinplanar heterojunction cells, the substances typically adopt crystalline phases. Here, a comparative charge transport simulation study is performed for the crystalline phases of the terthiophene (dcv3t), a methylated derivative (dcv3t-m), the quaterthiophene (dcv4t), and a methylated derivative (dcv4t-m), all of which are shown in Figure 5.1. First, attention is directed to single crystals of all four compounds (ection 5.1). Comparing the charge transport behavior leads to the conclusion that crystal structures characterized by a well-defined π-stacking are disadvantageous for an efficient transport. The microscopic origins of this counterintuitive finding are elucidated and provide an explanation for a similar trend observed in experimental measurements. econd, the study is extended to thin films of the two quaterthiophenes (ection 5.2). uch thin film layers can exhibit molecular packings different than the single crystals and thus altered charge transport capabilities. Taking into account the thin film molecular packings yields charge carrier mobilities with a systematically improved agreement with experimental device measurements. NC NC dcv3t NC dcv3t-m Me Me NC dcv4t NC NC dcv4t-m NC Me Me Me Me NC Figure 5.1. Chemical structures of terminally dicyanovinyl-substituted terthiophene (dcv3t), a methylated derivative compound (dcv3t-m), dicyanovinyl-substituted quaterthiophene (dcv4t), and a methylated derivative compound (dcv4t-m).

109 108 Chapter 5. Charge Transport imulations in Organic Crystals 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals Within the collaborating group of Bäuerle at the Institute of Organic Chemistry II and Advanced Materials in Ulm, Germany, the dicyanovinyl-substituted terthiophenes, dcv3t and dcv3t-m, and quaterthiophenes, dcv4t and dcv4t-m, have been purified via vacuum sublimation at temperatures well below their melting points [202, 203]. In these processes, highly pure molecular crystals of up to macroscopic size (linear dimensions of roughly 0.1 mm) could be grown, which is sufficiently large for an investigation by single crystal x-ray crystallography. In this technique, a crystal sample is mounted in a goniometer for gradual rotation around the three axes, while being exposed to an incident x-ray beam. From the scattered radiation one obtains a diffraction pattern with angular dependence on the crystal orientation, which allows one to infer the periodic structure of the crystal on an atomic length scale. This crystal structure includes both the unit cell, i.e., the parallelepiped defining the crystal lattice, and the repeating motif, i.e., the precise molecular arrangement at each lattice point. uch crystallographic information is publicly available for dcv4t and dcv4t-m in Reference [202] and [203], respectively, and is expected to be published in the near future for dcv3t and dcv3t-m. In what follows, the four crystal structures are taken as a basis to perform microscopic charge transport simulations (ection 4.2). First, mesoscopic morphologies at room temperature are generated and the inherent disorder is quantified (ection 5.1.1). ubsequently, the charge transfer parameters between the neighboring molecules are evaluated (ection 5.1.2). And finally, charge transport is investigated by calculating the charge carrier mobility tensors and relating the results to the energetic disorder and the charge transporting networks built from the transfer integrals (ection 5.1.3). The reported scientific results are the subject of prior publication [1], and are presented below in more detail using partially similar terms and illustrations Morphological Disorder As starting configurations for molecular dynamics simulations, supercells with 2880 molecules in crystalline arrangement (and periodic boundary conditions) are created by the multiplication of the crystal repeating motifs. To this end, the motifs are translated multiple times along the base vectors a, b, andc of the (primitive) unit cells: fordcv3t,24 5 6for dcv4t, fordcv3t-m, for dcv4t-m. Note that additional orthogonal coordinate systems are introduced which are defined

110 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals 109 a 0 c b x -y z Figure 5.2. Morphology of the dcv4t molecular crystal at room temperature, truncated in the y-direction for clarity. In the left third, the molecules are drawn explicitly, in the middle third, only their molecular backbones are shown, and in the right third, only their centers of mass are depicted. The parallelepiped indicates the (primitive) unit cell of the dcv4t crystal, which is nearly orthogonal and contains a repeating motif of four molecules. by the x-axis along the a-vector,the y-axis within the ab-plane, and the z-axis orthogonal to the prior two axes. Then, the prepared supercells are subjected to molecular dynamics simulations (ection 4.2.2) at a pressure of 1 bar, a temperature of 300 K and for a duration of 10 ns. This yields realistically disordered morphologies, as is illustrated for dcv4t in Figure 5.2. The view shows a slice of the y-direction and indicates molecular backbones and centers of mass by gray tubes and black spheres, respectively. In addition, the dcv4t unit cell, which contains four molecules, is shown. For all compounds, an overview of unit cells is provided in Figure 5.12 a on Page 122/123. Molecular Conformations The molecular conformations within the room-temperature morphologies are analyzed by evaluating the dcv-thiophene (cc c ca s) and thiophene-thiophene dihedral angles (s ca ca s) as defined in the force field (ection 4.2.1). The ensemble distributions of all dihedral angles in one molecular dynamics snapshot are shown in Figure 5.3. As can be seen, the molecules are almost flat on average,whichisinagreementwiththeexperimentallyobtainedcrystalstructures.

111 110 Chapter 5. Charge Transport imulations in Organic Crystals Probability Probability Prob. 0 Prob N Dip. Mom. (d µ /D) N N N dcv3t d x d y d z 90 N Dihedral Angle (θ/ ) N dcv3t-m N N N Dip. Mom. (d µ /D) Prob. 0 Prob N Dip. Mom. (d µ /D) N N dcv4t 90 dcv4t-m N Dip. Mom. (d µ /D) Dihedral Angle (θ/ ) N N N 180 Figure 5.3. Distributions of the dcv-thiophene and thiophene-thiophene dihedral angles in molecular crystals at 300 K. The dcv-thiophene dihedrals are shown in red and yellow, while thiophene-thiophene dihedrals are shown in blue, green, and cyan. The insets illustrate distributions of the molecular dipole moment components. Adapted from Reference [1]. This observation is expected for the bare oligomers, dcv3t and dcv4t, since ab initio calculations with geometry optimization performed on single molecules in vacuum result in planar conformations. The methylated compounds, dcv3t-m and dcv4t-m, however, have significantly twisted vacuum-optimized geometries and, consequently, the conformational planarization in the crystal arises due to non-bonded interactions with the neighboring molecules. As will be discussed below, this planarization effect has an important consequence for the reorganization energies. The deviations from the planar average geometries of the molecules, as a consequence of thermal fluctuations, result in conformational disorder. This conformational disorder can be quantified by the widths of the dihedral angle distributions. As seen in Figure 5.3, the widths are of the order of 30 for the dcv-thiophene and 20 for the thiophene-thiophene dihedral angles.

112 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals 111 x -y z Figure 5.4. Alternately oriented molecular dipole moments in a dcv3t molecular crystal, shown as arrows at the molecular centers of mass (black spheres). The dipole moments are color-coded according to their direction, that is, opposite directions have complementary colors. The hue represents the azimuth angle, while lightness reflects the elevation angle (hsv color model). Adapted from Reference [1]. Molecular Dipole Moments Based on the atomic partial charges of the dcvnt molecules, as defined in the force field (ection 4.2.1), permanent molecular dipole moments can be evaluated for the molecules in the room-temperature morphologies. The insets of Figure 5.3 show the ensemble distributions of molecular dipole moments, resolved into their components (in the orthogonal xyz-coordinates). For the terthiophenes, dcv3t and dcv3t-m, one observes bimodal distributions, which indicates that the dipole moments are non-zero on average. Indeed, the planar average conformations of the terthiophenes have an asymmetric or mirror-symmetric alignment of the electron-withdrawing dcv substituents, which leads to non-vanishing dipole moments. In contrast, the quaterthiophenes, dcv4t and dcv4t-m, are centro-symmetric and thus have zero dipole moments on average. This observation is in accordance with the experimental crystal structures: the unit cells of the terthiophenes contain an even number of molecules with alternating orientation, such that the dipole moments are mutually compensating. The deviations of the molecular dipole moments from their average values and orientations, which result from the conformational disorder, give rise to dipolar disorder. Figure 5.4 illustrates for the bare terthiophene, dcv3t, the fluctuating alternately oriented dipole moments in a small section of a molecular dynamics snapshot. Figure 5.5 shows similarly for the bare quaterthiophene, dcv4t, the fluctuating randomly oriented dipole moments. As will be discussed below, the dipolar disorder is closely related to disorder in the site energies.

113 112 Chapter 5. Charge Transport imulations in Organic Crystals x -y z Figure 5.5. Randomly oriented molecular dipole moments in a dcv4t molecular crystal, shown as arrows at the molecular centers of mass (black spheres). The dipole moments are color-coded according to their direction, that is, opposite directions have complementary colors. The hue represents the azimuth angle, while lightness reflects the elevation angle (hsv color model). Adapted from Reference [1] Charge Transfer Parameters In order to construct the charge transport Hamiltonian, the morphologies are partitioned on conjugated segments and rigid fragments (ection 4.2.3). Accordingly, the π-conjugated dcvnt moleculesarerepresentedbyseparaterigidfragmentsforeach dcv or (methylated) thiophene unit, and by single conjugated segments, whose centers of mass correspond to the hopping sites i. Theneighborlistofadjacenthopping sites i and j is created from all sufficiently close molecular pairs, characterized in that any pair of associated mutual rigid fragments is below a cutoff of 0.7 nm. Reorganization Energies The computation of reorganization energies is based on ab initio calculations for single charged and neutral molecules, once with and once without geometry optimization (ection 4.2.6). In the present case of one-component systems, it is sufficient to compute for each system a universal reorganization energy λ. To account for the conformational planarization of the molecules in the bulk, the calculations with geometry optimization are carried out while the dihedral angles are constrained to their average values in the molecular crystals. As a result, similar reorganization energies λ for all four compounds are obtained, which are listed in Table 5.1. The calculations are conducted by dft, the b3lyp functional, and the 6-311g basis set.

114 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals 113 dcv3t dcv4t dcv3t-m dcv4t-m λ σ Table 5.1. Reorganization energy λ and energetic disorder σ for pairs from the neighbor list, both given in units of ev. Note that unconstrained geometry optimization calculations, as are usually carried out, lead to an increase of the reorganization energies of the methylated compounds, dcv3t-m and dcv4t-m, by 0.15 ev and 0.13 ev, respectively. This is a consequence of their twisted neutral, but planar charged geometry, contrary to the bare compounds, which are flat in either charge state. These higher values of the reorganization energies would lead to an order of magnitude decrease in the charge carrier mobilities, which are determined later on. ite Energies The observed dipolar disorder is an indication of disorder in site energies: for instance, in phenomenological Gaussian disorder models (ection 4.1.2), randomly oriented dipoles are known to cause correlated energetic disorder [ ]. To quantify the energetic disorder in the molecular crystals, the site energies E i of all molecules i are determined based on the electrostatic interactions between atomic partial charges and induced dipole moments according to the Thole model (ection 4.2.5). As an example, Figure 5.6 visualizes the resulting site energies in a dcv4t morphology by means of color-coded hopping sites, i.e., molecular centers of mass. The range from low to high site energies is indicated by a color range from blue to red. Obviously, site energies are spatially correlated, but energetic defects also exist, as can be seen by singular blue or red hopping sites. For all compounds, the distributions of the site energy differences, E ij = E i E j, for pairs from the neighbor list are shown in Figure 5.7. In addition, the insets provide the spatial and temporal correlation functions, C E = (E i E )(E j E ) (E i E ) 2, (5.1) which are one if site energies are fully correlated and zero if they are uncorrelated. The widths of Gaussian functions fitted to the distributions yield the neighbor-listbased energetic disorder σ, as defined in (4.40). The results, summarized in Table 5.1, indicate that the four systems are characterized by fairly similar energetic disorder, which is substantial and approximately comparable to many amorphous systems [4].

115 114 Chapter 5. Charge Transport imulations in Organic Crystals x -y z Figure 5.6. Hopping sites (molecular centers of mass) in a dcv4t molecular crystal, colored according to their site energy. The range from low to high site energies is indicated by a color range from blue to red. Evidently, the site energies are spatially correlated and occasionally exhibit energetic defects. Adapted from Reference [1]. uch significant energetic disorder is unexpected for crystalline systems, but is a consequence of the molecular architecture of the dcv-substituted oligothiophenes, having the electron-withdrawing dcv groups attached in non-axial orientation and the dcv and thiophene units connected by single bonds. uch an acceptor-donor-acceptor architecture combined with fluctuations of the dihedral angles leads to substantial variations of the local electric fields. As mentioned before, the terthiophenes, dcv3t and dcv3t-m, have equilibrium conformations with non-zero permanent dipole moments, which are compensated in a unit cell. If the energetic disorder in the terthiophene molecular crystals is calculated based on the atomic partial charges only, thus neglecting effects of molecular polarization, the width of the site energy distribution becomes 0.17 ev and 0.30 ev, respectively. These significantly higher values reveal that not only the compensation of dipole moments in a unit cell [204], but also the large polarizabilities of the molecules (ection 4.2.5) play an important role in reducing the energetic disorder. For the quaterthiophenes, dcv4t and dcv4t-m, which have zero dipole moments in their equilibrium conformations, the neglect of polarization effects has no significant influence on the energetic disorder.

116 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals 115 Probability Probability r ij /nm t/ps r ij /nm r ij /nm dcv3t 1 1 CE(rij) CE(rij) dcv3t-m CE(t) CE(t) r ij /nm ite Energy Difference ( E ij /ev) ite Energy Difference ( E ij /ev) t/ps CE(rij) CE(rij) dcv4t 0 dcv4t-m CE(t) CE(t) t/ps t/ps Figure 5.7. Distributions of the site energy differences for pairs from the neighbor list. The width of Gaussian fit functions yields values of energetic disorder σ, listed in Table 5.1. The insets show spatial correlation functions and time correlation functions. The latter are calculated from molecular dynamics simulations of 0.5 ps and bare Coulomb interactions. Adapted from Reference [1]. Transfer Integrals The transfer integrals J ij between pairs of molecules i and j from the neighbor list are determined by means of the zindo method (ection 4.2.4). The set of transfer integrals constitutes a percolating network, which provides information about the directionality and dimensionality of charge transport. In Figure 5.8, this charge transporting network is illustrated for a dcv4t snapshot. To highlight its topological connectivity pattern, an enlarged view of a ten unit cell volume is also shown. Here, hopping sites are depicted by black spheres and transfer integrals are represented by intermediate colored bonds; their thickness displays the electronic coupling strength, while their color indicates the crystallographic direction. The discrimination of the distinct crystallographic directions is accomplished by shifting all bond

117 116 Chapter 5. Charge Transport imulations in Organic Crystals x -y z Figure 5.8. Charge transporting network of a dcv4t molecular crystal, truncated in the y-direction for clarity. Hopping sites (molecular centers of mass) are depicted as black spheres, while transfer integrals between molecular pairs are illustrated as colored bonds. The bond thickness reflects the electronic coupling strength, while the bond color indicates the crystallographic direction. The enlarged view corresponds to a volume of ten unit cells.

118 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals 117 x -y z Figure 5.9. Vectors between all hopping site pairs from the neighbor list, shifted to the origin (at the center). Groups are defined by the smallest cuboids enclosing clusters of vector endpoints. The groups are identified with distinct crystallographic directions, indicated by different cuboid colors. If directions are equivalent by the crystal symmetry, groups are united. Adapted from Reference [1]. vectors to the origin and finding groups of vector endpoints by an agglomerative hierarchical cluster algorithm [205]. This scheme is illustrated for the dcv4t system in Figure 5.9: all bond vectors of the system are drawn simultaneously pointing from the center to the exterior. The vector endpoints form clusters, which are enclosed by cuboids, whose coloring corresponds to the respective crystallographic directions. For all molecular crystals, an overview of charge transporting networks is shown in Figure 5.12 b on Page 122/123. For a quantitative analysis, Figure 5.10 depicts the associated total and direction-resolved transfer integral distributions, drawn in black and the above-defined colors, respectively. The total distributions comprise several peaks, which are clearly attributable to the underlying directional distributions, characterized by single pronounced peaks (the only exception is the direction in dcv3t-m marked in red, along which electronic coupling alternates). In addition, Table 5.2 lists the average directional transfer integrals as well as the average hopping site separations. One recognizes that the molecular crystals of the bare compounds, dcv3t and dcv4t, have a well-defined π-stacking direction (red) with an average coupling of 10 3 ev 2 and center-of-mass separation of 0.4 nm, while other crystallographic directions (blue and green) have several orders of magnitude lower coupling strengths. This demonstrates that the presence of the strongly coupled one-dimensional π-stacking direction is at the expense of the coupling strength of other directions. Moreover, these other directions are oriented almost perpendicular to the π-stacking, as seen in Figure 5.12 b on Page 122/123. In contrast, analyzing the methylated compounds, dcv3t-m and dcv4t-m, reveals the absence of a comparably pronounced π-stacking direction; instead, their charge transporting networks comprise multiple crystallographic directions of moderate coupling strengths.

119 118 Chapter 5. Charge Transport imulations in Organic Crystals 75 dcv3t 50 dcv4t Probability Probability dcv3t-m dcv4t-m Transfer Integral (J ij /ev) Transfer Integral (J ij /ev) Figure Total and direction-resolved distributions of the transfer integrals, drawn in black and distinct colors, respectively. The coloring of directional transfer integrals corresponds to the topological connectivity patterns in Figure 5.12 b on Page 122/123. Average values of the distributions are listed in Table 5.2. Adapted from Reference [1]. dcv3t dcv4t dcv3t-m dcv4t-m J 2 r J 2 r J 2 r J 2 r Table 5.2. Average direction-resolved transfer integrals J 2, given in ev 2, and associated average hopping site separations r, given in nm. The coloring of crystallographic directions corresponds to the topological connectivity patterns in Figure 5.12 b on Page 122/123.

120 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals 119 dcv3t dcv4t dcv3t-m dcv4t-m µ 1 σ= µ 2 σ= µ 3 σ= µ µ µ µ exp Table 5.3. Eigenvalues of the zerofield mobility tensors µ σ=0 and µ, calculated in the absence and presence of energetic disorder, respectively. Experimentally measured organic field-effect transistor mobilities µ exp are available for dcv4t and dcv4t-m [203]. All values are given in units of cm 2 /Vs Charge Carrier Mobility In crystalline phases, charge transport generally exhibits anisotropic behavior, i.e., the charge carrier mobility is dependent on the crystallographic direction. In order to determine the main transport directions and link them to the packing motifs, charge diffusion without an externally applied electric field is studied first. Based on the kinetic Monte Carlo method (ection 4.2.7), the diffusion tensor (ection 4.2.8) of each molecular crystal is determined by averaging over 2000 charge displacement vectors for time intervals of 10 3 s. Then, the zero-field mobility tensor follows from the Einstein relation (4.47). In order to assess the effects of the energetic disorder on charge transport, the zerofield mobility tensor is determined for each crystalline system twice. On the one hand it is calculated for a simplified model system where energetic disorder is turned off, i.e., the site energies are set to zero, and on the other hand it is calculated for the full model system including the site energies. These two zero-field tensors are denoted as µ σ=0 and µ,respectively,andtheireigenvaluesareprovidedintable5.3.thebounding eigenvalues, i.e., the first and the third, specify the highest and lowest mobilities and define the mobility ranges which are displayed in the left panels of Figure The striped and filled areas correspond to the mobility ranges in the absence and presence of energetic disorder, respectively. To link the topological connectivity of the systems to the anisotropy of the charge carrier mobility, the mobility tensors of the four molecular crystals are depicted as ellipsoids in Figure 5.12 c on Page 122/123; the colored arrows originating from the center of the ellipsoids indicate the crystallographic directions which were identified previously (Figure 5.12 b). Each ellipsoid is defined by its three principle axes pointing along the eigenvectors of the associated mobility tensor and by the associated equatorial radii given by the square root of the tensor eigenvalues (Table 5.3). Hence, the

121 120 Chapter 5. Charge Transport imulations in Organic Crystals 10 1 dcv3t 10 1 dcv4t cm2 ) Mobility (µ/ ) Vs dcv3t-m dcv4t-m Mobility (µ/ cm2 Vs Electric Field (F/ V cm )1/2 Electric Field (F/ V cm )1/2 Figure Left panels: ranges of the zero-field hole mobilities with and without energetic disorder, shown by filled and striped areas, respectively. The ranges are bounded by the largest and smallest eigenvalues of the mobility tensors. For dcv4t and dcv4t-m, experimentally measured organic field effect transistor mobilities are indicated by stars [203]. Right panels: hole mobilities as a function of an electric field applied along different directions. Directions of the largest and smallest principal axes of the zero-field tensors are indicated by dashed, and other crystallographic directions by solid lines, colored according to Figure 5.12 b, c on Page 122/123. Adapted from Reference [1]. transport direction with the maximum mobility, i.e., the first eigenvalue, is given by the longest ellipsoidal principle axis. The ellipsoids drawn in red correspond to the disregard of site energy variations, while those in gray include the effects of energetic disorder. If red ellipsoids are omitted, their orientation is similar to gray ones. First, the results obtained in the absence of energetic disorder (striped areas in Figure 5.11 and red ellipsoids in Figure 5.12 c) are analyzed and the different crystal structures compared. As regards the main transport directions, i.e.,thelongestellipsoidal axes, the following observation can be made. The crystals of the bare compounds, dcv3t and dcv4t, exhibit their maximum mobility along the π-stacking direction where the electronic coupling is the strongest (red). The crystals of the methylated compounds, dcv3t-m and dcv4t-m, where π-stacking is not well defined, show the

122 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals 121 highest mobility along a superposition direction of multiple crystallographic directions with moderate couplings (red and blue). When comparing the mobility values along the main transport direction, a rather counterintuitive result is observed. In spite of weaker electronic couplings, the crystals of the methylated dcv3t-m and dcv4t-m have an order of magnitude higher (maximum) mobilities as compared to the bare dcv3t and dcv4t. ince energetic disorder is switched off, this effect can only be attributed to a favorable connectivity pattern of the dcv3t-m and dcv4t-m charge transporting networks. In other words, the pronounced π-stacking with strong electronic couplings, found for the bare dcv3t and dcv4t, is disadvantageous for efficient charge transport. When analyzing the transfer integrals, it was discovered that this onedimensional π-stacking of the bare compounds leads to poorly coupled perpendicular transport directions, and inhibits electronic coupling in other directions completely (Figure 5.12 b). Thus, the topological connectivity of the methylated compounds with weaker couplings but better interconnection is favorable. 1 econd, charge transport is analyzed in the presence of energetic disorder (filled areas in Figure 5.11 and gray ellipsoids in Figure 5.12 c). As expected, the incorporation of energetic disorder leads to a reduction of the mobility values for all crystals. However, the mobility of the bare dcv3t and dcv4t crystals decreases by two orders of magnitude, while for the methylated dcv3t-m and dcv4t-m systems, it is only reduced by one order of magnitude. This discrepancy cannot be completely attributed to the slightly smaller energetic disorder present in the methylated systems, but mainly results from the aforementioned topological differences of their charge transporting networks. Indeed, considering the longest ellipsoidal principal axis, one realizes that the main transport directions of the bare dcv3t and dcv4t crystals are no longer aligned with the π-stacks. Instead, the maximum mobility is now along a superposition direction of the strongly coupled π-stacking (red) and other crystallographic directions with poor coupling (green). This reorientation of the mobility tensors is a consequence of energetic defects in the one-dimensional π-stacks, as seen in Figure 5.6. To bypass defective sites, a charge carrier has to escape to an adjacent π-stack by hopping along one of the perpendicular side directions with poor coupling. To illustrate this behavior, an external electric field of 800 (V/cm) 1/2 is applied along the main 1 Another aspect of the favorable topological connectivity for the methylated dcv3t-m and dcv4t-m crystals are the larger average separations of hopping sites. Often, the transfer integrals (and hence the rates) decay exponentially with the separation, thus, larger site-site separations lead to lower mobilities (which are proportional to the rate times the separation). However, molecules with extended π-conjugation can have relative lateral shifts which barely change the transfer integral, but increase the center-of-mass (site-site) separation, resulting in longer charge hops and higher mobilities.

123 122 Chapter 5. Charge Transport imulations in Organic Crystals dcv3t NC NC dcv4t NC NC (a) a 0 a 0 c b c (b) b x -y x -y (c) x z x z -y -y (d) z z x x -y -y z z Figure (a) Unit cells and crystal repeating motifs. (b) Charge transporting networks based on the transfer integrals. The transfer integrals between the hopping sites are shown as bonds connecting black spheres. The bond color indicates the crystallographic direction, while the bond thickness reflects the electronic coupling strength. (c) Zero-field mobility tensors. Without energetic disorder, tensors are depicted as red ellipsoids, with energetic disorder as gray ellipsoids. For dcv3t-m and dcv4t-m, red

124 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals 123 dcv3t-m NC Me Me NC dcv4t-m NC Me Me Me Me NC (a) 0 a b b a (b) -z c y 0 c -y x z x (c) -z y -y x z x (d) -z y -y x z x and gray ellipsoids have similar orientations, hence only the gray ones are shown. The colored arrows indicate the crystallographic directions of the charge transporting network. (d) Edge currents under an electric field of 800 (V/cm) 1/2 applied in the direction of the longest ellipsoidal principle axis, i.e., the main transport direction. The arrowheads indicate the direction of the current, while the thickness and color of the arrowshafts reflect the current amplitude. Adapted from Reference [1].

125 124 Chapter 5. Charge Transport imulations in Organic Crystals transport directions, and the edge currents betweenneighboring hopping sites are calculated (ection 4.2.8). The flow of electric current is visualized in Figure 5.12 d on Page 122/123, where the arrowheads indicate the direction of the current, while the thickness and color of the arrowshafts are proportional to the logarithm of the current amplitude. As expected, the average charge carrier in the bare dcv3t and dcv4t systems follows the one-dimensional π-stacking direction until a defect is reached, enforcing an escape hop to a neighboring π-stack. In contrast, the methylated dcv3t-m and dcv4t-m crystals show pathways consisting of hops along multiple crystallographic directions of moderate couplings. These multi-dimensional composite pathways enable the charge carrier to easily avoid energetic defects. To summarize, substantial energetic disorder combined with a strong π-stacking is detrimental to an efficient transport. The methylation preventing the molecules from a strong π-stacking leads to a beneficial topology of the charge transporting network with slightly lower couplings but a better interconnection of the sites. Another typical implication of spatially correlated energetic disorder is a non-linear dependence of the mobility µ on an externally applied electric field F. For example,in phenomenological Gaussian disorder models (ection 4.1.2), the Poole-Frenkel relationship ln µ F 1/2 (4.9) is reproduced [ ]. To validate this effect in the present molecular crystals, charge dynamics is studied for a charge carrier drift-diffusing under the influence of external electric fields. To this end, the kinetic Monte Carlo method (ection 4.2.7) is employed while applying field magnitudes in steps of 100 from 200 to 1000 (V/cm) 1/2 along several directions. The directional mobility for a given system with specified field vector is determined from the projection of the charge carrier velocity (ection 4.2.8) and is averaged over three independent time intervals of 10 3 s for both the forward and backward directions. The field-dependent charge carrier mobility is shown in the Poole-Frenkel plots, provided in the right panels of Figure The solid lines correspond to field directions along the distinct crystallographic directions (indicated by their colors), while the dotted lines are the directions of the two extremal axes of the zero-field mobility tensor. For all systems, the onset of a Poole-Frenkel dependence can be observed at moderate fields of approximately 400 (V/cm) 1/2, which is typical for correlated energetic disorder. Finally, the calculated mobilities are compared with experimental mobilities, which are obtained from measurements performed on organic field-effect transistors (ofet). Experimental ofet mobilities are available for the bare and methylated quaterthiophenes, dcv4t and dcv4t-m [203]. They are listed in Table 5.3 and indicated by stars

126 5.1. Dicyanovinyl-ubstituted Oligothiophenes: ingle Crystals 125 in the left panels of Figure Qualitatively, the experimental measurements support the simulation results in yielding a higher hole mobility for the methylated compound. A direct quantitative comparison cannot be performed for two reasons. First, the ofet devices comprise thin film layers with polycrystalline rather than single crystalline order [203]. Hence, grain boundaries between adjoining crystallites impair charge transport, which is a reason for the systematically lower measured mobilities as compared to the calculated ones. econd, and more importantly, crystallites in thin films can have different molecular packing than the single crystal structure and thus inherently altered charge transport capabilities. To address this issue, thin films need to be crystallographically analyzed and the implications on charge transport studied (ection 5.2).

127 126 Chapter 5. Charge Transport imulations in Organic Crystals 5.2. Dicyanovinyl-ubstituted Quaterthiophenes: Thin Films The fabrication of organic photovoltaic cells based on dicyanovinyl-substituted oligothiophenes (dcvnt) involves several processing steps of the dcvnt rawsubstances. In an initial step, raw substances are purified by vacuum sublimation leading to crystalline powder materials. Previously, attention was focused on the single crystals identified within these powders, i.e., the few crystallites coming close to millimeter size (ection 5.1). Active layers of photovoltaic cells are, however, based on thin films of nanometer thickness, which are prepared by further processing the powders and depositing the materials on specific substrates. Hence, in situ, one might find molecular ordering different from the single crystal structure, and therefore altered charge transport capabilities. Based on the dicyanovinyl-substituted quaterthiophene (dcv4t) and its methylated derivative (dcv4t-m), thin film active layers have been produced and crystallographically analyzed [3, 203]. While this analysis indicates that thin films of methylated dcv4t-m possess the same crystalline arrangement as its single crystals, it is indeed found that thin films of the bare dcv4t exhibit a packing motif different than its single crystal structure. This new dcv4t thin film crystal structure can be reconstructed based on x-ray diffraction measurements, as is addressed initially (ection 5.2.1). ubsequently, the previous single crystal study is complemented by a comparative analysis of charge transport in thin films and its impact on solar cell efficiencies (ection 5.2.2). The reported scientific results are the subject of prior publication [3], and are presented below in more detail using partially similar terms and illustrations Crystal tructure Analysis The following crystallographic investigations were performed within the collaborating group of Leo at the Institute for Applied Photo Physics in Dresden, Germany. The studies involve two experimental methods of x-ray diffraction (xrd), which, contrary to the single crystal x-ray crystallography (ection 5.1), can only provide limited information about the molecular packing. The first method is Bragg-Brentano xrd, which is applied to crystalline powders. imilar to the Debye-cherrer method, an incident x-ray beam is directed on a powder sample containing randomly oriented crystallites. Hence, the Bragg reflections for varying lattice orientations are simultaneously observed in the resulting diffraction pattern. The second method is grazing incidence

128 5.2. Dicyanovinyl-ubstituted Quaterthiophenes: Thin Films 127 Intensity dcv4t thin film dcv4t powder dcv4t single crystal cattering Angle / Intensity Figure X-ray diffraction patterns of the dcv4t thin film, the dcv4t powder, and the dcv4t single crystal. If the two dominant reflections of the dcv4t powder pattern are Gaussian-broadened and superposed, the red dotted curve is obtained. For clarity the ordinate axis has a different scale for low and high scattering angles. The inset has an adjusted ordinate axis to visualize small peaks of the dcv4t powder pattern. Adapted from Reference [3]. x-ray diffraction (gixrd), which is applied to thin films. Here,an incident x-ray beam strikes nearly parallel onto the surface of the thin film in order to increase the surface sensitivity and decrease the substrate sensitivity. For not too large scattering angles, only lattice spacings which are almost perpendicular to the surface (out-of-plane) contribute to the resulting diffraction pattern. This means in particular that patterns comprising the information on all lattice orientations are only obtained from polycrystalline thin films of randomly oriented crystallites. Given this brief methodology, an outline of the crystal structure analysis of the dcv4t thin film is provided below. X-Ray Diffraction Measurements The analysis of dcv4t thin films by gixrd measurements results in the diffraction pattern displayed in Figure 5.13 (black curve). The pattern shows a strong Bragg reflection at a scattering angle of 12.8 and several weaker ones at 26.0 and While the absence of diffuse scattering indicates a high crystallinity of the thin film, the small number of Bragg reflections reveals that crystallites have a preferred orientation and no random alignment. The broadness of the reflections is a consequence of the small layer thickness of 50 nm and the resulting small crystallite sizes. 2 Notice that the weakness of the reflection at 26.0 is an artifact of the gixrd method. This is because the corresponding lattice spacing is parallel to the thin film surface (in-plane) and does therefore not contribute to the diffraction pattern. In fact, two-dimensional giwaxs measurements on thin films of dcv4t and C 60 show not only a strong Bragg reflection at 12.8 in the out-of-plane direction but also a strong Bragg reflection at 26.0 in the in-plane direction [206].

129 128 Chapter 5. Charge Transport imulations in Organic Crystals In addition to the dcv4t thin film, the dcv4t powder is investigated using Bragg- Brentano xrd resulting in the diffraction pattern depicted in Figure 5.13 (red curve). The pattern shows not only more Bragg reflections, as expected, but also much sharper ones, which is due to larger crystallite sizes. Apparently, the intense reflections of the dcv4t powder are at similar scattering angles as the visible reflections of the dcv4t thin film. Moreover, if the two dominant powder reflections at 12.7 and 13.3 are Gaussian-broadened and superposed (red dotted curve), the dominant thin film reflection can be reproduced. From these findings, it is concluded that the dcv4t powder has the same crystal structure as the dcv4t thin film. Lastly, the diffraction patterns of both the dcv4t thin film and the dcv4t powder are compared to the pattern of the dcv4t single crystal studied previously (ection 5.1). This single crystal pattern, which is obtained from a powder xrd simulation, is visualized in Figure 5.13 (gray curve). It is obviously very different: it does not match the dominant thin film reflection (or the superposition of the two powder reflections), but instead exhibits four prominent Bragg reflections for scattering angles below Upon closer inspection, three of these reflections are found to have a very low intensity in the powder pattern, as is seen in the inset of Figure 5.13, which has an adjusted scale of the ordinate axis. This suggests that the dcv4t powder actually consists of two crystallographic phases, with the dcv4t single crystal phase being the minor phase. In fact, this minor phase and the major dcv4t thin film phase have a weight content of 2wt.%and98wt.%,respectively[3]. Crystal tructure Reconstruction Now it is briefly outlined how the crystal structure of the dcv4t thin film is reconstructed from the measured diffraction data. To this end, one should remember that a diffraction pattern is given in reciprocal space and reflects the long-range periodicity of molecular packing (by Bragg reflections). Also, by performing a Fourier transformation into real space, one obtains an atomic pair distribution function (pdf), which directly represents the short-range order of the molecular arrangement. For crystal structure analysis on an atomic scale, the data representation as a pdf may be considered the more appropriate representation [207]. Here, a pdf-based algorithm for crystal structure reconstruction is employed which was developed at the Institute for Applied Photo Physics in Dresden, Germany [208]. Without going into details, the algorithm aims at reproducing a reference pdf by a best-fitting model pdf, which is obtained by iteratively refining atomic model crystal structures. The model crystal structures are generated based on predefined triclinic

130 5.2. Dicyanovinyl-ubstituted Quaterthiophenes: Thin Films 129 Pair Distribution Function dcv4t powder dcv4t model Distance / nm Figure Atomic pair distribution functions (pdf) of the dcv4t powder and the dcv4t model crystal structure. The dcv4t powder pdf is obtained from the corresponding xrd pattern in Figure 5.13 by a Fourier transformation. The dcv4t model pdf is calculated from a model crystal structure. Adapted from Reference [3]. lattices combined with a repeating motif of a single molecule with rigid molecular geometry. After convergence of the iteration, the final unit cell parameters are additionally refined using the Rietveld analysis [209]. In order to reconstruct the dcv4t thin film crystal structure, the dcv4t powder diffraction pattern, which has a significantly higher resolution, is Fourier transformed into a reference pdf. This dcv4t powder pdf, shown in Figure 5.14 (red curve), is, after applying the reconstruction algorithm, well reproduced by the best-fitting dcv4t model pdf (black curve). Note that there are several reasons why the two functions do not perfectly agree [3]. These include in particular that the dcv4t powder is actually a mixture of two phases and that the algorithm treats the repeating motif as a rigid body. The reconstructed dcv4t thin film crystal structure is depicted in Figure 5.16 a c on Page 132/133, together with the previously studied single crystal structures (ection 5.1). Obviously, the new dcv4t thin film crystal structure (second column) is very different from the dcv4t single crystal structure (first column). Interestingly, however, it is remarkably similar to the dcv4t-m single crystal structure (third column) Charge Carrier Mobility The identified dcv4t thin film crystal structure is used to model a mesoscopic molecular crystal and examine its charge transport properties by microscopic charge transport simulations (ection 4.2). The starting point are molecular dynamics simulations performed on a morphology of 2880 molecules (with periodic boundary conditions), prepared bya fold translation of the crystal repeating motif. This is followed by evaluating the charge transfer parameters (i.e., the reorganization energies, site energy differences, and transfer integrals) and the Marcus rates for neighboring molecules.

131 130 Chapter 5. Charge Transport imulations in Organic Crystals dcv4t single crystal dcv4t thin film dcv4t-m single crystal µ µ µ µ µ µ exp Table 5.4. Eigenvalues of the zerofield mobility tensor, i.e., values for the principle axis directions, and values for the direction perpendicular to the dcv4t and dcv4t-m thin film substrates. Also provided are experimentally measured organic-field effect transistor mobilities µ exp [203]. All values are given in units of cm 2 /Vs. Finally, kinetic Monte Carlo simulations are employed to calculate the zero-field mobility tensor and a set of directional mobilities with electric field dependence. Unless otherwise stated, all methods and parameters employed are analogous to those used in the previous study on single crystals (ection 5.1). Before linking the efficiency of charge transport to the microscopic material properties, the simulation results of the dcv4t thin film are summarized. The energetic disorder, i.e., the width of the site energy difference distribution, is σ = 0.10 ev, which is identical to that of the dcv4t single crystal (ection 5.1). The topology of the charge transporting network, i.e., the three-dimensional assembly of the transfer integrals, and the zero-field mobility tensor are displayed in Figure 5.16 d, e (second column) on Page 132/133. The eigenvalues of the zero-field mobility tensor, i.e., the mobilities along the ellipsoidal principle axes, are provided in Table 5.4, together with values for the previously studied single crystal structures. The bounding eigenvalues define the mobility ranges indicated by the gray areas in the left panels of Figure The right panels show corresponding field-dependent mobilities. A first important observation can be made by comparing the charge transport in the two different crystal phases of dcv4t, i.e., the single crystal (first column) and the thin film (second column). The dcv4t single crystal exhibits a charge transporting network characterized by a distinct transport direction with strong electronic couplings (red), as discussed previously (ection 5.1). Counterintuitively, this network topology turned out to be disadvantageous. On the one hand, the strong electronic couplings in one direction lead to poorly coupled perpendicular directions and suppress couplings in other directions completely, and on the other hand, the one-dimensional strongly coupled direction is prone to defects due to the substantial energetic disorder. Contrary, the dcv4t thin film has a network topology with less pronounced one-dimensional character, hence the system is less sensitive to the energetic disorder. This has an imme-

132 5.2. Dicyanovinyl-ubstituted Quaterthiophenes: Thin Films 131 dcv4t single crystal dcv4t thin film dcv4t-m single crystal Mobility (µ/ cm2 ) Vs Electric Field (F/ V cm )1/2 Electric Field (F/ V cm )1/2 Electric Field (F/ V cm )1/2 Figure Left panels: ranges of the zero-field hole mobilities, bounded by the largest and smallest eigenvalues of the mobility tensors. Zero-field mobilities for the direction perpendicular to the dcv4t and dcv4t-m thin film substrates are indicated by red stars. Experimentally measured organic field effect transistor mobilities are indicated by black stars [203]. Right panels: hole mobilities as a function of an electric field applied along the largest and smallest principal axes of the zero-field tensors. diate implication on the charge carrier mobility: the maximum mobility of the dcv4t thin film is more than an order of magnitude higher than that of the dcv4t single crystal. It should be emphasized that this higher mobility is the exclusive result of the advantageous network topology since energetic disorder is identical in both systems. A second important finding follows from comparing the charge transport for the two different compounds with similar crystal structures, i.e., the dcv4t thin film (second column) and the dcv4t-m single crystal (third column). The dcv4t-m single crystal provides a charge transporting network characterized by multiple directions of moderate couplings, as discussed previously (ection 5.1). The longest principle axis of the mobility ellipsoid in Figure 5.16 e indicates that the main transport direction is a superposition of two completely symmetrical directions (both in red). Therefore, a charge carrier moving in a plane spanned by these two directions has two equally coupled hopping options at each site. This truly two-dimensional freedom of hopping allows for an effective bypassing of energetic defects. In contrast, the dcv4t thin film has its main transport direction along a superposition of two directions with asymmetric coupling (red and blue). Hence, this network topology adopts a hybrid character: although it is not strongly one-dimensional, it is not truly two-dimensional either. As a consequence it is more affected by the energetic disorder. This is reflected in the charge carrier mobility: the maximum mobility of the dcv4t-m single crystal is nearly an order of magnitude higher than that of the dcv4t thin film.

133 132 Chapter 5. Charge Transport imulations in Organic Crystals dcv4t single crystal NC NC (a) 0 c b (b) a 0 c (c) a 0 c b (d) -y x Figure (a c) Unit cells and crystal repeating motifs. The orientations are chosen to emphasize the similarity between the two crystal structures of the dcv4t thin film and the dcv4t-m single crystal. (d) Charge transporting networks based on the transfer integrals. The transfer integrals between the hopping sites are shown as bonds connecting black spheres. The bond color indicates the crystallographic direction while the bond thickness reflects the electronic coupling strength. (e) Zero-field mobility tensors depicted as ellipsoids. The colored arrows indicate the crystallographic directions of the charge transporting network. Adapted from Reference [3]. (e) x -y z z

134 5.2. Dicyanovinyl-ubstituted Quaterthiophenes: Thin Films 133 dcv4t thin film NC NC dcv4t-m single crystal NC Me Me Me Me NC (a) a b c 0 0 c (b) a b b 0 a 0 (c) c b 0 a a b 0 c (d) x y y -z x z (e) x y y -z x z

135 134 Chapter 5. Charge Transport imulations in Organic Crystals For the aforementioned dcv4t thin film (second column) and dcv4t-m single crystal structure (third column) remember that the latter equals its thin film structure, charge transport is additionally analyzed in the respective direction perpendicular to the film substrate, i.e., the direction of interest in an organic photovoltaic cell (opvc). ince Bragg reflections in the gixrd patterns correspond to lattice spacings perpendicular to the substrate, each Bragg reflection is linked to a crystallographic direction of interest. The gixrd pattern of the dcv4t thin film, shown in Figure 5.13 (black curve), has only a single dominant reflection at 12.8, which is, however, a superposition of two reflections at 12.7 and 13.3 observable in the dcv4t powder pattern (red curve). The lattice spacing of the stronger and the weaker of these two reflections can be identified with the and crystallographic directions of the reconstructed dcv4t thin film structure, respectively. In other words, the dcv4t thin film comprises two preferred orientations of crystallites: a major one, defined by the direction, and aminorone,definedbythe102directionbeingperpendiculartothesubstrate.note that the direction corresponds to the a-vector in Figure 5.16 a c (second column). The zero-field mobilities along both directions are listed in Table 5.4 as µ and µ, respectively, and are indicated by red stars in the left panels of Figure The situation is very similar for the dcv4t-m thin film: again, the gixrd pattern is dominated by a superposed Bragg reflection, which can be identified with two crystallographic directions, here the and directions of the dcv4t-m single crystal structure [3]. Note that the former corresponds to the b-vector in Figure 5.16 a c (third column). The zero-field mobilities of both directions are analogously presented in Table 5.4 and Figure Comparing the results for both systems reveals a similar trend as observed for the tensorial maximum mobility: the dcv4t-m system has nearly an order of magnitude higher mobility perpendicular to the substrate than the dcv4t system. Finally, the calculated mobilities for the dcv4t thin film (second column) and the dcv4t-m single crystal structure (third column) which equals its thin film structure are compared to experimentally measured thin film mobilities. The experimental mobilities originate from organic field-effect transistors (ofet) [203] and are listed in Table 5.4 and indicated by black stars in the left panels of Figure As expected and previously discussed (ection 5.1), the experimental mobilities are systematically lower than the calculated ones since they are measured on polycrystalline samples containing grain boundaries, which are an impediment to charge transport. In fact, this mobility reduction is particularly strong in ofet devices, where charge transport occurs parallel to the substrate over long distances and thus many grain boundaries.

136 5.2. Dicyanovinyl-ubstituted Quaterthiophenes: Thin Films 135 In opvc devices, in contrast, charge transport occurs perpendicular to the substrate over short distances, which corresponds more closely to the simulated situation. Apart from the expected offset, the experiments and simulations agree in yielding a higher thin film mobility for the methylated dcv4t-m as compared to the bare dcv4t compound. Moreover, opvc devices based on the methylated dcv4t-m show power conversion efficiencies of 3.8%, while those based on the bare dcv4t exhibit 1.5% [203], which also correlates with the calculated mobilities.

137 136 Chapter 5. Charge Transport imulations in Organic Crystals

138 Chapter 6. Charge Transport imulations in Organic Mesophases Among the most notable materials for building organic photovoltaic cells (Chapter 1) is the combination of dicyanovinyl-substituted oligothiophenes (dcvnt) and buckminster fullerene (C 60 ) as a donor-acceptor pair. For instance, devices prepared from amethylateddcv5t-mandc 60 have shown power conversion efficiencies of 6.9% [28] and with proprietary derivatives record efficiencies of 12% have been established [27]. An important aspect for achieving these high efficiencies is the use of bulk heterojunction device architectures where the donor and acceptor compounds are mixed in a blend layer. However, contrary to pure layers, which typically show crystalline order, these blend layers often exhibit increased disorder [206, 210]. Due to the molecular geometries, such disorder particularly affects the dcvnt donor domain of the blend and deteriorates its hole transport capabilities. To study this effect, the previous charge transport simulations in organic crystals (Chapter 5) are now complemented by simulations in organic mesophases of dicyanovinyl-substituted thiophene (dcv1t) up to sexithiophene (dcv6t). This series of compounds has been experimentally synthesized [202, 211] and is depicted in Figure 6.1. dcv1t dcv3t dcv5t NC NC NC NC NC NC dcv2t dcv4t dcv6t NC NC NC NC NC NC Figure 6.1. Chemical structures of terminally dicyanovinyl-substituted oligothiophenes.

139 138 Chapter 6. Charge Transport imulations in Organic Mesophases 6.1. Dicyanovinyl-ubstituted Oligothiophenes: Amorphous/mectic The present charge transport simulation study covers on the one hand amorphous mesophases of dcv1t to dcv6t, that is, systems with completely disordered molecular alignment, and on the other hand a more ordered smectic mesophase of dcv6t. An important result for the series of amorphous systems is that increasing oligomer length leads to increasing energetic disorder and hence to decreasing charge carrier mobility (anti-correlation). The smectic dcv6t mesophase, however, exhibits not only a lower energetic disorder than the amorphous dcv6t phase, but also a significantly lower mobility (correlation). This finding is not only inconsistent but also contradicts the common belief that a higher mesophase order should promote charge transport. The microscopic origins of this inconsistency are elucidated by analyzing the energetic landscapes of site energies and their interrelations to the charge transporting networks of transfer integrals. The study is performed by microscopic charge transport simulations (ection 4.2). As a starting point, mesoscopic morphologies are generated as well as characterized with respect to their disorder (ection 6.1.1). This is followed by the evaluation of the charge transfer parameters between neighboring molecules (ection 6.1.2). Based on the resulting Marcus rates, charge carrier dynamics is simulated leading to the charge carrier mobilities (ection 6.1.3). The discussion concludes with the identification and analysis of electric current pathways (ection 6.1.4). The reported scientific results are the subject of prior publication [2], and are presented below in more detail using partially similar terms and illustrations Morphological Disorder The morphologies, comprising 4096 molecules, are prepared by first distributing rigid molecules with random orientations in boxes (with periodic boundary conditions) and then running molecular dynamics simulations (ection 4.2.2) in the isobaric-isothermal ensemble. For each compound of the dcv1t to dcv6t series, an amorphous mesophase is generated by starting with a preliminary equilibration process well above the glass transition temperature, performed at a temperature of 800 K, a pressure of 1bar,andforadurationof10ns,andtheninitiatingaquenchingstepto300K,after which a further equilibration for 10 ns takes place. A final room temperature morphology of the dcv4t system is illustrated in Figure 6.2.

140 6.1. Dicyanovinyl-ubstituted Oligothiophenes: Amorphous/mectic 139 Figure 6.2. Amorphous morphology of dcv4t molecules at room temperature, obtained by a molecular dynamics simulation. For the longest oligomer, i.e., the sexithiophene dcv6t, the system is isotropic at 800 K, but it spontaneously transitions into a smectic mesophase upon reducing the temperature to 700 K. In order to avoid having defects in the smectic layers, a nonrandom initial configuration is created where molecules are arranged in regular layers. This configuration is equilibrated at 700 K for 10 ns and then quenched to 300 K, where it is equilibrated for another 10 ns. The resulting molecular arrangement has the same spacing between the smectic layers as the one emanating from the isotropic phase. Molecular Conformations During the equilibration of the systems at high temperatures, cis-trans isomerization of the dcvntmoleculesoccursduetotherotationof dcv-thiophene (cc c ca s) and thiophene-thiophene dihedral angles (s ca ca s). This is a consequence of the energetic barrier between cis and trans states of these dihedrals, which is 0.45 ev and 0.15 ev, respectively, as determined within the force field development (ection 4.2.1). In fact, each amorphous or smectic system contains all possible molecular rotamers arising from the combinations of the two dihedral states. For a dcv4t snapshot, the rotamer distribution of the 4096 molecules is exemplified in Figure 6.3 (top panel). 1 This conformational, or rotameric disorder leads, due to the acceptor-donor-acceptor molecular architecture of dcvntmolecules,toastrong dipolar disorder. As can be seen in Figure 6.3 (bottom panel), the different conformers exhibit fluctuations in their permanent molecular dipole moments of up to 15 D, as estimated from the atomic partial charges. This significant dipolar disorder is closely related to energetic disorder, as will be discussed below. 1 Note that the incidence of rotamers is roughly anti-correlated with their internal site energy (4.29) [2].

141 NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC NC 140 Chapter 6. Charge Transport imulations in Organic Mesophases Count Dip. Mom. (d/d) Figure 6.3. Top panel: distribution of dcv4t rotamers in an amorphous morphology of 4096 molecules. Bottom panel: permanent molecular dipole moments of dcv4t rotamers, calculated from atomic partial charges, given in units of Debye. Adapted from Reference [2]. Molecular Orientational Order In order to illustrate the differences in molecular ordering of the two different dcv6t mesophases, i.e., the amorphous and smectic dcv6t mesophases, Figure 6.4 a visualizes the molecular backbones within slices of their morphologies. For all systems, the molecular orientational order is quantified by means of the orientational order parameter, which is zero for ideal isotropic and one for perfect parallel alignment of the molecules. The orientational order parameter s of a system is defined as the largest eigenvalue of its order tensor Q,givenbythecomponents Q µν = 1 2 3uµ u ν δ µν, (6.1) where the unit vectors u,withcomponentsu µ, point along the direction of the molecular backbones and the averaging is performed over all molecules of the system. Here, the unit vector is defined by the two branching carbons of the dcv groups (c), as defined in the force field (ection 4.2.1). For the amorphous systems of dcv1t to dcv6t, the evaluation of the order tensor (6.1) yields values s < 0.1, which confirms almost complete isotropy of the molecular orientations. For the smectic mesophase of dcv6t, a value of s > 0.8 is obtained, which reflects the strong orientational order. The preferred molecular orientation in the smectic mesophase is along the x-axis, as follows from the director, that is, the eigenvector associated with the eigenvalue s. The exact order parameter values s are summarized in Table 6.1, together with the mass densities d of all systems.

142 6.1. Dicyanovinyl-ubstituted Oligothiophenes: Amorphous/mectic 141 dcv6t dcv6t smectic (a) z z x x y (b) 1eV 1eV z z x 1eV x y 1eV (c) z z x x y Figure 6.4. (a) Morphologies of amorphous and smectic dcv6t mesophases at room temperature. The illustrations show molecular backbones within 2 nm thick slices oriented perpendicular to the y-direction, as well as to the x-direction for the smectic mesophase. (b) Cross-sections of the site energy landscapes, visualized as density plots. The underlying hopping sites (molecular centers of mass) are indicated by black dots. (c) Charge transporting networks based on the transfer integrals. The transfer integrals between the hopping sites are shown as bonds connecting black dots. The bond color and thickness reflects the electronic coupling strength. Adapted from Reference [2].

143 142 Chapter 6. Charge Transport imulations in Organic Mesophases dcv1t dcv2t dcv3t dcv4t dcv5t dcv6t dcv6t smectic s d Table 6.1. Orientational order parameters s, obtained from the orientational order tensors, and mass densities d,giveninunitsofg/cm Charge Transfer Parameters To construct the charge transport Hamiltonians, the morphologies are partitioned on conjugated segments and rigid fragments (ection 4.2.3). Accordingly, the dcvntmolecules are represented by separate rigid fragments for each dcv or thiophene unit, and by single conjugated segments, acting as the hopping sites i. The neighbor list of hopping sites i and j contains molecular pairs, subject to the condition that the distance between any pair of associated mutual rigid fragments is below a cutoff of 0.8 nm. Reorganization Energies Evaluating the reorganization energies (ection 4.2.6) is performed by dft calculations on isolated molecules using the b3lyp functional and the 6-311g basis set. ince the systems under consideration are homogeneous, a universal reorganization energy λ is determined for each compound. For the series of dcv1t to dcv6t, the resulting values range from 0.17 ev to 0.23 ev, as is summarized in Table 6.3 on Page 148. In order to estimate the error resulting from the neglect of geometry variations in the bulk, reorganization energies are calculated for all different dcv4t rotamers yielding a standard deviation of 0.01 ev, which is negligible. ite Energies The observed fluctuations in molecular dipole moments are a clear indication of variations of the site energies. To quantify this energetic disorder in the mesomorphic systems, the site energies E i of all molecules i are explicitly calculated from the electrostatic interaction energy including polarization effects. To this end, electrostatic interactions are determined self-consistently between atomic partial charges as well as induced atomic dipole moments (enabled by atomic polarizabilities) using the Thole model (ection 4.2.5). To achieve viable runtimes for these calculations, a spherical interaction cutoff of 3.5 nm is employed. For the amorphous and smectic mesophases of dcv6t, cross-sections of the energetic landscapes are visualized in Figure 6.4 b. Here, the range from low to high site energies is indicated by a color range from blue to red. As can be seen, the site energies are spatially correlated, which is a result of the long-range nature of electrostatic interactions.

144 6.1. Dicyanovinyl-ubstituted Oligothiophenes: Amorphous/mectic 143 dcv6t z dcv6t smectic z x y x Figure 6.5. Direction-resolved distributions of neighboring hopping sites. This sketch corresponds to the set of endpoints of all coupling vectors shifted to the origin. Adapted from Reference [2]. For all studied systems, the distributions of site energy differences, E ij = E i E j, for molecular pairs from the neighbor list are displayed in Figure 6.9 a on Page 148. The widths of Gaussian functions fitted to these distributions yield the neighbor-listbased values of energetic disorder σ, as defined in (4.40), which are listed in Table 6.3. Obviously, all systems possess substantial energetic disorder, which increases with the number of thiophene units per molecule. As might be expected, the more ordered smectic mesophase of dcv6t exhibits a lower energetic disorder than the amorphous phase of dcv6t. Transfer Integrals The transfer integrals J ij between pairs of molecules i and j from the neighbor list are determined by means of the zindo method (ection 4.2.4). The topological graph defined by the transfer integrals represents a percolating network for charge transport which is characteristic for each system. For the amorphous and smectic mesophases of dcv6t, representative slices of the charge transporting networks are visualized in Figure 6.4 c. Here, hopping sites (molecular centers of mass) are drawn as black spheres, while transfer integrals are represented by intermediate bonds, whose color and thickness reflect the magnitude of the electronic coupling strength. As can be seen, the amorphous mesophase has a spatially uniform distribution of hopping sites. Contrary to this, the smectic mesophase comprises a set of twodimensional layers parallel to the yz-plane with strong intra-layer, but weak inter-layer couplings. For both systems, Figure 6.5 illustrates the direction-resolved distribution of neighboring hopping sites, which corresponds to the set of endpoints of all coupling vectors shifted to the origin. For a charge carrier located at the central site, the average set of accessible sites is a spherical shell in the case of the amorphous dcv6t, and aringwithtwoislandsinthecaseofthesmecticdcv6t. Whiletheringrepresents the hopping options within the same layer, the islands correspond to the neighboring layers. Apparently the layer formation in the smectic system leads to a dimensionality reduction of the charge transporting network.

145 144 Chapter 6. Charge Transport imulations in Organic Mesophases dcv1t dcv2t dcv3t dcv4t dcv5t dcv6t dcv6t smectic µ σ= / 0.11 µ Table 6.2. Zero-field mobilities µ σ=0 and µ, calculated in the absence and presence of energetic disorder, respectively. For the smectic dcv6t phase, it is distinguished between the directions along the x-axis and within the yz-plane (first/second value). This distinction is waived in the presence of energetic disorder, where the mobility becomes isotropic. All values are given in units of cm 2 /Vs. For all systems, the distributions of transfer integrals are displayed in Figure 6.9 c on Page 148. The average values J 2 are listed in Table 6.3, together with the average center-of-mass separations r. For the smectic dcv6t system, the inter-layer and intra-layer transfer integrals are distinguished between (first/second value reported). Ageneralobservationisthattheaverageelectroniccouplingstrengthisanti-correlated with the average separations. For the amorphous dcv1t to dcv6t systems, increasing molecular chain length causes an increase in separations and a decrease in electronic couplings. For the smectic dcv6t system, there are smaller separations and larger electronic couplings within the layers than in between layers Charge Carrier Mobility Based on the charge transfer parameters, the charge carrier dynamics can be studied. In order to explore the effects of energetic disorder on charge transport, the zero-field mobility is evaluated for each system twice. First, energetic disorder is turned off and the zero-field mobility µ σ=0 is calculated from charge diffusion without external fields. econd, energetic disorder is included and the zero-field mobility µ is extrapolated from the field-dependent charge carrier mobility, determined from charge drift-diffusion under externally applied electric fields. Note that an external field is required to accelerate the charge carrier motion since energetic disorder is particularly strong. Technically, the evaluation of charge carrier mobilities is based on kinetic Monte Carlo simulations (ection 4.2.7) and the relations to macroscopic observables (ection 4.2.8). Without an external field, averages of 2000 charge displacements for diffusion times of 10 3 sarecarriedout.whenapplyinganelectricfield,itsmagnitudeisincreasedin steps of 100 from 200 to 1000 (V/cm) 1/2.Foreachfieldmagnitude,anaveragingover independent time intervals (of 0.1 s for dcv1t and dcv2t, 1 s for dcv3t and dcv4t, 10 s for dcv5t and dcv6t, and 10 5 s for the smectic dcv6t) for 14 field directions (along and against the axes and main diagonals of the simulation box) is performed.

146 6.1. Dicyanovinyl-ubstituted Oligothiophenes: Amorphous/mectic 145 Mobility (µ/ cm2 ) Vs dcv1t dcv6t dcv6t smectic Figure 6.6. Right panel: hole mobilities as function of an external electric field, plotted in a logarithmic scale and given as a square root, respectively (Poole-Frenkel chart). Linear fits are employed to extract the zero-field mobilities. Left panel: zero-field mobilities. Adapted from Reference [2]. Electric Field (F/ V cm )1/2 In the absence of energetic disorder,the zero-field mobilitiesof all systemslieinthe range of 0.40 to 0.11 cm 2 /Vs, as summarized in Table 6.2. Hence, the mobility is almost independent of both the molecular length and the mesophase ordering. When analyzing the transfer integrals, it was found that the electronic coupling strength anticorrelates with the average separations between molecular centers of mass. ince this relationship is roughly linear, variations in these two quantities balance each other in the expression for the charge carrier mobility tensor (4.48), leading to uniform mobilities. It is particular noteworthy that this applies equally to the smectic system, which does not show a substantially increased mobility as compared to the amorphous systems, not even within the smectic layers. This is due to the still disordered, i.e., liquidlike, molecular arrangement in the smectic layers. One can generalize this finding by stating that the charge transport efficiency depends on the local molecular ordering, e.g., the crystallization within smectic layers, just as on the long-range ordering, e.g., the formation of a smectic mesophase. In the presence of energetic disorder, thechargecarriermobilitiesdependexponentially on the square root of the external electric field F,asisseenintherightpanelof Figure 6.6. This relationship, known as the Poole-Frenkel effect, is used to extrapolate the zero-field mobilities, provided in the left panel as well as in Table 6.2. In the amorphous systems it is expectedly found that the strong energetic disorder leads to a mobility reduction by several orders of magnitude, along with a systematic mobility decrease as the energetic disorder increases. This is in agreement with the results of Gaussian disorder models (ection 4.1.2), where higher energetic disorder leads to lower charge carrier mobilities. The situation is, however, markedly different in the smectic mesophase of dcv6t. Here, the reduction in mobility is exceptionally large, resulting in values substantially lower than those of the amorphous dcv6t, in spite of the lower energetic disorder. This inconsistency, apparent in Figure 6.10 a on Page 149, contradicts

147 146 Chapter 6. Charge Transport imulations in Organic Mesophases Probability Edge Current (c ij / e m s ) dcv1t dcv2t dcv3t dcv4t dcv5t dcv6t dcv6t smectic Figure 6.7. Distributions of absolute values of edge currents. Each filled area indicates the partial distribution of highest edge currents which contribute to 90% of the total current. Adapted from Reference [2]. both the Gaussian disorder models and the common belief that higher mesophase order should improve charge transport efficiency. In order to understand this behavior, the microscopic channels of electric current flow are analyzed below Electric Current Pathways To identify the preferred channels of charge flow in a system, the first step is the evaluation of the edge currents c ij, as defined in (4.46), for all molecular pairs i and j from the neighbor list. ubsequently, the absolute values c ij are sorted and the molecular pairs with small currents are removed until the remaining sum of currents reaches 90% of the original total current. In this way, a subset of the neighbor list is created, containing only the molecular pairs which conduct 90% of the total current, i.e., that contribute to the charge transport to an appreciable extent. This subset of pairs, which is actually used by a drift-diffusing charge carrier, is only a small fraction of the total neighbor list, as can be seen in Figure 6.7. Here, the contributing edge currents are indicated by a filled area of the total current distribution. These contributing edge currents are visualized in Figure 6.8 a for an amorphous dcv4t morphology, where it is clearly seen that they form filamentary pathways of electric current. The significant reduction of molecular pairs participating in charge transport can be attributed to the spatial correlations of site energies, caused by the long-range nature of electrostatic interactions. These correlations lead to extended areas of low site energies in which the charge carrier migration takes place. To illustrate this behavior, Figure 6.8 b shows the current filament in combination with a cross-section of the site energy landscape. The filament is clearly percolating within the energetic valleys, colored in blue, while avoiding the energetically unfavorable regions, colored in red.

148 6.1. Dicyanovinyl-ubstituted Oligothiophenes: Amorphous/mectic 147 (a) (b) z z x y x y Figure 6.8. (a) Edge currents contributing to 90% of the total current in an amorphous dcv4t system. An electric field of 1000 (V/cm) 1/2 is applied in the z-direction. The arrowheads indicate the direction of the current while the thickness and color of the arrowshafts reflect the current amplitude. (b) Current filament and cross-section of the site energy landscape. Adapted from Reference [2]. On the basis of the reduced neighbor list, which contains the molecular pairs k and l conducting 90% of the total current, the distributions of the site energy differences, E kl = E k E l, are evaluated again. These distributions are displayed for all systems in Figure 6.9 b. The widths of Gaussian functions fitted to these distributions provide values of an effective energetic disorder σ eff,whicharesummarizedintable6.3. Evidently, this effective energetic disorder is significantly lower than the previously determined energetic disorder σ,which was evaluated based on the entire neighbor list. Hence, a first conclusion is that, even for amorphous systems, the characterization of energetic disorder by site energies alone can be misleading.