Simulation of Temperature-Dependent Charge Transport in Organic

Article pubs.acs.org/JCTC

Simulation of Temperature-Dependent Charge Transport in Organic Semiconductors with Various Degrees of Disorder Alexander Heck,†,‡ Julian J. Kranz,† and Marcus Elstner*,†,‡ †

Department of Chemistry, Karlsruhe Institute of Technology, Kaiserstr. 12, 76131 Karlsruhe, Germany HEiKA - Heidelberg Karlsruhe Research Partnership, Heidelberg University, Karlsruhe Institute of Technology (KIT), Germany

‡

ABSTRACT: Different trends in the temperature dependence of the mobility can be observed in organic semiconductors, which constitutes a serious challenge for theoretical approaches. In this work, we apply an atomistic bottom-up simulation for the calculation of temperaturedependent mobilities of a broad selection of materials, ranging from single crystal to amorphous solid. We evaluate how well the method is able to distinguish temperature dependences of different materials and how the findings relate to experimental observations. The applied method is able to cover the full range of temperature dependencies from activated transport in amorphous materials to band-like transport in crystals. In well-characterized materials, we find good agreement with the experiment and a band-like temperature dependence. In less-ordered materials, we find discrepancies from the experiment that indicated that experimentally studied materials possess a higher degree of disorder than do the simulated defect-free morphologies. ⎧ ⎛ T ⎞2 ⎫ μ(T ) = μ∞ exp⎨−⎜ 0 ⎟ ⎬ ⎩ ⎝T ⎠ ⎭

1. INTRODUCTION The correct reproduction of the temperature dependence of the mobility is a formidable challenge because completely opposite trends can be observed in different materials. In single crystals, the mobility usually decreases with increasing temperature following a power law

μ(T ) = C·T −n

⎪

(2)

The second common rate equation was derived by Marcus,5−7 and it features a more involved temperature dependence with temperature intervals where the rate may also decrease with increasing temperature. When applying Marcus rates to derive mobilities, its temperature dependence therefore strongly depends on the system parameters entering the equation, and no universal trend of μ(T) can be given. However, it was shown that the assumptions made in hopping models are not applicable to organic semiconductors with high mobilities.8,9 Furthermore, the experimental observation of the Hall effect in single crystals of rubrene indicates the band-like motion of the charge carrier.10 It is also still not clear whether charge carriers in high-mobility organic semiconductors are best described in terms of extended Bloch waves or localized on single molecules.11 It seems that there is an original regime of transport in organic semiconductors where dynamic fluctuations between the molecules can transiently localize the charge carrier.12,13 Therefore, it was argued that novel approaches are needed that go beyond the commonly used models.14 Promising approaches are semiclassical simulations where the charge-carrier wave function is directly propagated, while the dynamic motion of the system is explicitly considered at the same time scale. Traditionally, this can be achieved by

(1)

with an exponent typically in the range between 0.5 and 3. This temperature dependence is typical for band transport and originates from the increased scattering of the charge carrier at higher temperatures. On the contrary, in disordered (e.g., solution-processed) materials, the mobility usually increases with temperature,1 which is typical for hopping transport. In these materials, the charge is trapped in a localized state and needs thermal activation for transfer to a neighboring site. The trapping can originate from (1) a disordered morphology, giving rise to a broad distribution of site energies, (2) guest molecules (i.e., impurities) with lower site energies acting as trapping centers, or (3) self-trapping, if the response of the nuclei to the charge causes a sufficiently large relaxation (λ ≫ Hab). Nevertheless, band-like temperature dependence was also found in solution-processed organic semiconductors.2 There are two widely used rate equations that are applied to disordered systems. Assuming a Gaussian disorder of site energies and applying Miller−Abrahams rates3 for the charge transfer (CT), the temperature-dependent mobility can be obtained as4 © 2016 American Chemical Society

⎪

Received: February 26, 2016 Published: May 25, 2016 3087

DOI: 10.1021/acs.jctc.6b00215 J. Chem. Theory Comput. 2016, 12, 3087−3096

Article

Journal of Chemical Theory and Computation Ehrenfest or surface hopping simulations.15 In this work, however, we are interested only in the dynamics of the excess charge carrier, which takes place in a small energy window around the Fermi level, whereas strongly bound electronic states are hardly involved. Therefore, in the special case of CT in organic semiconductors, one has to ask whether the knowledge about all electronic degrees of freedom is actually needed to describe the dynamics of the excess charge. Furthermore, the weak electronic coupling between the molecular building blocks of organic semiconductors allows for a fragmentation of the electronic problem, which is usually not exploited in conventional semiclassical simulations. To reduce the high computational cost, it is beneficial to apply model Hamiltonians, which describe the system in terms of charged molecular sites, in Ehrenfest16,17 or surface hopping18,19 simulations. A drawback of these model Hamiltonians, however, is that one has to derive model parameters in advance, which may be straightforward for periodic (crystalline) systems but gets less clear-cut for disordered (amorphous) systems. Furthermore, simplifications are introduced such as considering only one (or a few) normal mode(s) for each site and describing the (nonlocal) electron−phonon coupling within a linear approximation. The latter approximation can become especially problematic, considering the large intermolecular motions and the high sensitivity of the electronic coupling on the relative molecular orientation. The method that we apply in this work can be seen as a kind of Ehrenfest simulation that falls in between the conventional full ab initio approach and the application of model Hamiltonians. On the one hand, we remain at an atomistic resolution and get all site energies, couplings, and forces on the fly. On the other hand, we propagate only the excess charge and describe the system within the framework of molecular sites. This method was applied in a similar way to CT in DNA20−22 and proteins.23−25 Recently, we applied it to CT in selfassembled monolayers (SAMs)26 and made further advancements for the simulation of organic semiconductors,27 where we found that the eigenstates of the charge carrier are localized on a few molecules because of thermal disorder and no stable polaron was formed. This article is organized as follows. In section 2, we give a brief overview of the applied simulation scheme and the calculation of mobilities. Afterwards, we investigate a broad range of organic semiconductor systems, which differ largely by the inherent degree of molecular order. The respective molecular building blocks are collectively shown in Figure 1. As a typical ordered system, we investigate in section 4.1 a single crystal of anthracene, which can be obtained at very high purity28 and for which robust data of the intrinsic mobility at different temperatures are available.29 As a typical disordered material, we study in section 4.2 the amorphous phase of N,N′bis(1-naphthyl)-N,N′-diphenyl-1,1′-biphenyl-4,4′-diamine (αNPD), which is widely used in organic light-emitting diodes as a hole-injection layer,30 hole-transport layer,31,32 or electronblocking layer.33 For these two materials, completely different temperature dependences of the mobility are observed, which we try to capture with our simulation method. Subsequently, we extend our simulations to polymers and liquid crystals, whose morphologies are in between completely ordered and disordered. As a representative of a polymer, we investigate in section 4.3 the interstrand transport in an ordered section of regioregular poly(3-hexylthiophen) (P3HT), which is prominent as a donor material in organic solar cells.34 As a

Figure 1. Molecular building blocks of the systems studied in this work: (a) anthracene, (b) α-NPD, (c) P3HT, (d) HBC-LC, and (e) HBC-SAM.

liquid crystal, we investigate in section 4.4 a hexabenzocoronen (HBC) derivate with racemic branched (3,7-dimethyloctanyl) side chains (HBC-LC).35 In the liquid crystalline phase of HBC, the molecules stack in an unfavorable conformation, in which the intermolecular electronic coupling is near its minimum.36 By modifying the side chains, the relative orientation of molecules in columnar liquid crystals can be changed, which can result in higher mobilities.37 Nevertheless, one drawback remains in these liquid crystals: the hardly known degree of impurities and defects in solvent-processed HBC derivatives, combined with their strong impact on the chargecarrier mobility, complicates the comparison between the experiment and the simulation. In combined experimental and theoretical studies,36 it was found that the mobility from pulseradiolysis time-resolved microwave conductivity (PR-TRMC) measurements and hopping simulations of a disordered columnar liquid crystal agree quite well, but as soon as the disorder is removed in the simulation by proper annealing, the mobility increases by nearly 2 orders of magnitude, indicating that structural defects may play an important role in these materials. Interesting systems that have the potential to reduce the amount of defects and allow the monitoring of intrinsic mobilities are SAMs of π-stacking organic molecules. After functionalizing the side chains of HBC derivatives with thiol groups, it is possible to obtain very ordered stacks, where the absence of defects can be probed with scanning tunneling microscope (STM) images.26,38 Band-like behavior was proposed for such a SAM,38 and room temperature mobilities of up to 6.8 cm2/(V·s) have been estimated experimentally.26 In this work, we will take a closer look at the temperature dependence of such a SAM (see HBC-SAM in Figure 1) and compare its performance to that of HBC-LC.

2. METHODS 2.1. Coupled Electron-Ion Dynamics Simulations. Here, we give a brief summary of the simulation method, but we direct the reader to ref 27 for an extensive discussion of this approach. The idea behind the coupled electron-ion dynamics (CEID) simulation scheme is that we do not need to propagate 3088


Article

Journal of Chemical Theory and Computation all electronic degrees of freedom, in contrast to conventional Ehrenfest simulations, because we are interested only in the dynamics of the excess charge carrier. We therefore divide the system into a quantum mechanical (QM) region, where the charge carrier is located, and the remainder of the system, which is described by molecular mechanics (MM). Furthermore, the charge transport takes place in a small energy window around the Fermi level, whereas strongly bound electronic states are hardly affected. We therefore approximate the energy of the cationic system with tot E + ≈ EMM − ⟨Ψ|H[ρ0 ]|Ψ⟩ + ΔEQM/MM

aṁ = −i∑ Hmnan − n

∑∑

am|ϕm⟩

(3)

ρ0 =

∑

mK R̈ k = −

ρA0

μ=

D=

⟨Δx 2(t )⟩ nt

K

⎛ Δq mq 0 ⎞ α K ⎟ | R − RK | ⎟⎠ α∈A ⎝ α

∑ ⎜⎜

(11)

∑ (xA(t ) − x0)2 pA (t )

(12)

where xA(t) is the center of mass of molecule A and pA(t) is its occupation given as pA(t) = ∑m∈A|am|2. The center of charge at t = 0 is given as x0 = ∑ ApA(t = 0)xA(t = 0). Note that it is also possible to derive the mobility from the drift that is induced by an electric field. The disadvantage of this approach is discussed in the appendix.

where we calculate the Hamilton matrix elements, H0mn = ⟨ϕm| H[ρ0]|ϕn⟩, with the density functional tight-binding (DFTB) method, because its good performance for the calculation of electronic couplings has been reported.39−41 Note, that we apply here the non-self-consistent variant of DFTB as discussed in ref 27. The QM/MM coupling term is given within the DFTB method as

A m∈A

(10)

A

(6)

|am|2 ∑

(9)

eD kBT

⟨Δx 2(t )⟩ =

∑ ∑ ∑ am*an⟨ϕm|H[ρ0 ]|ϕn⟩ + ΔEQM/MM

∑∑

mn

∂ ⟨ϕ |H[ρ0 ]|ϕn⟩ ∂R k m

where n = 2, 4, and 6 for one-, two- and three-dimensional transport directions, respectively. The mean squared displacement, ⟨Δx2(t)⟩, of the hole-wave function at time t is defined as

AB m ∈ A n ∈ B

ΔEQM/MM =

∑ am*an

where kB is the Boltzmann constant, e denotes the elementary charge, and diffusion constant D is given as

Although it is often sufficient to include only the highest occupied molecular orbital (HOMO) of each fragment in the expansion, it becomes necessary to also include lower-lying orbitals if they are energetically not well separated from the HOMO. The final total energy expression in the FO basis thus reads tot E + ≈ EMM −

(8)

where k denotes x, y, and z of atom K. For an efficient treatment of the QM/MM coupling terms, the instantaneous hole wave function is projected on the partial charges of the force field. 2.2. Derivation of Mobilities. The mobility, μ, of the hole is obtained via the Einstein−Smoluchowski equation:

(5)

A

tot ∂EMM + ∂R k

∂ − ΔE QM/MM ∂R k

(4)

A m∈A

n

where Hmn are the DFTB Hamilton matrix elements including the QM/MM interactions. In a spatial fixed basis, the last term of eq 8 vanishes, whereas in an atom-centered basis, there are some contributions due to the time dependence of the atomic positions.42 In this work, the last term is omitted except in systems where multiple FOs on a single fragment are included. There it allows transitions between these orbitals, which is handled by projecting the wave function on the new basis as described in ref 27. Parallel to the quantum propagation of the hole wave function, a classical propagation of the nuclei is performed. The corresponding equation of motion reads27

where |Ψ⟩ is the single-particle wave function of the electron hole, H[ρ0] is the Kohn−Sham Hamiltonian of the neutral system, Etot MM is the potential energy of the force field of the whole neutral system, including all atoms of the MM and the QM regions, and ΔEQM/MM is the change in the QM/MM interaction with respect to the neutral reference system. Because of the weak intermolecular interactions in molecular semiconductors, it is feasible to partition the system via a fragment orbital (FO) formalism. The total QM zone is divided into single molecules (fragments), for which parallel quantum chemical calculation can be performed. The single-particle wave function, |Ψ⟩, of the electron hole can then be expressed in the basis of (orthogonalized) molecular orbitals, |ϕm⟩, of individual molecule A. |Ψ⟩ =

∑ an⟨ϕm|ϕṅ ⟩

3. COMPUTATIONAL DETAILS Bonded and van der Waals interactions were modeled in all systems with the GAFF43 force-field parameters, if not stated otherwise. Note that the partial charges of the various systems are derived with different methods. However, their values are quite small in the nonpolar molecules that we have studied in this work. Therefore, only a minor impact on the molecular arrangements is expected. Periodic boundary conditions were used with particle-mesh-Ewald (PME)44 electrostatics in equilibration runs and in the CEID simulations. The temperature and pressure equilibration of the neutral systems was performed with the Berendsen scheme45 and lasted for several nanoseconds in all cases. During the CEID runs, the Parrinello−Rahman barostat46,47 and the Nosé−Hoover

(7)

q0K

where is the partial charge on atom K of the classical environment and Δqmα is the change in the Mulliken charge of atom α resulting from the missing electron in orbital m. From the total energy expression, coupled equations of motions for the charge carrier and the nuclei can be derived. By applying the Lagrangian formalism, one obtains the timedependent Schrödinger equation for the electronic degrees of freedom:27 3089


Article

Journal of Chemical Theory and Computation thermostat48,49 were used to yield a correct canonical ensemble at 1 bar and the respective temperature for each system. The CEID method is implemented in a local version of Gromacs 4.6,50 which was also used for the equilibration runs. The global time step of the molecular dynamics (MD) simulation was 1 fs, but note that the Runge−Kutta method used to integrate eq 8 performs a large number of shorter time steps internally to achieve a sufficiently accurate integration. The mobility was calculated according to eq 10. The mean square displacement, ⟨Δx2(t)⟩, was taken as an average over 100 independent simulations to achieve sufficient statistical sampling. The starting wave function in each run was the lowest eigenvector of the FO Hamiltonian. The feasible duration of a single run was limited by the time it takes for the excess charge to reach the edge of the quantum region and therefore depends on the mobility of the material and the size of the QM zone. Therefore, the simulations normally lasted for 500 fs, but because of the low mobility of α-NPD, the individual charge transport runs could be performed for 10 ps. The mobility was then obtained from a fit of ⟨Δx2(t)⟩ versus time in the linear region before excessive boundary effects set in. If not stated otherwise, only the HOMO of each fragment was used to express the hole wave function (eq 4). 3.1. Anthracene. The anthracene single crystal was taken from ref 27, where atomic partial charges were fitted according to the RESP protocol.51 The initial coordinates of a unit cell were obtained from the crystal structure,52 and the complete system was modeled as a supercell consisting of 20 × 20 × 20 unit cells. For each crystallographic direction, a QM zone was constructed from 18 sequential molecules along the respective axis. Simulations were performed for temperatures in the range of 150−400 K.

Figure 3. QM zone of α-NPD (green). From the 300 molecules in the box, the QM zone was constructed as a cluster of 27 molecules around the molecule with the lowest site energy.

energy (lowest site energy) by including the closest 26 neighbors. 3.3. P3HT. The simulation box of P3HT consisted of 400 parallel strands, each built of 20 monomers. Force-field parameters were taken from ref 53, where OPLS-AA54,55 force-field parameters were used with partial charges from ref 56 and dihedral parameters were fitted to a B3LYP/6311+g(d,p) potential energy scan. The initial structure corresponds to morphology I′, as described in ref 53. A single fragment in the CEID calculations was formed of 10 sequential monomers on a single strand, and for each fragment, the 10 highest occupied orbitals were taken to form the basis of the hole wave function. This ensures that the whole length of the fragment may be covered with, potentially localized, states. The remainders of the strands were cut off from the QM zones and were capped by hydrogen atoms. The total quantum region consisted of 20 sequentially π-stacked fragments (see Figure 4). Simulations were performed for temperatures in the range of 150−350 K. 3.4. HBC-LC. The HBC liquid crystal was adopted from ref 27, where partial charges were derived from DFTB calculations and corrected according to the CM3 method.57 The system consists of nine columns of 110 molecules, which are arranged in a hexagonal lattice. The molecules within a column stack

Figure 2. Three QM zones of anthracene along the a (red), b (green), and c (blue) crystallographic axes. A total of 18 sequential molecules from the 8000 molecules in the box were selected.

3.2. α-NPD. The force-field partial charges of α-NPD were obtained from PM6 calculations. For the generation of the amorphous starting structures, an annealing procedure was performed, where a simulation box of 300 molecules was first melted at 700 K and then successively cooled down. The cooling steps consisted of a temperature reduction of 50 K over 1 ns, followed by a 1 ns simulation at constant temperature, from which equidistant starting structures for the CEID simulations were extracted. A subset of 27 molecules were included in the QM region, as shown in Figure 3. To depict the energetic disorder realistically and capture the rate-limiting detrapping steps, the HOMO energies of all 300 molecules were scanned at the start of each simulation. The QM zone was then constructed around the molecule with the highest HOMO

Figure 4. QM zone (green) of P3HT. The central 10 monomers of 20 strands were selected from the 400 strands in the box. 3090


Article

Journal of Chemical Theory and Computation with a relative twist of 30°. The simulations were performed in the temperature range from 400 to 600 K because a phase transition from the liquid crystalline to the crystalline phase occurs below 350 K.35 As fragments in the quantum calculations, the aromatic HBC cores were used, with side chains being replaced by hydrogen capping atoms. The QM zone consisted of 30 sequential molecules on the same column, and the two degenerate HOMOs of HBC were considered in the CEID simulation. 3.5. HBC-SAM. A SAM of HBC-SAM was modeled as 110 molecules in a single stack. The last atom of both anchoring side chains was kept fixed during the simulation to mimic the anchoring to the surface. Like in HBC-LC, a random configuration at the stereogenic centers of the side chains was used. The gold surface was not modeled explicitly because in such an upright orientation the molecules stand on their anchoring side chains, leaning on each other, and therefore no significant interaction with the surface can be expected. The angle between the molecular plane and the surface correlates with the degree of surface coverage. Several degrees of molecular coverage were tested, and the experimental NEXAFS angle of 65° was reproduced best with an intermolecular spacing of 3.9 Å along the stacking direction.26 No pressure coupling was applied, with fixed box dimensions in the surface plane of 42.9 and 3 nm along and perpendicular to the stacking direction, respectively. Three-dimensional periodic boundary conditions were applied, enabling PME electrostatics. However, the box axis perpendicular to the surface was increased to 6 nm to avoid steric interactions, whereas force and potential corrections were applied in the PME calculations to produce a pseudo-2D summation. The QM zone consisted of 30 sequential molecules on the same column. Like for HBC-LC, the fragments consisted of the HBC cores and the two degenerate HOMOs of each molecule were used in the CEID simulation.

Figure 7. Exemplary time evolution of the excess charge in anthracene. The initially localized charge oscillates between neighboring sites and delocalizes quickly until it is spread over the complete system.

to broaden. Fast oscillations between sites with similar energy are observed until the charge spreads over the whole system. As can be seen in Figure 8, the simulated hole mobilities range from 0.3 to 5.8 cm2/(V·s) and decrease with temperature.

Figure 8. Hole mobilities in anthracene at different temperatures in the direction of the three crystallographic axes as shown in Figure 2 and fits of μ = C·T−n to the simulated data.

The temperature dependence is well described by a power law, and the least-square fitting of μ = C·T−n to the simulated data yields exponent n between 1 and 1.5, depending on the crystallographic axis with coefficient of determination for the fits of r2 > 0.988. Such a thermal dependence of the mobility is characteristic of large-polaron, band-like behavior. Because anthracene is an ultrapure organic crystal, this is the expected result. The simulated mobilities are in very good agreement with the experiment, as can be seen in Table 1. Compared to timeof-flight (TOF) measurements on highly purified single crystals, the simulated mobilities deviate by less than a factor of 2. The same holds true for the exponent in the power law, which describes the temperature dependence of the mobility. The slightly underestimated mobility in anthracene could be caused by the artificial restriction of the QM zone to one dimension, whereas in an actual three-dimensional system, the charge carrier could circumvent bottlenecks, introduced by thermal fluctuations in the CT parameters. 4.2. Study of a Disordered Material: α-NPD. As can be seen in Figure 9a, a completely different temperature dependence compared to the previous system is found. The mobility increases with temperature, which indicates activated hopping, and is with 2.3−4.5 × 10−2 cm2/(V·s) about 2 orders

Figure 5. QM zone (green) of HBC. From the nine columns in the box, the QM zone was constructed as 30 sequential molecules in a single column.

Figure 6. Thirty molecules in the QM zone (green) of HBC-SAM. The monolayer was modeled as a single column of 110 molecules with an increased box dimension perpendicular to the virtual gold surface.

4. RESULTS AND DISCUSSION 4.1. Study of an Ordered Material: Anthracene. An exemplary trajectory showing the charge dynamic of a single run of anthracene is shown in Figure 7. Because of thermal vibrations, the wave function is still quite localized at the beginning of the simulation, even in this crystalline system. No stable polaron is formed, and the charge starts instantaneously 3091


Article

Journal of Chemical Theory and Computation

preparation conditions but by applying admittance spectroscopy,61 whereas thin-film transistor (TFT) measurements yielded mobilities 1 order of magnitude lower.60 Therefore, compared to experimental values, our calculated room temperature mobility is about 2 orders of magnitude too high. However, the room temperature mobility (μRT) is strongly affected by impurities and defects inside the material and is therefore not a very good estimate for the intrinsic transport properties of an organic material. For amorphous materials such as α-NPD, the purification is more difficult than for crystalline materials, and the reduction of μRT due to defects and impurities can be expected to be a dominant factor. It was argued that the hightemperature limit of the mobility (μ∞) is better suited for describing the intrinsic mobility because at high temperatures the activated detrapping of the excess charge ceases to be the limiting factor.62 This was illustrated for α-NPD, doped with various compounds.62 Nearly the same μ∞ was found for all doped systems, whereas their μRT differed by more than 1 order of magnitude. Our simulated high-temperature limit of the mobility, μ∞ = 4.8 × 10−2, is very similar to the experimental value of μ∞ = 1.7 × 10−2 cm2/(V·s). The weaker decrease of μ at lower temperature results thus from a weaker trapping (smaller σ) in our simulations compared to the experiment. One factor may be the presence of impurities in the experiments, but furthermore, it is surprising that the σ that we obtain from the fit is quite small (18 meV). Note that the actual standard deviation of the site energies, obtained from snapshots along our MD simulations, is ∼150 meV, in good agreement with experimental values.63 The discrepancy between actual σ of site energies and the apparent σ resulting from the fit points to a problem of the Ehrenfest simulation scheme. In contrast to surface hopping simulations, there is no decoherence that allows the charge-carrier wave function to fully relax to the lowest adiabatic state. Therefore, when the energetic disorder is derived from the dynamics of the charge, the depth of the traps appears to be smaller than it actually is. An exemplary trajectory, representative of the charge dynamic in α-NPD, is shown in Figure 10. The charge is for about 2.5 ps localized on the initial site. Then, it hops to a neighboring site, where it stays for 1 ps before successive hops occur. Additionally, one can observe a continuous broadening

Table 1. Comparison of the Simulated Hole Mobilities of Anthracene with Experimental Values Obtained from TOF Measurementsa a-direction 200 K 250 K 300 K n

b-direction

c-directionb

sim.

exp.

sim.

exp.

sim.

exp.

1.37 1.14 0.87 1.09

2.22 1.51 1.14 1.50

3.95 3.09 2.01 1.45

5.05 3.74 2.93 1.34

1.41 0.96 0.70 1.46

1.98 1.20 0.85 2.38

a

Mobilities are shown only for three selected temperatures at which the mobility tensor was given in ref 29 together with exponent n after the fit of μ = C·T−n to all simulated and experimental data. b Orthogonalized coordinate system a⊥b⊥c′ was used in the experiment, where c′ differs by about 35° from the crystallographic c-direction used in our simulations.

Figure 9. Simulated mobilities of α-NPD for different temperatures in the range of 150−400 K. (a) Absolute value of the mobility for different temperatures. (b) Fit of μ = μ∞·exp[−(2σ/3kT)2] to the simulated mobilities of α-NPD.

of magnitude smaller. Such a temperature dependence can be expected considering the strong disorder in the amorphous material. The temperature dependence can be explained by the Gaussian disorder model and is fitted quite well by μ = μ∞· exp[−(2σ/3kT)2] with r2 > 0.892 (see Figure 9b). The fit yields a high-temperature limit of the mobility of μ∞ = 4.8 × 10−2 cm2/(V·s) and a standard deviation of site energies of 18 meV. TOF measurements at room temperature on vacuum-vapordeposited thin films yield hole mobilities of 3−4 × 10−4 cm2/ (V·s),58−60 and similar results are also obtained under the same

Figure 10. Exemplary time evolution of the excess charge in α-NPD. In contrast to Figure 7, there is a significant residing time between intermolecular transfers. Note that because of the amorphous threedimensional QM zone of α-NPD the molecular indices shown on the x axis do not reflect the molecular positions, in contrast to the case in Figure 7. 3092


Article


found in webs of these fibers, in contrast to the band-like temperature dependence of our simulations. However, reducing the polarity of the SiO 2 surface by treatment with hexamethyldisilazane reduces the activation barrier from 108 to 65 meV,71 which indicates a strong influence of the surface on the intrinsic mobility. Furthermore, the polymer in our simulation is 100% regioregular, which is not achievable in the experiment. The best commercially (Sigma-Aldrich, October 2015) obtainable P3HT has a regioregularity of >98%, but values around 95% are also quite common. 4.4. Enforcing Order: HBC. In this section, we simulate the temperature-dependent mobility of HBC-LC, which we previously investigated in detail at constant temperature.27 Additionally, we study a SAM derived from HBC,26 where two side chains are functionalized with thiol groups that anchor the molecule to a gold surface and compare how the altered morphology influences the charge transport characteristics. As can be seen in Figure 12, we find decreasing mobility with increasing temperature for both systems. In HBC-LC, the

of the charge-carrier wave function, which is caused by the missing decoherence in Ehrenfest simulations, as discussed earlier. 4.3. Study of a Semiordered Material: P3HT. Many conjugated polymers exhibit complex morphologies, where crystalline domains are embedded in an amorphous matrix.64 The simulation of such an involved structure is, however, beyond the scope of this work. Therefore, we focus on a single aspect of CT in polymers, namely the interchain transport in crystalline domains. As can be seen from Figure 11, we find mobilities in the range of 1.5−4 cm2/(V·s), which decrease with increasing

Figure 11. Simulated mobility of P3HT at different temperatures, which can be fitted by a power law.

temperature. Again, the mobility follows a power law (μ ∼ T−n) with an exponent of n = 1.05, indicating large-polaron, bandlike behavior. The coefficient of determination of the fit is r2 = 0.981. In carefully prepared field-effect transistors (FETs), mobilities as high as 0.1 cm2/(V·s) are obtainable,65 which is 1 order of magnitude smaller than our simulated mobility. Direct comparison of these values to the experiment is difficult, however, because we face again uncertainties regarding the purity of the material and its influence on the mobility, as discussed already in section 4.2. This difficulty is illustrated by the fact that the mobility of regioregular P3HT can change by 3 orders of magnitude, depending on the applied solvent and the preparation conditions in FET measurements.66 Additionally, the morphology of P3HT is dependent on its molecular weight and can change after annealing.67 We simulate a single crystalline domain, whereas standard experiments such as TOF or PR-TRMC probe polycrystalline materials with grain boundaries between the crystallites, which reduce the observed direct current mobility. The formation of grain boundaries appears to be strongly dependent on the molecular weight of the polymer, with low-molecular-weight P3HT layers being more crystalline than their high-molecular-weight counterparts.68 It was further proposed that long polymer chains might bridge crystallites in high-molecular-weight P3HT, leading to higher mobilities.69 One possibility to get rid of grain boundaries and measure the pure interchain transport in crystalline domains of P3HT is the study of nanofibers. In these fibers, the P3HT chains align perpendicular to the fiber direction and stack on top of each other.70 Contacting opposing ends of a fiber thus allows the measurement of interchain transport along the fiber. Fourpoint-FET measurements yield a mobility of 0.06 cm2/(V·s) for the transport along a single fiber.71 Activated transport was

Figure 12. Simulated mobility of the liquid crystalline derivative of HBC (HBC-LC) and of the self-assembling monolayer of HBC (HBC-SAM) for different temperatures fitted by a power law (μ = C· T−n).

mobility ranges from 1.81 cm2/(V·s) at 400 K to 1.19 cm2/(V· s) at 600 K and follows a power law (μ ∝ T−n) with an exponent of n = 1.06. Band-like behavior is less expected in the HBC liquid crystal compared to the anthracene single crystal because of the increased disorder. However, the good fit (r2 > 0.965) of μ = C·T−n to the simulated data and the exponent of n = 1.06, which is in the typical range for band-like transport, strongly indicate that the large-polaron, band-like transport is the predominant charge transport regime for the defect-free liquid crystal. Because of the different morphology, significantly higher mobility is observed in HBC-SAM ranging from 3.99 to 4.93 cm2/(V·s). At 400 K, the mobility is larger by a factor of 2, compared to HBC-LC, and also the temperature dependence is significantly weaker with an exponent of n = 0.45. Compared to PR-TRMC measurements on the liquidcrystalline phase of HBC derivatives, we overestimate the mobility roughly by a factor of 5.72 Furthermore, the experimental mobility is observed to be nearly temperature independent in this phase. As already mentioned in Section 4.2, missing impurities and defects would explain an overestimation of the mobility. Furthermore, because the detrapping of the charge is temperature-activated, the inclusion of defects in the simulation would counteract the trend of decreasing mobility at higher temperatures, yielding a temperature dependence in 3093


Article


experimental observations, we found decreasing mobility with increasing temperature also in semicrystalline materials such as the liquid crystal HBC-LC and the polymer P3HT. From the extrapolation of the mobility to the high-temperature limit in αNPD, we concluded that we are underestimating the depth of the charge trapping. This can result either from a too weak molecular relaxation or from the absence of impurities. However, there are also several drawbacks of the applied method, which were discussed in detail in ref 27. Without any restrictions, the propagation of the charge carrier leads to occupation of energetically higher adiabatic states, in violation with a Boltzmann distribution. This overheating of the charge carrier is a common problem, which may be corrected by additional terms in the propagation.76 The detrapping of the charge carrier in α-NPD is affected probably the most by this shortcoming, whereas the effect on anthracene is expected to be much smaller. Furthermore, we confirmed that the preparation of a SAM from HBC derivatives can lead to higher mobilities compared to their liquid crystalline counterpart, and our simulations showed that the transport is also less temperature sensitive, which is a desirable feature for high-performance devices. Atomistic bottom-up simulations of charge-carrier dynamics with such a multiscale method are thus a valuable tool for complementing experimental results and establishing structure−property relations and might be used for the rational design of new compounds in future applications.

better agreement with the experiment. Another observation, strengthening the assumption that defects may be responsible for the differences in the mobility between experiment and simulation, is that even in the crystalline phase of various HBC derivatives an increasing mobility with higher temperatures was found.72,73 It was argued that such a temperature dependence in crystalline materials points to insufficient purification and does not reflect the intrinsic mobility of the material.14 Furthermore, it is known that rotations of the molecular disks around the columnar axis take place in the liquid crystalline phase.74 Considering the strong dependence of the electronic coupling on the relative orientation of two molecules,36 such structural defects can have a strong impact on the charge-carrier dynamics of the system. In our MD simulations, however, these slow structural changes are not observable on the nanosecond time scale. To get better insight into these kind of defects, more elaborate methods such as umbrella sampling75 have to be applied. The calculation of relative activation barriers for the rotation of a single molecule versus the concerted rotation of a larger section of the columnar stack could be a first step in the search for the most prominent types of defects. For the single molecular layer in HBC-SAM, it is not possible to derive mobilities with standard techniques such as TOF, PRTRMC, or TFT measurements. However, it is possible to derive some estimates from STM data.26,38 For such a system, a mobility of 6.8 cm2/(V·s) was estimated,26 significantly higher than that of the PR-TRMC experiments of the liquid crystal and in good agreement with our calculation. Furthermore, for a similar system with a single anchoring side chain, an increasing mobility at lower temperatures was derived from the STM data,38 which is in agreement with our predicted temperature dependence for pristine HBC-SAMs and liquid crystals.

■

APPENDIX Besides eq 10, mobility μ can also be derived from the drift that is induced by electric field E v = μE

5. CONCLUSIONS We can predict mobilities in a range of 3 orders of magnitude as summarized in Figure 13 and are able to reproduce the

(A.1)

where v is the drift velocity of the center-of-charge. There is no technical problem to derive the mobilities from the fieldinduced drift. Because we are using a finite QM zone without periodic boundary conditions, we can add the electronic potential to each atomic orbital energy and then run the simulations with this altered Hamiltonian. However, in the experiment, it is more practical to measure the mobility by applying electric fields, whereas it is more practical to take the diffusion in the simulation as we outline in the following. In a material with a mobility of 1 cm2/(V·s), exposed to an electric field of 50 kV/cm, the charge carrier travels with 5 Å/ ps, which means only to the nearest neighboring molecule within 1 ps. After the same time, we obtain from the diffusive motion (eq 10) at T = 300 K a root-mean-square displacement of ⟨Δx 2⟩ ≈ 40 Å for a three-dimensional system. This means that the distance of the charge carrier from its original position caused by random diffusion is about 1 order of magnitude larger than that of the directed movement along the electric field. Note that for smaller mobilities, which are very common in organic semiconductors, the diffusive motion has even a larger contribution to the total traveled distance because the

Figure 13. Summary of the temperature-dependent mobilities of all systems studied in this work. System-specific temperature dependences are obtained, and mobilities in a range of 3 orders of magnitude are accessible.

opposing temperature dependencies that are significant for band-like transport in single crystals such as anthracene, as well as activated hopping in amorphous materials such as α-NPD. In well-characterized materials such as anthracene, even a quantitative agreement with the experiment was found, whereas differences between simulation and experiment emerge in materials with more involved morphologies. In contrast to

root-mean-square displacement, ⟨Δx 2⟩ , of the diffusive motion decreases only with μ , whereas the drift distance in an electric field decreases linearly with μ. It is therefore more practical to derive μ directly from the diffusive motion, which would otherwise appear as a strong noise in the simulations of directed transport in an electric field. 3094


Article


■

(32) Son, K. S.; Yahiro, M.; Imai, T.; Yoshizaki, H.; Adachi, C. Chem. Mater. 2008, 20, 4439−4446. (33) Lee, S.; Chung, C.-H.; Cho, S. M. Synth. Met. 2002, 126, 269− 273. (34) Dang, M. T.; Hirsch, L.; Wantz, G. Adv. Mater. 2011, 23, 3597− 3602. (35) Fechtenkötter, A.; Tchebotareva, N.; Watson, M.; Müllen, K. Tetrahedron 2001, 57, 3769−3783. (36) Feng, X.; Marcon, V.; Pisula, W.; Hansen, M. R.; Kirkpatrick, J.; Grozema, F.; Andrienko, D.; Kremer, K.; Müllen, K. Nat. Mater. 2009, 8, 421−426. (37) Pisula, W.; Feng, X.; Müllen, K. Adv. Mater. 2010, 22, 3634− 3649. (38) Käfer, D.; Bashir, A.; Dou, X.; Witte, G.; Müllen, K.; Wöll, C. Adv. Mater. 2010, 22, 384−388. (39) Kubas, A.; Gajdos, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Phys. Chem. Chem. Phys. 2015, 17, 14342−14354. (40) Kubas, A.; Hoffmann, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. J. Chem. Phys. 2014, 140, 104105. (41) Kubas, A.; Hoffmann, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. J. Chem. Phys. 2015, 142, 129905. (42) Deumens, E.; Diz, A.; Longo, R.; Ö hrn, Y. Rev. Mod. Phys. 1994, 66, 917−983. (43) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174. (44) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089−10092. (45) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690. (46) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182−7190. (47) Nosé, S.; Klein, M. L. Mol. Phys. 1983, 50, 1055−1076. (48) Nosé, S. Mol. Phys. 1984, 52, 255−268. (49) Hoover, W. G. Phys. Rev. A 1985, 31, 1695−1697. (50) Páll, S.; Abraham, M.; Kutzner, C.; Hess, B.; Lindahl, E. Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS. In Solving Software Challenges for Exascale; Markidis, S., Laure, E., Eds.; Springer International Publishing: Switzerland, 2015; Vol. 8759, pp 3−27. (51) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269−10280. (52) Brock, C. P.; Dunitz, J. D. Acta Crystallogr., Sect. B: Struct. Sci. 1990, 46, 795−806. (53) Poelking, C.; Andrienko, D. Macromolecules 2013, 46, 8941− 8956. (54) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657−1666. (55) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225−11236. (56) Moreno, M.; Casalegno, M.; Raos, G.; Meille, S. V.; Po, R. J. Phys. Chem. B 2010, 114, 1591−1602. (57) Kalinowski, J. A.; Lesyng, B.; Thompson, J. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A 2004, 108, 2545−2549. (58) Naka, S.; Okada, H.; Onnagawa, H.; Yamaguchi, Y.; Tsutsui, T. Synth. Met. 2000, 111−112, 331−333. (59) Lin, C.; Chen, Y.; Chen, H.; Fang, F.; Lin, Y.; Hung, W.; Wong, K.; Kwong, R.; Xia, S. Org. Electron. 2009, 10, 181−188. (60) Cheung, C. H.; Tsung, K. K.; Kwok, K. C.; So, S. K. Appl. Phys. Lett. 2008, 93, 083307. (61) Nguyen, N. D.; Schmeits, M.; Loebl, H. P. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 075307. (62) Tsung, K. K.; So, S. K. Org. Electron. 2009, 10, 661−665. (63) van Mensfoort, S. L. M.; Shabro, V.; de Vries, R. J.; Janssen, R. A. J.; Coehoorn, R. J. Appl. Phys. 2010, 107, 113710. (64) Samuelsen, E. J.; Mardalen, J. Handbook of Organic Conductive Molecules and Polymers; Wiley: Chichester, UK, 1997; Vol. 3, pp 87− 120. (65) Sirringhaus, H.; Tessler, N.; Friend, R. H. Science 1998, 280, 1741−1744.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +49 721 608-45705. Fax: +49 721 608-45710. Author Contributions

A.H. and J.J.K. contributed equally to this study. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Denis Danilov and Carl Poelking for providing forcefield parameters and starting structures of α-NPD and P3HT, respectively.

■

REFERENCES

(1) Sirringhaus, H. Adv. Mater. 2005, 17, 2411−2425. (2) Sakanoue, T.; Sirringhaus, H. Nat. Mater. 2010, 9, 736−740. (3) Miller, A.; Abrahams, E. Phys. Rev. 1960, 120, 745−755. (4) Bässler, H. Phys. Status Solidi B 1993, 175, 15−56. (5) Marcus, R. A. J. Chem. Phys. 1956, 24, 966−978. (6) Marcus, R. A. J. Chem. Phys. 1956, 24, 979−989. (7) Marcus, R. A. Annu. Rev. Phys. Chem. 1964, 15, 155−196. (8) Gajdos, F.; Oberhofer, H.; Dupuis, M.; Blumberger, J. J. Phys. Chem. Lett. 2013, 4, 1012−1017. (9) Troisi, A. Org. Electron. 2011, 12, 1988−1991. (10) Podzorov, V.; Menard, E.; Rogers, J. A.; Gershenson, M. E. Phys. Rev. Lett. 2005, 95, 226601. (11) Sirringhaus, H.; Sakanoue, T.; Chang, J.-F. Phys. Status Solidi B 2012, 249, 1655−1676. (12) Fratini, S.; Mayou, D.; Ciuchi, S. Adv. Funct. Mater. 2016, 26, 2292−2315. (13) Ciuchi, S.; Fratini, S.; Mayou, D. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, No. 081202, DOI: 10.1103/PhysRevB.83.081202. (14) Troisi, A. Chem. Soc. Rev. 2011, 40, 2347−2358. (15) Drukker, K. J. Comput. Phys. 1999, 153, 225−272. (16) Wang, L.; Beljonne, D.; Chen, L.; Shi, Q. J. Chem. Phys. 2011, 134, 244116. (17) Troisi, A.; Orlandi, G. Phys. Rev. Lett. 2006, 96, 086601. (18) Wang, L.; Beljonne, D. J. Chem. Phys. 2013, 139, 064316. (19) Wang, L.; Beljonne, D. J. Phys. Chem. Lett. 2013, 4, 1888−1894. (20) Kubař, T.; Elstner, M. Phys. Chem. Chem. Phys. 2013, 15, 5794− 5813. (21) Kubař, T.; Elstner, M. J. R. Soc., Interface 2013, 10, 20130415. (22) Kubař, T.; Elstner, M. J. Phys. Chem. B 2010, 114, 11221− 11240. (23) Lüdemann, G.; Solov’yov, I. A.; Kubař, T.; Elstner, M. J. Am. Chem. Soc. 2015, 137, 1147−1156. (24) Lüdemann, G.; Woiczikowski, P. B.; Kubař, T.; Elstner, M.; Steinbrecher, T. B. J. Phys. Chem. B 2013, 117, 10769−10778. (25) Woiczikowski, P. B.; Steinbrecher, T.; Kubař, T.; Elstner, M. J. Phys. Chem. B 2011, 115, 9846−9863. (26) Bashir, A.; Heck, A.; Narita, A.; Feng, X.; Nefedov, A.; Rohwerder, M.; Müllen, K.; Elstner, M.; Wöll, C. Phys. Chem. Chem. Phys. 2015, 17, 21988−21996. (27) Heck, A.; Kranz, J. J.; Kubař, T.; Elstner, M. J. Chem. Theory Comput. 2015, 11, 5068−5082. (28) Karl, N. High Purity Organic Molecular Crystals. In Organic Crystals, Germanates, Semiconductors; Freyhardt, H. C., Ed.; Springer: Berlin, 1980; Vol. 4, pp 1−100. (29) Karl, N.; Marktanner, J. Mol. Cryst. Liq. Cryst. Sci. Technol., Sect. A 2001, 355, 149−173. (30) Gao, W.; Kahn, A. J. Phys.: Condens. Matter 2003, 15, S2757− S2770. (31) Kang, J.-W.; Lee, D.-S.; Park, H.-D.; Kim, J. W.; Jeong, W.-I.; Park, Y.-S.; Lee, S.-H.; Go, K.; Lee, J.-S.; Kim, J.-J. Org. Electron. 2008, 9, 452−460. 3095


Article

Journal of Chemical Theory and Computation (66) Bao, Z.; Dodabalapur, A.; Lovinger, A. J. Appl. Phys. Lett. 1996, 69, 4108−4110. (67) Erb, T.; Zhokhavets, U.; Gobsch, G.; Raleva, S.; Stühn, B.; Schilinsky, P.; Waldauf, C.; Brabec, C. J. Adv. Funct. Mater. 2005, 15, 1193−1196. (68) Kline, R. J.; McGehee, M. D.; Kadnikova, E. N.; Liu, J.; Fréchet, J. M. J. Adv. Mater. 2003, 15, 1519−1522. (69) Kline, R. J.; McGehee, M. D.; Kadnikova, E. N.; Liu, J.; Fréchet, J. M. J.; Toney, M. F. Macromolecules 2005, 38, 3312−3319. (70) Merlo, J. A.; Frisbie, C. D. J. Phys. Chem. B 2004, 108, 19169− 19179. (71) Merlo, J. A.; Frisbie, C. D. J. Polym. Sci., Part B: Polym. Phys. 2003, 41, 2674−2680. (72) van de Craats, A. M.; Warman, J. M.; Fechtenkötter, A.; Brand, J. D.; Harbison, M. A.; Müllen, K. Adv. Mater. 1999, 11, 1469−1472. (73) Warman, J. M.; Van De Craats, A. M. Mol. Cryst. Liq. Cryst. 2003, 396, 41−72. (74) Fischbach, I.; Pakula, T.; Minkin, P.; Fechtenkötter, A.; Müllen, K.; Spiess, H. W.; Saalwächter, K. J. Phys. Chem. B 2002, 106, 6408− 6418. (75) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. J. Comput. Chem. 1992, 13, 1011−1021. (76) Ren, J.; Vukmirović, N.; Wang, L.-W. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 205117.

3096


Simulation of Temperature-Dependent Charge Transport in Organic

Recommend Documents