Charge Transport in Organic Materials: Norm-Conserving Imaginary

Subscriber access provided by University of Florida | Smathers Libraries

Article

Charge Transport in Organic Materials: Norm-Conserving Imaginary Time Propagation with the Local Ionization Energy as External Potential Maximilian Kriebel, Dmitry Sharapa, and Timothy Clark J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00568 • Publication Date (Web): 19 Oct 2017 Downloaded from http://pubs.acs.org on October 21, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

1

Charge Transport in Organic Materials: NormConserving Imaginary Time Propagation with the Local Ionization Energy as External Potential Maximilian Kriebel, Dmitry Sharapa, and Timothy Clark*

Computer-Chemie-Centrum and Interdisciplinary Center for Molecular Materials, Department of Chemistry and Pharmacy, Friedrich-Alexander-Universität Erlangen-Nürnberg, Nägelsbachstrasse 25, 91052 Erlangen, Germany.

Abstract: An additional charge carrier described as its wavefunction is propagated in imaginary time using stepwise matrix multiplication and a correction to ensure that the simulation is norm conserving. The propagation Hamilton operator uses the local ionization energy of a rubrene single crystal, calculated with semiempirical MO-theory, as an external potential for holes to model the interaction with the underlying molecular structure. Virtual electrodes are modeled by setting the potentials in the appropriate areas to constant values with the difference corresponding to the source-drain voltage. Although imaginary time cannot be interpreted directly as time, the simulated gate-dependent imaginary transfer rate is in acceptable qualitative agreement with the experimentally measured gate-dependent hole-transfer rate through a rubrene single crystal.

E-Mail: [email protected] ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 30

2

1. INTRODUCTION Organic and hybrid inorganic/organic electronic devices and photovoltaic cells represent a promising alternative to conventional solid-state semiconductor technology. Among others, their advantages include cheap and easy fabrication, often by self-assembly, and low power consumption. A major advantage is that a wealth of experience in organic synthesis can be used to assemble the component molecules. For instance, molecules for self-assembled monolayers (SAMs) can be designed with anchoring, isolating and conducting groups and their self-assembly properties used to build devices quickly and economically.1 Similarly, organic molecules can be designed specifically for ink in printable electronic devices2 to build flexible electronic devices or solar cells by attractive roll-to-roll printing processes. Additionally, organic crystals offer one of the few opportunities to investigate features such as the direction anisotropy of charge-carrier mobility.3 Transport mechanisms in solids are sometimes categorized according to the way they are treated theoretically. Delocalized band transport in which, for example, Monte-Carlo algorithms in momentum space are used,4 or the average mean free path of a charge carrier before an elastic scattering event occurs5 can be evaluated. Band transport considers global effects and is therefore able to describe tunneling and long-range resonances. This contrasts with localized hopping transport, whose treatment in its different variations is often based on Marcus theory.6 Hopping transport treatments are well suited for local anomalies or micro-structured geometries and generally less well for homogeneous systems.1,7 Experiment8 suggests charge carriers in organic molecular crystals to be delocalized over several molecules. Charge delocalization in organic aggregates in general seems to depend on the order of the system but charge delocalization over multiple molecules has been observed in several cases9 and has 2

ACS Paragon Plus Environment

Page 3 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3

been demonstrated in solution10 for C60-fullerene van der Waals dimers. Therefore, the charge-carrier behavior for organic crystals can neither be described adequately by fully delocalized band transport nor by fully localized hopping transport.11 This is sometimes called the intermediate regime between band and hopping transport and can be targeted by propagating the excess charge-carrier wavefunction, for example with the mean-field dynamics approach, which is based on the Ehrenfest theorem,12 or by Tully’s surface-hopping method.13-15 The applicability of these methods, however, is currently limited by their high computational cost. Other established methods commonly used to describe charge transport are non-equilibrium Greens functions (NEGF)16 also in a time dependent manner (TDNEGF)17,18 and Landauer transport.19,20 In the hopping picture, small finite temperatures increase the probability of charge carriers escaping from traps21 and therefore increase the overall mobility. If the temperature becomes too high, the electron-phonon coupling increases,

21,22

lowering the mobility. However, weak electron-phonon

coupling leads to strong charge delocalization, so that the temperature dependence of charge transport makes the choice of band or hopping regime even more difficult. Putting the nuclear and electronic dynamics on the same footing by explicit time dependency, as mixed quantum-classical dynamics (MQCD) aims to do, represents a further step in the right direction.23 Therefore, charge transport in organic materials, or more generally in inhomogeneous molecular aggregates or polymers, is not as well understood as that in conventional solid-state devices. Because of this, the design and construction of high-performance organic field-effect transistors (OFETs) with high mobility and fast on/off ratios is currently predominantly a matter of trial and error. Rational, knowledge-based design requires a fast and easy way to determine the charge-transfer properties of layers and aggregates in order to characterize the transfer of energy by the flow of charges or even the evolution of single charge carriers. 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 30

4

It is not easy to access the intrinsic charge-conduction properties, neither experimentally nor by simulation, for several reasons. Another very important factor that affects the measured device performance is the interface between the electronically active structure and the electrode, a subject of growing interest.24,25 Experimental mobilities of single-crystal OFETs also depend strongly on the purity,11,21 of the crystal. In microcrystalline OFETs, the complete morphology of the conducting layer is even more important21,26 because inhomogeneities can interrupt charge-transport pathways1 or trap charge carriers at grain boundaries.26,27 A direct comparison of charge-carrier mobility measurements on single crystals with those on amorphous thin-film transistors is therefore not suitable for model validation. Most practical applications of OFETs rely on amorphous structures. Because, however, so many external factors can influence measurements on such structures, we focus here on the better-defined situation in single crystals. Our ultimate purpose, however, is to simulate charge-transport in amorphous, microcrystalline or defect-rich structures. The geometry of the electrodes used28 and even the insulating layer29 also influence experimentally determined mobilities strongly, both for self-assembled devices and for single crystals. In the latter, the band gap itself is more important at the interface of the electrodes applied to the crystal and for inducing electron-hole pairs in the crystal, while the shape of the conduction band is

1 4


Page 5 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5

most important for charge propagation inside the crystal. Rubrene, 1, single crystal organic field-effect transistors (OFETs) have been widely studied experimentally.30-33 Although they are influenced by many parameters, the reported gate-dependent mobility measurements34 of these crystals provide a first opportunity to validate the charge-carrier propagation technique described below. The propagation Hamiltonian contains the interaction with the structure in the form of a local molecular property, the local electron affinity35,36 for electrons and, in

Figure 1: View of a rubrene single-crystal structure showing the local ionization energy at the left with an isovalue surface of 400 kcal mol−1 in blue and the local electron affinity on the right with an isovalue of -40 kcal mol−1. The crystal axis from left to right corresponds to the y-axis of the simulation space. The charge-transport channels are much more evident for the local ionization energy, which can be thought of as a potential for holes Vh( xr ) = I l ( xr ) , than for the local electron affinity, which can be thought of as potential for electrons Vh( xr ) = − EAl ( xr ) . The figures indicate that Rubrene is a better hole- than electron-conductor.

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 30

6

this work, the local ionization energy37 for hole-propagation. These local properties represent an approximation in their present form, 38 as they are strictly speaking not pure potentials but also contain a kinetic-energy contribution. We, however, have previously used them as potential hypersurfaces for a multi-agent quantum Monte Carlo model for charge transport with some success to simulate an organic self-assembled monolayer system.1,39 The simulated charge-carrier pathway densities found in this work agree with the ones calculated for the same system using the Monte-Carlo technique. Using these 3Denergy maps as external potentials in the Hamilton operator allows the interaction of charge carriers with a complicated quantum mechanical system to be represented as simple scalar potentials embedded in three-dimensional space. The unique aspect of this work is to combine the use of these 3D energy potentials in combination with the propagation of a single charge carrier wavefunction itself directly in 3D real space without a stochastic component. The essential point of this approach is that it reduces the number of dynamic objects to be treated but can still reveal information about the conduction characteristics of the system. The approximation of treating charge transport against the background of a geometrically and electronically frozen substrate is, of course, drastic. Nonetheless, we believe that the technique described can provide useful semi-quantitative information at low computational cost and form the basis for more sophisticated treatments. In this work, we focus on the injection and extraction of a single hole to and from the rubrene crystal and its transport within the crystal itself. The calculations do not yet represent simulations of complete electronic components.

6


Page 7 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


7

Figure 2: The local ionization energy isosurface (400 kcal mol−1) of a rubrene crystal is shown in opaque red and the probability density of a hole (2.5 × 10-6) in a given timeframe snapshot in blue. The molecular orbitals and their energies are generated with Neglect of Differential Diatomic Overlap (NDDO) based semiempirical Molecular Orbital (MO) theory, which allows very large systems to be calculated without local approximations.40-42 The local ionization energy38 is calculated from the NDDO wavefunction. Hole-transport paths are visible simply by plotting isovalue surfaces of the local ionization energy, as shown in Figure 1, or by projection onto an isodensity molecular surface.43 Imaginary time evolution is usually used to find the ground state of quantum systems, for example by solving static Hartree Fock problems,44 or evolving real-space wavefunctions45 to their ground state. We now describe a technique for calculating the motion of charge carriers using the imaginary time evolution of an unobserved hole described as a wavefunction in real space in order to gain insight into the charge injection, transport and extraction of a molecular crystalline system. An example is shown in Figure 2. This technique is analogous to the multi-agent quantum Monte-Carlo methodology described earlier39 but treats a single charge carrier and provides a time axis, albeit imaginary. The convergence of pathway calculations with previous work that used the local electron affinity as external potential in 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 30

8

Monte Carlo charge-carrier hopping1,39 and with direct wave propagation on a self-assembled monolayer system will be the content of future work.

2. METHOD An additional hole associated with a molecular system is described as a wavefunction in real space

r Ψ ( x , t ) and evolved with imaginary time propagation. The interaction of this additional charge carrier with the neutral molecular structure with its electrons and nuclei is described using suitable local properties as external potential terms in real space in the propagation Hamiltonian of the charge-carrier wave function. In atomic units, e = m = h = 4πε 0 = 1 , the time-dependent Schrödinger equation describes the change in time of a quantum system:

r r r i∂t Ψ ( x , t ) = Hˆ ( x ) Ψ ( x , t )

(1)

r r Where Ψ ( x , t ) describes the wavefunction of the hole, Hˆ ( x ) the Hamiltonian of this hole on the system, ∂ t the time derivative operator and i the complex unit. This leads to an operator that allows the wavefunction after a small time-step dt to be calculated: (2)

If now dt is chosen to be imaginary in the following way:

(3)

,which is called Wick rotation,46 the operation (2) becomes real and is used in this work. This is usually known as imaginary time propagation and is mathematically strongly related to the Diffusion Equation:47 8


Page 9 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


9

(4)

The interaction partner described by the wave is the hole, the other interaction partner is the molecule. In this work the atomic structure is completely fixed and no explicit phonon treatment is performed. It therefore cannot take temperature dependence, which manifests itself through movements of the atomic grid, into account. However, we do not exclude the possibility of considering phonon coupling in later work. All molecular orbitals are based on restricted AM148 self-consistent field wavefunctions for the neutral molecule, which are calculated using the NDDO approximation.49,50 The AM1 calculations use a

r

minimal basis set of Slater orbitals. The local ionization energy IL is given37 at position x by:

r IL ( x ) =

−

HOMO

r

∑

ε n ρn ( x )

∑

r ρm ( x )

n =1 HOMO

(5)

m =1

Where HOMO is the highest occupied molecular orbital, εn is the eigenvalue of molecular orbital n

r r and ρ n ( x ) is the electron density at position x . The Hamiltonian for an additional hole on a molecular structure is then given by: (6)

The hole potential in the electrodes is modeled as follows using the source-drain voltage, VSD , (see also Figure 2): (7)

(8)

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 30

10

Because the propagation is performed in real space, not only the potentials, but also the starting wavefunction itself are arbitrary in the sense that the complexity of the calculations does not have to be reduced artificially by assumed symmetries or periodicities. Also, no split-operator technique is used. Therefore, numerical application of the double derivative operator on the wave is necessary to calculate the kinetic energy. The discretization scheme is forward in time and central in space. With the number of grid points Nx in the x-direction, Ny in the y-direction and Nz in the z-direction, the second spatial derivative operators of the respective axis are calculated with a central-in-space finite-difference scheme. The first derivative is constructed by adding intermediate grid points: (9)

The second derivative then does not depend on the intermediate grid points:

(10)

With dx = 1 leads to the double derivative matrix:

(11) Because these double-derivative operators are used in the program directly, the exponential function is truncated after the third order to reduce the numerical error and to allow larger time steps to be used:

10


Page 11 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


11

(12) Then, the imaginary time propagation of the hole-wavefunction is given by:

(13) To conserve the norm, the imaginary time propagation operation is modified as follows: (14)

Which can easily be shown, beginning with a normalized distribution: (15)

The norm of the wave propagated a small imaginary time step

is given by:

(16)

While for small

:

(17)

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 30

12

Imaginary time propagation is not energy conserving. For a single particle with a Boltzmann-like energy distribution, this operation can be seen as a temperature reduction, where d%t is connected to a temperature change when assuming a Boltzmann distributed momentum density that depends on a temperature :

with

(18)

with

(19)

While

The wavefunction can be written as: (20)

, which is equivalent to: (21)

if the operation e

r − dt% Hˆ ( x )

%ˆ

e − dtH represents a change in temperature. The resulting wavefunction after the

operation is expected to be Boltzmann distributed and to depend on a new temperature

:

(22) Because we know that operation (11) is norm conserving, it then must be true that:

12


Page 13 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


13

(23)

and with: (24)

or: (25)

One can calculate the new temperature

:

(26)

and the change in temperature, which depends on the current temperature

and the imaginary time dt% :

(27)

A second derivation is given in the appendix. With an applied electric field accelerating a charge-carrier, mobility is defined by the maximal average drift of the charge carrier when an equilibrium is reached, in which the charge carrier transforms as much potential energy to kinetic energy as it loses through interaction with the system (phonon scattering) and radiative cooling (photon scattering). In our case, we use the energy loss inherent in imaginary time propagation as a surrogate for the radiative energy loss mechanism. Whether it can also serve as a surrogate for phonon scattering depends on whether there is a correlation between the

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 30

14

underlying potential, which is given by the electronic structure, and the coupling between the charge carrier and lattice vibrations. An "imaginary" current for the final imaginary time period t% is calculated by summing over the density given by the absolute square of the hole-wavefunction in the drain electrode: (28)

We refer to Q& imag as imaginary current because we use imaginary time in the denominator. While in this work: (29)

is denoted the imaginary mobility, where d is the distance traveled, which is calculated from the center of the modeled source to the center of the modeled drain electrodes.

3. SIMULATIONS The molecular orbitals (ϕn , ε n ) were calculated using the restricted Hartree-Fock (RHF) formalism on molecular clusters using the standard AM1 Hamiltonian48 with the massively parallel semiempirical MO-program EMPIRE.40,41 Periodic SCF calculations to obtain the electronic structure of the neutral rubrene crystal were performed without an applied gate field. Calculations were based on the X-ray structure of rubrene obtained from the Crystallography Open Database, COD ID 2100410.3 A supercell containing 6,720 atoms was generated (2 × 4 × 3 unit-cells or 53.578 × 28.680 × 42.633 Å). Localproperty grids (local electron affinity EAl ( xr ) 35,36 and local ionization energy I l ( xr ) 37) were generated for a 50 × 250 × 50 Å box (spacing 0.5 Å) using the in-house eh5cube utility. Thus, the simulation space 14


Page 15 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


15

contains more than 200 unit cells. Isovalue plots of I l ( xr ) and EAl ( xr ) for the rubrene crystal are shown in Figure 1. In general, one can assume that a local property that enables more continuous conduction paths and less traps has a positive influence to the conductivity of the structure.1,51,52 Determined by the calculated resolution of the underlying local-property data file, the simulation space on which the propagation is performed consists of a cubic cell of 100 × 500 × 100 grid points with a spatial grid for each simulation. In this work, the electrodes are modeled by setting their potentials to constant values as follows:

A hole as a charge carrier located somewhere in the source electrode without assuming a concrete localization is represented by the starting wavefunction Ψ( xr,t=0) as a function with a constant value in the source electrode, and zero everywhere else. The gate effect is modeled by shifting the local ionization energy between the electrodes (VGate in Figure 3). The geometry on which the calculation of the local ionization energy is based describes a perfectly pure, defect-free single crystal. Using the EMPIRE program, however, would also allow us to consider defects or impurities in very large periodic cells.40,41,53 Potential differences from source to drain are approximately 5, 10 and 15 V. We carried out the simulations with a custom Fortran 90 program that reads the grid of local-property values as a cube file.54 The propagation involves stepwise matrix multiplications, which can easily be performed on multiple cores.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 30

16

Figure 3: Schematic diagram of the potentials and the differences between them. The source potential is raised by the half of the source-drain potential, VSD and lowered correspondingly for the 2

drain. The gate voltage is modeled in the simulations by shifting the hole potential on each grid point between the electrodes by a constant value, Vgate . The averaged rubrene potential for a chargecarrying hole I L is 17.47V.

VRS is the potential difference between the hole source and the averaged

rubrene system local ionization energy. −VRS is used for the horizontal axis in Figure 4(B) in order to be able to compare the simulated results with experimental ones more easily. Figure 4 Shows the results of a grid convergence study in which simulations for the conditions that gave the highest mobility were performed in simulation grids of the same dimensions but different resolution (200×1000×200, 100×500×100, 50×250×50, and 25×125×25 with grid spacings of 0.25, 0.5, 1.0 and 2.0 Å, respectively). All other parameters of the calculation were kept fixed. The Figure shows that the results depend on the grid spacing but that they are essentially converged at spacings smaller than 0.5 Å.

16


Page 17 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


17

Figure 4: A grid convergence analysis with varying grid spacing was performed with all other parameters fixed. The calculated mobilities converge at small grid spacing.

4. RESULTS AND DISCUSSION Figure 5(A) shows plots of the calculated imaginary hole-mobility as a function of Vg as obtained from the simulations. Figure 5(B) shows imaginary hole-mobility as a function of VSG(VG). The sharp onset occurs at the correct source drain voltage in the simulations. The curves exhibit the correct general form but, in contrast to the experimental curves shown in Figure 5(C), the height of the mobility maximum decreases and the gate voltage at which it occurs shifts to more negative values with increasing sourcedrain voltage. The plots shown in Figure 5(B) are those from Figure 5(A) scaled so that the positions of the mobility maxima coincide, as in the experiment. Despite the quantitative differences, our simulations reproduce the small dependence of the single-crystal mobilities on the gate field at strongly negative gate potentials and the very sharp field-effect onset.26,34 The peak positions and shapes are close to experiment. Apart from this shift, the peak becomes flatter and broader and the onset becomes flatter, both in agreement with experiment. The maximum imaginary mobility values decrease for higher source-drain fields, which suggests that the ratio of imaginary time passed in comparison to real time depends on the potential. In other words, the rate of energy loss in the imaginary time propagation increases too strongly with higher source-drain field strengths. The overall mobilities in experiments are also influenced by impurities,11,21 the metal organic interface28 and even the insulating layer.29 Therefore, experimental mobilities are expected to be lower than simulated. Nevertheless, the qualitative agreement of the shapes of the curves with experiment drew our interest. Real-time propagation is energy conserving and because the mobility of a system is characterized by a saturated charge carrier 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 30

18

velocity, it would be not able to describe the charge-transport process considered here. Here, the simulations were prepared in a manner that fits well to the experimental situation of ref. 34 by modeling r the electrodes and choosing the starting wavefunction Ψ ( x , t% = 0 ) as a constant function inside the hole drain electrode. Rubrene, however, also exhibits direction dependent mobilities,55 so that we have conducted some calculations to investigate the anisotropy of hole mobility. We find a mobility of 46.8 ± 1.0 cm2 V-1 s-1 in the crystal direction studied in ref. 34, which is perpendicular to the stacking direction r of rubrene, and 6.4 ± 0.5 cm2 V-1 s-1 in the stacking direction. The wavefunction Ψ ( x , t% = 0 ) in these calculations was initially placed inside the system and the applied electric field modeled with a linear potential term in the Hamiltonian. Resolving the directional dependency of the rubrene mobility requires specific preparation of the simulated situation and the choice of a starting wavefunction

that is as

meaningful as possible. More detailed anisotropic organic-crystal mobilities will be presented in future work.

18


Page 19 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


19

Figure 5: The imaginary mobility µimag cm 2V −1s −1  for different source-drain voltages dependent on the gate voltage. (A) The curves obtained from the simulations, where Vg refers to the constant holepotential shift on every grid point of the rubrene crystal Il between the electrodes. (B) The same curves dependent on VSG(VG) (see Figure 3) and multiplied by a constant factor to bring the maxima of the peaks into agreement. (C) Experimental curves taken from reference 34. These data were obtained from a two-point measurement. Figure 3(C) is reprinted by permission of AIP Publishing LLC.

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 30

20

5. SUMMARY As outlined above, imaginary time propagation does not conserve energy. It can thus be a useful tool, not only to saturate to ground states, but also as an approximate energy-loss mechanism for charge transport through structures. It can thus serve as a surrogate for the electron-photon coupling energy-loss mechanism but does not take into account electron-phonon coupling energy loss, which would require a more explicit treatment of the relevant lattice vibrations and their energies.56-58 The question arises, however, as to whether it is possible to find a connection between energy loss per-time

in the

simulations and in experiment. This would make it possible to map the imaginary time scale to real time or to perform a realistic combination of real and imaginary time propagation. Our results and the qualitative agreement with experimental data34 (Figure 5) suggest that this is possible.

6. ACKNOWLEDGMENTS This work was supported by the Deutsche Forschungsgemeinschaft as part of SFB953 “Synthetic Carbon Allotropes” and EXC315 “Engineering of Advanced Materials”. We thank Michael Thoss for fruitful discussions.

20


Page 21 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


21

APPENDIX A: Imaginary time propagation as energy-loss operator Only two simple Gaussian formulae are needed to calculate all integrals: and

In the following, a free particle, with

(30)

, is considered and time independent plane waves

are combined linearly in such a way that the momentum distribution for a single particle is equal to the Boltzmann distribution, where

is the Boltzmann constant. Then, the unitary evolution

operator is applied to this linear combination:

(31) and its complex conjugate:

(32) 21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 30

22

Which leads to the analytic time-dependent solution, where the density has the form of a Gaussian wave packet with a time dependent width:

(33) Which is normalized at all times: (34) Considering the energy:

(35) We also find that the energy is conserved for all times . If we now apply an imaginary time propagation operator to the same wavepacket, we find:

(36) Which leads to the density:

(37) Where the norm is not conserved: (38) 22


Page 23 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


23

To achieve norm conservation, with

at

, the imaginary time propagation operator can

be modified in the following way:

(39) Then one has:

(40) And the new density:

(41) Which then is normalized for all : (42) If we now equate the density created by unitary evolution and that created with imaginary time propagation: (43) Which reads:

(44) Therefore the same density can be created if: (45) or: 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 30

24

(46)

Obviously one can create the same densities with unitary evolution and with imaginary time propagation applied to a linear combination of plane waves, describing a free particle. However, considering the energy:

(47) If we now use the relationship (46):

We see that the energy is not conserved even if the same density is produced as it is with unitary evolution following (46). Therefore, the imaginary time propagation operator can be interpreted as an energy loss or cooling operator if modified to be norm conserving.

REFERENCES Jaeger, C. M.; Schmaltz, T.; Novak, M.; Khassanov, A.; Vorobiev, A.; Hennemann, M.; (1) Krause, A.; Dietrich, H.; Zahn, D.; Hirsch, A.; Halik, M.; Clark, T. Improving the Charge Transport in Self-Assembled Monolayer Field-Effect Transistors: From Theory to Devices. J. Am. Chem. Soc. 2013, 135, 4893-4900. (2)

Organic Electronics II: More Materials and Applications, Klauk, H., (Ed.), Wiley-VCH, Weinheim, 2012. 24


Page 25 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


25

(3)

Lian, J. X.; Lherbier, A.; Wang, L. J.; Charlier, J.-C.; Beljonne, D.; Oliver, Y. Electronic Structure and Charge Transport in Nanostripped Graphene. J. Phys. Chem.C 2016, 120, 2002420032.

(4)

Jacoboni, C.; Reggiani, L. The Monte Carlo Method for the Solution of Charge Transport in Semiconductors with Applications to Covalent Materials. Rev. Mod. Phys. 1983, 55, 645-705.

(5)

Principles of the Theory of Solids, Ziman, J. M., Cambridge University Press, Cambridge, 1972.

(6)

Marcus, R. A. On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. II. Applications to Data on the Rates of Isotopic Exchange Reactions. J. Chem. Phys 1957, 26, 867-871.

(7)

Beeler, J. R. Displacement Spikes in Cubic Metals. I. α-Iron, Copper, and Tungsten. Phys. Rev. 1966, 150, 470-487.

(8)

Brown, P. J.; Sirringhaus, H.; Harrison, M.; Shkunov, M.; Friend, R. H. Optical Spectroscopy of Field-Induced Charge in Self-Organized High Mobility Poly(3-hexylthiophene). Phys. Rev. B 2001, 63, 125204/125201-125204/125211.

(9)

D’Avino, G.; Olivier, Y.; Muccioli, L.; Beljonne, D. Do charges delocalize over multiple molecules in fullerene derivatives? J.Mater.Chem.C 2016, 4, 3747

(10)

Shubina, T. E.; Sharapa, D. I.; Schubert, C.; Zahn, D.; Halik, M.; Keller, P. A.; Pyne, S. G.; Jennepalli, S.; Guldi D. M.; Clark, T. Fullerene van der Waals Oligomers as Electron Traps, J. Am. Chem. Soc. 2014, 136, 10890–10893.

(11)

Coropceanu, V.; Cornil, J.; Filho, D. A. D. S.; Olivier, Y.; Silbey, R.; Bredas, J. L. Charge Transport in Organic Semiconductors. Chem. Rev. 2007, 107, 926-952.

(12)

Li, X.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab Initio Ehrenfest Dynamics. J. Chem. Phys 2005, 123, 084106.

(13)

Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys 1990, 93, 10611071.

(14)

Tully, J. C. Perspective: Nonadiabatic Dynamics Theory. J. Chem. Phys 2012, 137, 22A301.

(15)

Tully, J. C; Berne, B. J.; Ciccotti, G.; Coker, D. F.; (Ed.), Classical and Quantum Dynamics in Condensed Phase Simulations. World Scientific, Singapore, 1998, 34-71.

(16)

Do, V. N. Non-equilibrium Green function method: Theory and Application in Simulation of Nanometer Electronic Devices Adv. Nat. Sci. 2014, 033001. 25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 30

26

(17) Croy, A.; Saalmann, U. Propagation scheme for nonequilibrium dynamics of electron transport in nanoscale devices Phys. Rev. B 2009, 80, 245311 (18) Leitherer, S.; Jäger, C. M.; Krause, A.; Halik, M.; Clark, T.; Thoss, M. Simulation of Charge Transport in Organic Semiconductors: A Time-Dependent Multiscale Method Based on NonEquilibrium Green’s Functions 2017, arXiv:1705.07323v1 (19)

Emberly, E. G.; Kirczenow, G. Landauer Theory, Inelastic Scattering, and Electron Transport in Molecular Wires Phys. Rev. B 1999, 61, 5740.

(20)

Leitherer, S.; Jaeger, C. M.; Halik, M.; Clark, T.; Thoss, M. Modeling Charge Transport in C60-Based Self-Assembled Monolayers for Applications in Field-Effect Transistors J. Chem. Phys. 2014, 140, 204702.

(21)

Jurchescu, O. D.; Baas, J.; Palstra, T. T. M. Effect of Impurities on the Mobility of Single Crystal Pentacene. Appl. Phys. Let. 2004, 84, 3061-3063.

(22)

Karl, N. Charge Carrier Transport in Organic Semiconductors. Synth. Met. 2003, 133-134, 649-657.

(23)

Wang, L.; Prezhdo, O., V.; Beljonne, D. Mixed quantum-classical dynamics for charge transport in organics. Phys. Chem. Chem. Phys. 2015, 17, 12395.

(24)

Dong, H.; Jiang, L.; Hu, W. Interface Engineering for High-Performance Organic Field-Effect Transistors. Phys. Chem. Chem. Phys. 2012, 14, 14165-14180.

(25)

Hwang, J.; Wan, A.; Kahn, A. Energetics of Metal–Organic Interfaces: New Experiments and Assessment of the Field. Mat. Sci. Eng.: R: Rep. 2009, 64, 1-31.

(26)

Goldmann, C.; Haas, S.; Krellner, C.; Pernstich, K. P.; Gundlach, D. J.; Batlogg, B. Hole Mobility in Organic Single Crystals Measured by a “Flip-Crystal” Field-Effect Technique. J. Appl. Phys. 2004, 96, 2080-2086.

(27)

Mladenovic, M.; Vukmirovic, N.; Stankovic, I. Electronic States at Low-Angle Grain Boundaries in Polycrystalline Naphthalene. J. Phys. Chem. C 2013, 117, 15741-15748.

(28)

Bisri, S. Z.; Takenobu, T.; Takahashi, T.; Iwasa, Y. Electron Transport in Rubrene SingleCrystal Transistors. Appl. Phys. Let. 2010, 96, 183304.

(29)

Etschel, S. H.; Waterloo, A. R.; Margraf, J. T.; Amin, A. Y.; Hampel, F.; Jaeger, C. M.; Clark, T.; Halik, M.; Tykwinski, R. R. An Unsymmetrical Pentacene Derivative with Ambipolar Behaviour in Organic Thin-Film Transistors. Chem. Commun. 2013, 6725-6727.

(30)

Boer, R. W. I. D.; Gershenson, M. E.; Morpurgo, A. F.; Podzorov, V. Organic Single-Crystal Field-Effect Transistors. Phys. Stat. Sol. A 2004, 201, 1302-1331. 26


Page 27 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


27

(31)

Reese, C.; Bao, Z. Organic Single-Crystal Field-Effect Transistors. Mat. Today 2007, 10, 2027.

(32)

Hasegawa, T.; Takeya, J. Organic Field-Effect Transistors Using Single Crystals. Sci. Tech. Adv. Mater. 2009, 10, 024314.

(33)

Jurchescu, O. D.; Meetsma, A.; Palstra, T. T. M. Low Temperature Structure of Rubrene Single Crystals Grown by Vapour Transport. Acta Cryst. B 2006, B62, 330-334.

(34)

Podzorov, V.; Sysoev, S. E.; Loginova, E.; Pudalov, V. M.; Gershenson, M. E. Single-Crystal Organic Field Effect Transistors with the Hole Mobility ~8 cm2/Vs. Appl. Phys. Let. 2003, 83, 3504-3506.

(35)

Ehresmann, B.; Martin, B.; Horn, A. H. C.; Clark, T. Local Molecular Properties and their Use in Predicting Reactivity. J. Mol. Model. 2003, 9, 342-347.

(36)

Clark, T. The Local Electron Affinity for Non-Minimal Basis Sets. J. Mol. Model. 2010, 16, 1231-1238.

(37)

Sjoberg, P.; Murray, J. S.; Brinck, T.; Politzer, P. Average Local Ionization Energies on the Molecular Surfaces of Aromatic Systems as Guides to Chemical Reactivity. Can. J. Chem. 1990, 68, 1440-1443.

(38)

Brinck, T.; Carlqvist, P.; Stenlidhem, J. H. Local Electron Attachment Energy and its Use for Predicting Nucleophilic Reactions and Halogen Bonding. J. Phys. Chem. A 2016, 120, 1002310032.

(39)

Bauer, T.; Jaeger, C. M.; Jordan, M. J. T.; Clark, T. A Multi-Agent Quantum Monte Carlo Model for Charge Transport: Application to Organic Field-Effect Transistors. J. Chem. Phys 2015, 143, 044114.

(40)

Hennemann, M.; Clark, T. EMPIRE: A Highly Parallel Semiempirical Molecular Orbital Program: 1: Self-Consistent Field Calculations. J. Mol. Model. 2014, 20, 2331.

(41)

Margraf, J. T.; Hennemann, M.; Meyer, B.; Clark, T. EMPIRE: A Highly Parallel Semiempirical Molecular Orbital Program: 2: Periodic Boundary Conditions. J. Mol. Model. 2015, 21, 144.

(42)

Wick, C. R.; Hennemann, M.; Stewart, J. J. P.; Clark, T. Self-Consistent Field Convergence for Proteins: A Comparison of Full and Localized-Molecular-Orbital Schemes, J. Mol. Model., 2014, 20, 2159

(43)

Atienza, C.; Martin, N.; Wielepolski, M.; Haworth, N.; Clark, T.; Guldi, D. Tuning Electron Transfer Through p-Phenyleneethynylene Molecular Wires. Chem. Commun. 2006, 3202-3204. 27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 30

28

(44)

Davies, K. T. R.; Flocard, H.; Krieger, S.; Weiss, M. S. Application of the Imaginary Time Step Method to the Solution of the Static Hartree-Fock Problem. Nucl. Phys. A 1980, 342, Nr. 1.

(45)

Truong, T. N.; Tanner, J. J.; Bala, P.; McCammon, J. A.; Lesyng, B.; Hoffman, D. K. A Comparative Study of Time Dependent Quantum Mechanical Wave Packet Evolution Methods. J. Chem. Phys 1992, 96, 2077-2084.

(46)

Wick, G. C. Properties of Bethe-Salpeter Wave Functions. Phys. Rev. 1954, 96, 1124–1134.

(47)

Favella, L. F. Brownian Motions and Quantum Mechanics. Annales de l'I.H.P. Physique théorique 1967, 7, 77-94.

(48)

Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. Development and Use of Quantum Mechanical Molecular Models. 76. AM1: A New General Purpose Quantum Mechanical Molecular Model. J. Am. Chem. Soc. 1985, 107, 3902-3909.

(49)

Pople, J. A.; Santry, D. P.; Segal, G. A. Approximate Self‐Consistent Molecular Orbital Theory. I. Invariant Procedures. J. Chem. Phys 1965, 43, 129-135.

(50)

Pople, J. A.; Beveridge, D. L.; Dobosh, P. A. Approximate Self-Consistent Molecular-Orbital Theory. V. Intermediate Neglect of Differential Overlap. J. Chem. Phys 1967, 47, 2026-2033.

(51)

Clark, T.; Halik, M.; Hennemann, M.; Jäger, C. M. Simulating “Soft” Electronic Devices. Molecular Engineering and Control (M. G. Hicks and C. Kettner, Eds), Logos Verlag 2013, 137-150.

(52)

Clark, T. Simulating Charge Transport in Flexible Systems. Perspectives in Science 2015, 6, 58-65.

(53)

Clark, T.; Stewart, J. J. P. MNDO-Like Semiempirical Molecular Orbital Theory and Its Application to Large Systems. in, Computational Methods for Large Systems, J. J. Reimers(ed.), Wiley, Chichester, 2011, Chapter 8.

(54)

Bourke, P. “Gaussian Cube Files.” 2003 http://paulbourke.net/dataformats/cube/. Accessed 1 Jun 2017.

(55)

Sundar, V., C.; Zaumseil, J.; Podzorov, V.; Menard, E.; Willett, R., L.; Someya, T.; Gershenson M. E.; Rogers J. A. Elastomeric Transistor Stamps: Reversible Probing of Charge Transport in Organic Crystals. Science 2004, 303, 1644-1646.

(56)

Bässler, H. Charge Transport in Disordered Organic Photoconductors. Phys. Stat. Sol. B 1993, 175, 15-56.

(57)

Tamura, H.; Tsukada, M.; Ishii, H.; Kobayashi, N.; Hirose, K. Roles of Intramolecular and Intermolecular Electron-Phonon Coupling on the Formation and Transport of Large Polarons in 28


Online

resource:

Page 29 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


29

Organic Semiconductors. Phys. Rev. B 2012, 86, 035208. (58)

Troisi, A.; Orlandi, G. Charge-Transport Regime of Crystalline Organic Semiconductors. Phys. Rev. Let. 2006, 96, 086601.

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

30

ToC Graphic

Snapshot of a hole (isodensity plot at 7×10-5 a.u.) taken from a simulation of the crystal.

30


Page 30 of 30

Charge Transport in Organic Materials: Norm-Conserving Imaginary

Recommend Documents