Describing charge carrier anisotropy in crystalline organic semiconductors with *ab initio* methods is challenging because of the weak intermolecular interactions that lead to both localized and delocalized charge transport mechanisms. Small polaron hopping models (localized) are generally used to describe materials with small charge carrier mobilities, while periodic band models (delocalized) are used to describe materials with high charge carrier mobilities. Here, we prove the advantage of applying the constant relaxation time approximation of the Boltzmann transport equation (BTE) to efficiently predict the anisotropic hole mobilities of several unsubstituted (anthracene, tetracene, pentacene, and hexacene) and substituted (2,6-diphenylanthracene, rubrene, and TIPS-pentacene) high-mobility *n*-acene single crystals. Several density functionals are used to optimize the crystals, and the composite density functional PBEsol0-3c/sol-def2-mSVP predicts the most experimentally similar geometries, adequate indirect bandgaps, and the theoretically consistent n-acene charge transport mobility trend. Similarities between BTE and Marcus mobilities are presented for each crystal. BTE and Marcus charge carrier mobilities computed at the same geometry result in similar mobility trends, differing mostly in materials with more substitutions or structurally complex substituents. By using a reduced number of calculations, BTE is able to predict anisotropic carrier mobilities efficiently and effectively for a range of high-mobility organic semiconductors.

## I. INTRODUCTION

Organic semiconductors (OSCs) are a promising class of materials for cutting edge electronic applications, including light emitting diodes (OLEDs), field effect transistors (OFETs), photovoltaics (OPV), and chemical sensors.^{1–5} Their reasonable cost, processability, flexibility, and synthetic and electronic tunability make them an attractive alternative to traditional inorganic semiconductors. One drawback of organic semiconductors is that the molecules comprising the crystalline lattice are weakly bound to one another through soft van der Waals interactions.^{6,7} These weak interactions lead to charge carrier anisotropy, packing defects or crystalline irregularities, and phonon-induced transient charge carrier localization.^{6,8–18} Combinations of these effects significantly limit the maximum charge carrier mobilities of OSCs to generally less than 10 cm^{2} V^{−1} s^{−1}.^{19}

Because of the complicated interplay of electronic structure and crystalline symmetry, predicting anisotropic conductivities and charge carrier mobilities in crystalline organic semiconductors is rather difficult and computationally expensive.^{10,14,20–23} It is common to assume that the charge carrier dynamics are dominated by a completely localized monomer-to-monomer hopping mechanism.^{9,24,25} In the Marcus formalism of small polaron theory, charge carriers hop between the frontier orbitals (HOMO to HOMO for holes, LUMO to LUMO for electrons) of nearest neighbor monomers in a lattice. Hopping rates, *W*_{i}, are usually predicted with the Marcus–Hush equation [Eq. (1)],

where *V*_{i} is the effective electronic coupling between the nearest neighbor monomers of the bulk lattice, *λ* is the total reorganization energy, which is commonly approximated as the intramolecular portion only, as the intermolecular reorganization energy is often assumed to be negligible, *T* is the temperature, and *ℏ* is the reduced Planck’s constant.^{9} The expression assumes that the free energy difference (Δ*G*) to exchange a charge carrier between two nearest neighbor monomers is zero. Additional details regarding how the parameters *V*_{i} and *λ* are computed with density functional theory are presented in Appendix A.

In assuming that charge carriers probabilistically hop through various channels in materials with some directional preference, mobility anisotropy [Eq. (2)] can be described by a Marcus hopping rate by assuming Einstein carrier diffusion.^{9} We have

where *e* is the fundamental charge of the carrier, *n* is the number of spatial dimensions, *r*_{i} is the hopping distance between the mass-weighted centers of two monomers that make up the *i*th dimer, *P*_{i} is the probability (*P*_{i} ≡ *W*_{i}/∑_{i}*W*_{i}) to hop in the *i*th direction, and *F*_{i}(*θ*) is an approximate decay function (usually cos^{2}) away from the mass centered vectors between two monomers.

Critically, this model assumes that isolated monomers and their corresponding dimer couplings correctly capture the behavior of the crystalline bulk. These assumptions are problematic for highly ordered organic crystals such as pentacene and rubrene whose charge carriers dynamically delocalize across many molecules of the lattice (∼17 in pentacene and ∼14 in rubrene).^{26,27} In such materials, the Marcus hopping assumption underestimates charge carrier mobilities by enforcing maximal localization. Additionally, hopping models become invalid when the electronic couplings are greater than half of the reorganization energy (*V*_{i} > *λ*/2).^{28} In this limit, the hopping activation energy approaches zero and therefore the small polaron does not exist. Thus, individual monomer localized states are a bad representation of the carrier in high-mobility materials, making the Marcus hopping mobilities strictly undefined. This would be the case for organic crystals that are conductive enough to be relevant for semiconductor applications, in fact tetracene, pentacene, hexacene, and rubrene already surpass this limit.

In an attempt to approximate the charge transport dynamics of more delocalized materials, affordable computational models beyond the Marcus hopping approximation are required. Band theory based models, like the constant relaxation time Boltzmann Transport Equation (BTE), provide a completely delocalized carrier picture. BTE may better describe the carriers in high-mobility organic crystals where carriers delocalize across multiple molecules. Band theory mobilities depend on the band velocity of a carrier until it is scattered. For organic crystals, the mean free path of charge carriers is significantly smaller ($\u223c1.0$ Å) than typical lattice vector magnitudes or the intermolecular distances, and thus the Mott–Ioffe–Regel (MIR) limit is violated.^{28,29} At this limit, band theory breaks down and the predicted BTE conductivities are probably not representative of the true carrier dynamics.

In recent years, more robust and physically relevant models, such as the transient localization theory, have been developed to describe charge transport in high-mobility organic semiconductor materials that are between the limits of *V*_{i} > *λ*/2 and MIR.^{11,23,26,27} These models predict mobilities that agree exceptionally well with what is observed experimentally. For screening purposes, though, these models are rather computationally demanding and are probably less practical because of this. As an alternative to Marcus hopping mobilities for screening organic semiconductors, this contribution benchmarks BTE as an initial screening tool for predicting the mobility trends of new organic semiconductor materials. Even though the MIR limit is surpassed in these materials, we show that BTE (*τ* = 10 fs) predicts the mobilities of a set of unsubstituted and substituted oligoacene single crystals (Fig. 1) in agreement with experiment and to those predicted by some higher level methods, providing an efficient way to benchmark new transport theories.

## II. THEORETICAL METHODOLOGY

Since the mobility of charge carriers in conductive or semiconductive materials is approximately proportional to the electrical conductivity divided by the charge carrier concentration, it is possible to predict mobilities from the semiclassical BTE conductivity.^{30–32} The BTE predicts the electrical conductivity (*σ*) through various paths of a lattice from a periodic density functional theory (DFT) band structure [Eq. (3)],

where Ξ is the k-point (**K**) dependent transport distribution function describing the *α* direction through the crystal. *f* is the Fermi–Dirac distribution, which depends on the energy (*E*), the Fermi-level $(F)$ or chemical potential, and the temperature (*T*). Here, the transport distribution function [Ξ, Eq. (4)] given by

is defined in terms of the mean charge carrier scattering time (*τ*) and the band velocities (*v*_{i,α}), which are the band energy derivatives along reciprocal space vectors that map to the Cartesian directions of a lattice (i.e., *α* = x, y, or z). Only calculating conductivities along Cartesian directions rather than any reciprocal space vector is a limitation of the current CRYSTAL17 implementation. Finally, the Fermi-level BTE conductivities [Eq. (3)] divided by carrier concentration [number of carriers ($NF,T$ at the Fermi-level) per volume (**V**)] provides a direct way to calculate the anisotropic hole mobilities [Eq. (5)] with periodic DFT. We have

Within the CRYSTAL17 code, the electrical conductivities are calculated with the constant relaxation time (*τ*) BTE.^{33,34} This approximation assumes that charge carrier scattering events, which are expected to be induced predominantly by lattice phonons in softly bound organic crystals, can be represented by an averaged, directionally invariant, constant scattering time.^{30} Directly computing approximate electron–phonon scattering rates (1/*τ*) is significantly more computationally demanding.^{31} To limit the computational cost of the present model, we further assume that charge carrier relaxation times are approximately equivalent among the selected set of crystals. One study approximated *τ* at 300 K to be $\u223c2fs$ in napthalene.^{29,35} It is common to affordably approximate *τ* with a deformation potential model.^{36,37} This approximation includes the effects from acoustic phonon scattering, but it ignores optical phonon effects. Using the deformation potential model at the B3LYP/6–31++G(d,p) level of theory, Kobayashi *et al.* predict *τ* to be 43 and 19 fs for pentacene and rubrene, respectively. As a middle ground, we set *τ* to be 10 fs in all the crystals, allowing for easy scaling.

In this study, the BTE was used to predict directional mobilities in a set of high-mobility single crystals, including four unsubstituted acenes (anthracene, tetracene, pentacene, and hexacene) and three substituted acenes with various axial substitutions (2,6-diphenylanthracene, rubrene, and TIPS-pentacene) (Fig. 1). First, the various crystal structures were independently optimized with several density functionals. The DFT optimized geometries are important since the electronic structure and the resulting transport properties are highly sensitive to the geometry of the system.^{6,7,38,39} In particular, to reproduce or understand experimental anisotropic observations, it is crucial that the optimized structures be close to that of the experimental geometries, symmetries, and lattice parameters, so that the predicted electronic structure and charge transport properties are physically relevant. Here, atom only optimizations with B3LYP/pob-TZVP-rev2 and full (i.e.*,* atoms plus lattice vectors) geometry optimizations with B3LYP/pob-TZVP-rev2, PBE-D3/pob-TZVP-rev2, PBE0-D3/pob-TZVP-rev2, PBEsol0-3c/sol-def2-mSVP, and HSEsol-3c/sol-def2-mSVP were performed.^{40–48} Grimme’s D3 corrected density functionals were included to correct for the weak but important intermolecular dispersion in all but the B3LYP structures.^{49} Additionally, solid-state composite functionals (3c) that combine small basis sets, geometrical counterpoise corrections, and dispersion corrections (two and three body dispersion) were tested because of their computational affordably and good solid-state electronic descriptions.^{48,49} For all of the optimizations, the DIIS convergence accelerator was used. An anisotropic *k*-mesh was defined separately for each crystal such that each real space primitive lattice vector times its corresponding SHRINK factor was greater than 30 Å. The anisotropic *k*-mesh was quadrupled for all of the BTE computations to increase the accuracy of the electronic structure for the band velocities. The tolerance on change in total energy (TOLDEE) was set to 10^{−8} for all calculations. Truncation criteria for the bi-electronic integrals were set with the keyword TOLINTEG to 8 8 8 8 16. Periodic vibrational frequency calculations were conducted to ensure that the optimized structures were true local minima. The BTE mobilities were then compared with Cartesian-axis projected Marcus–Hush–Goddard mobilities to reveal the similarities and differences between the models for each crystal by closely inspecting the mobilities and their trends.

## III. RESULTS

### A. Relevant geometric structure

The mean absolute relative error [MARE, Eq. (6)] in the primitive cell volume (*V*) shows how each crystalline optimization procedure compares with the experimental lattice volumes (Fig. 2 and Table I). Mean relative errors (MREs) are included in Table I to portray the sign errors across the set. It is clear that the solid-state functionals (HSEsol and PBEsol0) deviate the least from the experimental lattice volumes, with MARE of 2.51% and 3.65%, respectively. Note that

. | . | . | Bandgap . | Bandgap . |
---|---|---|---|---|

Method . | Vol. [MARE (%)]^{a}
. | Vol. [MRE (%)]^{b}
. | [MAE (eV)]^{c}
. | [ME (eV)]^{d}
. |

B3LYP (atom only) | ⋯ | ⋯ | 0.30 | −0.30 |

B3LYP | 13.77 | +13.77 | 0.18 | −0.11 |

PBE-D3 | 9.55 | −9.55 | 1.46 | −1.46 |

PBE0-D3 | 10.91 | −10.91 | 0.36 | −0.36 |

PBEsol0-3c | 3.65 | +3.65 | 0.21 | +0.003 |

HSEsol-3c | 2.51 | +2.38 | 0.61 | −0.61 |

. | . | . | Bandgap . | Bandgap . |
---|---|---|---|---|

Method . | Vol. [MARE (%)]^{a}
. | Vol. [MRE (%)]^{b}
. | [MAE (eV)]^{c}
. | [ME (eV)]^{d}
. |

B3LYP (atom only) | ⋯ | ⋯ | 0.30 | −0.30 |

B3LYP | 13.77 | +13.77 | 0.18 | −0.11 |

PBE-D3 | 9.55 | −9.55 | 1.46 | −1.46 |

PBE0-D3 | 10.91 | −10.91 | 0.36 | −0.36 |

PBEsol0-3c | 3.65 | +3.65 | 0.21 | +0.003 |

HSEsol-3c | 2.51 | +2.38 | 0.61 | −0.61 |

^{a}

Mean absolute relative error.

^{b}

Mean relative error.

^{c}

Mean absolute error.

^{d}

Mean error.

B3LYP full optimizations (nuclear coordinates and crystalline lattice vectors) resulted in consistent volume expansions (MRE = +13.77), while the dispersion corrected (D3) PBE (MRE = −9.55) and PBE0 (MRE = −10.91) functionals resulted in geometries whose unit cells were consistently contracted. The B3LYP expansions were expected since functionals without dispersion corrections are known to underestimate the intermolecular London-dispersion attractive forces.^{49} However, dispersion corrections lead to contraction of the unit cell due to an overestimation of intermolecular attractions, as in PBE-D3 and PBE0-D3. Over contraction may be due to DFT neglecting temperature dependent lattice expansion effects.^{50} These compacted unit cells likely overestimate the conductivity, making them unsuitable for experimental comparison. HSEsol-3c (MRE = +2.38) and PBEsol0-3c (MRE = +3.65) only slightly overestimate the volumes. Thus overall, the solid-state composite functionals (3c) provide the most experimentally relevant optimized unit cells.

### B. Capturing experimental electronic structure

Since the transport properties are directly derived from the DFT band structure, a correct electronic structure is required. Predicted indirect optical bandgaps provide a good assessment of the overall electronic structure for each computational method. Here, experimental solution phase optical bandgaps of anthracene, tetracene, pentacene, rubrene, and TIPS-pentacene were used for comparison to avoid inconsistencies in reported solid-state bandgaps due to variations in preparation.^{51,52} Larger mean absolute errors [MAEs, Eq. (7)] of the bandgaps (*BGs*) indicate that the method used provides unrepresentative electronic structure. Note that

Of all the functionals tested, the fully optimized B3LYP (*MAE* = 0.18 eV) and PBEsol0-3c (*MAE* = 0.21 eV) resulted in the smallest mean absolute errors in the bandgap. The PBE-D3 (*MAE* = 1.46 eV) and HSEsol0-3c (*MAE* = 0.61 eV) functionals consistently underestimated the bandgap for all of the crystals in the test set. While the fully optimized B3LYP has the smallest bandgap MAE, its large geometric MARE makes it a poor method overall. Thus, fully optimized PBEsol0-3c provides the most experimentally relevant geometric and electronic structure of all the DFT methods tested.

### C. Predicting charge carrier mobility anisotropy with BTE

Previously reported theoretical hole mobilities for the unsubstituted acene series suggest that the maximum charge carrier mobilities should increase as a function of acene backbone lengthening both for Marcus models^{53} and fragment orbital-based surface hopping (FOB-SH) models.^{26} After determining a computational method that provides both good geometric and electronic structure, the mobility-lengthening trend indicates a functional that provides representative BTE mobilities. The largest Fermi-level hole mobilities of any Cartesian direction (i.e.*,* of the diagonal matrix elements of the conductivity tensor) for each method mentioned above (Fig. 3) compare reasonably well with those reported by Deng and Goddard.^{53}

The B3LYP atom only, PBEsol0-3c, and HSEsol-3c functionals all predict rubrene to have the highest hole mobility followed by hexacene, TIPS-pentacene, or DPA. All of the density functionals tested, except for the B3LYP full optimizations (purple), correctly describe the increasing mobility as a function of the lengthening of the unsubstituted acene backbone. Unsurprisingly, the expansion and contraction of unit cells are correlated with either a decrease or increase, respectively, in the magnitude of calculated mobilities.

Unit cell contraction effects became particularly problematic and obvious for the non-hybrid PBE-D3, whose structural compaction suggested that hexacene would behave like a semimetal, with nonzero conductivity through the bandgap (Fig. 4). This unphysical prediction suggests that non-hybrid density functionals should be avoided for small bandgap materials. The fully optimized B3LYP unsubstituted acene BTE mobilities are the only method we tested that did not have the acene lengthening trend. However, the mobilities do match reported maximum hole mobilities from a deformation potential model at the LDA level of theory, which predicts mobilities of 42.2, 92.5, and 55.6 cm^{2} V^{−1} s^{−1} for anthracene, tetracene, and pentacene, respectively.^{36} Since neither LDA nor B3LYP includes dispersion corrections to account for known intermolecular interactions, it is likely that this trend is unrepresentative.

The PBEsol0-3c largest Fermi-level hole mobilities of the three Cartesian directions for pentacene (*μ*_{x} = 6.525 cm^{2} V^{−1} s^{−1}) and rubrene (*μ*_{y} = 14.524 cm^{2} V^{−1} s^{−1}) compare reasonably well with those reported by Kobayashi *et al.*^{37} They report the largest unit cell vector mobilities to be 58 and 51 cm^{2} V^{−1} s^{−1} for pentacene and rubrene, respectively. If we scale our *τ* = 10 fs to match their deformation potential model derived *τ* (*τ*_{PEN} = 43 fs, *τ*_{RUB} = 19 fs), then this implementation of BTE predicts mobilities that are on the same order of magnitude and ordering, 28.058 and 27.596 cm^{2} V^{−1} s^{−1}, respectively. They are also consistent with maximum 2D hole mobilities along the herringbone plane predicted with FOB-SH by Giannini *et al.*^{26} The FOB-SH mobilities for naphthalene, anthracene, and pentacene (2.1, 3.5, and 9.6 cm^{2} V^{−1} s^{−1}, respectively)^{26} are similar in magnitude to the PBEsol0-3c BTE mobilities (Table II). Additionally, FOB-SH predicts rubrene to have a higher maximum 2D mobility than pentacene, in line with BTE predictions (Table II).

. | $\mu xBTE$ . | $\mu yBTE$ . | $\mu zBTE$ . |
---|---|---|---|

Crystal . | (cm^{2}/V s)
. | (cm^{2}/V s)
. | (cm^{2}/V s)
. |

Anthracene | 2.802 | 2.974 | 2.176 |

Tetracene | 4.728 | 2.430 | 0.669 |

Pentacene | 6.525 | 2.687 | 0.454 |

Hexacene | 3.508 | 10.602 | 0.867 |

2,6-Diphenylanthracene | 0.324 | 8.883 | 4.413 |

Rubrene | 0.354 | 14.524 | 5.161 |

TIPS-pentacene | 10.306 | 4.872 | 0.004 |

. | $\mu xBTE$ . | $\mu yBTE$ . | $\mu zBTE$ . |
---|---|---|---|

Crystal . | (cm^{2}/V s)
. | (cm^{2}/V s)
. | (cm^{2}/V s)
. |

Anthracene | 2.802 | 2.974 | 2.176 |

Tetracene | 4.728 | 2.430 | 0.669 |

Pentacene | 6.525 | 2.687 | 0.454 |

Hexacene | 3.508 | 10.602 | 0.867 |

2,6-Diphenylanthracene | 0.324 | 8.883 | 4.413 |

Rubrene | 0.354 | 14.524 | 5.161 |

TIPS-pentacene | 10.306 | 4.872 | 0.004 |

Crystal . | |V_{12}| (eV)
. | |V_{13}| (eV)
. | |V_{23}| (eV)
. | |V_{34}| (eV)
. | |V_{15}| (eV)
. | |V_{26}| (eV)
. | |V_{37}| (eV)
. | |V_{48}| (eV)
. |
---|---|---|---|---|---|---|---|---|

Anthracene | 0.0389 | 0.0286 | 0.0286 | 0.0389 | 0.0004 | 0.0004 | 0.0004 | 0.0004 |

Tetracene | 0.0074 | 0.0119 | 0.0702 | 0.0119 | 0.0015 | 0.0015 | 0.0005 | 0.0005 |

Pentacene | 0.0148 | 0.0288 | 0.0725 | 0.0017 | 0.0011 | 0.0011 | 0.0005 | 0.0005 |

Hexacene | 0.0493 | 0.0755 | 0.1001 | 0.0482 | 0.0005 | 0.0005 | 0.0005 | 0.0005 |

2,6-DPA^{a} | 0.0082 | 0.0562 | 0.0562 | 0.0082 | 0.0061 | 0.0061 | 0.0061 | 0.0061 |

Rubrene | 0.0699 | 0.0170 | 0.0170 | 0.0699 | 0.0025 | 0.0025 | 0.0025 | 0.0025 |

TIPS-pentacene | 0.0169 | 0.0137 | |V_{14}| = 0.0007 | 0.0169 | 0.0001 | 0.0001 | 0.0001 | 0.0001 |

Crystal . | |V_{12}| (eV)
. | |V_{13}| (eV)
. | |V_{23}| (eV)
. | |V_{34}| (eV)
. | |V_{15}| (eV)
. | |V_{26}| (eV)
. | |V_{37}| (eV)
. | |V_{48}| (eV)
. |
---|---|---|---|---|---|---|---|---|

Anthracene | 0.0389 | 0.0286 | 0.0286 | 0.0389 | 0.0004 | 0.0004 | 0.0004 | 0.0004 |

Tetracene | 0.0074 | 0.0119 | 0.0702 | 0.0119 | 0.0015 | 0.0015 | 0.0005 | 0.0005 |

Pentacene | 0.0148 | 0.0288 | 0.0725 | 0.0017 | 0.0011 | 0.0011 | 0.0005 | 0.0005 |

Hexacene | 0.0493 | 0.0755 | 0.1001 | 0.0482 | 0.0005 | 0.0005 | 0.0005 | 0.0005 |

2,6-DPA^{a} | 0.0082 | 0.0562 | 0.0562 | 0.0082 | 0.0061 | 0.0061 | 0.0061 | 0.0061 |

Rubrene | 0.0699 | 0.0170 | 0.0170 | 0.0699 | 0.0025 | 0.0025 | 0.0025 | 0.0025 |

TIPS-pentacene | 0.0169 | 0.0137 | |V_{14}| = 0.0007 | 0.0169 | 0.0001 | 0.0001 | 0.0001 | 0.0001 |

^{a}

Diphenylanthracene.

Since the PBEsol0-3c functional resulted in the combination of small errors in the geometries and bandgaps and the correct directional mobility trends of the unsubstituted acenes, we chose to study its optimized crystal structures more closely. The PBEsol0-3c Cartesian directional conductivities (Fig. 5) show a range of variation across the test set. In particular, the y direction in rubrene or anthracene and the x direction in tetracene or pentacene are predicted to be significantly more conductive than the other directions in those crystals. The largest Fermi-level hole mobilities [gray-dashed line, computed with Eq. (5) and Table II] highlight the varying anisotropy among the crystals, even those with similar herringbone structures.

The directionally dependent mobility ratios provide a quantitative metric of conduction/mobility anisotropy, where ratios close to one indicate isotropic carrier diffusion and ratios much larger than one indicate anisotropic carrier diffusion. Of all the groups, anthracene (*μ*_{y}/*μ*_{x} = 1.061, *μ*_{y}/*μ*_{z} = 1.367) is predicted to be the most isotropic of the unsubstituted acenes according to the Cartesian directions measured. Both tetracene (*μ*_{x}/*μ*_{y} = 1.946, *μ*_{x}/*μ*_{z} = 7.067) and pentacene (*μ*_{x}/*μ*_{y} = 2.428, *μ*_{x}/*μ*_{z} = 14.37) share similar preferential anisotropic transport in the x direction of the lattice, while hexacene (*μ*_{y}/*μ*_{x} = 3.022, *μ*_{y}/*μ*_{z} = 12.23) prefers conducting in the y direction. However, it is important to note that the y direction of hexacene is symmetrically equivalent to the x direction in tetracene and pentacene, making these results directly comparable. For the substituted acenes, 2,6-diphenylanthracene (*μ*_{y}/*μ*_{z} = 2.013, *μ*_{y}/*μ*_{x} = 27.417) and rubrene (*μ*_{y}/*μ*_{z} = 2.814, *μ*_{y}/*μ*_{x} = 41.03) prefer conducting in the y direction, while TIPS-pentacene prefers conducting in the x direction (*μ*_{x}/*μ*_{y} = 2.115, *μ*_{x}/*μ*_{z} = 2576.5). Of the set of crystals, TIPS-pentacene and rubrene show the most anisotropic charge transport.

### D. Comparison of BTE and Marcus–Hush–Goddard mobilities

The BTE band derivatives can be interpreted as an average velocity of a charge carrier (i.e., hole) in some *i*th band (i.e., highest valence band) along a direction (*α*) in the absence of external forces. In any material at finite temperature, phonons are expected to disrupt (or scatter) charge carriers and therefore govern their average lifetime(s) (which BTE defines as constant *τ*). In contrast, Marcus theory predicts hopping rates that depend on treating each carrier hop as a separate event with some finite probability of occurring. The Marcus–Hush–Goddard mobility then comes from a random walk of uncorrelated hopping events between neighboring molecules in a material or diffusion of the carrier.

Conceptually, BTE assumes complete delocalization, while Marcus hopping assumes complete localization of the carrier. Both theories calculate a maximum carrier mobility that is reduced through scattering or geometric reorganization, respectively. One might expect that as charge carrier lifetimes decrease (i.e., the constant scattering limit) or as the hopping rates increase (i.e., at a constant measurable velocity), the mobilities predicted by the two models would converge. However, because both models break down at their respective limits (*V*_{i} > *λ*/2 and MIR, *vide supra*), the error in the predicted mobilities may be very large. As the mobility trends of the selected materials have been well described by Marcus rates previously, we use these mobilities to compare our BTE mobilities against. We keep the geometry consistent in both models to make the comparison as appropriate as possible.

First, the hopping rates for each unique dimer pair of a crystal are calculated according to the effective electronic couplings (Table III, computed via Gaussian16 and CATNIP) and the isolated monomer reorganization energy ( Appendix A).^{54,55}

Crystal . | W_{avg(12,34)} (THz)
. | W_{13} (THz)
. | W_{23} (THz)
. | W_{ortho} (THz)
. |
---|---|---|---|---|

Anthracene | 15.350 | 8.342 | 8.342 | 0.002 |

Tetracene | 1.685 | 2.033 | 71.004 | 0.031 |

Pentacene | 4.035 | 15.419 | 97.988 | 0.019 |

Hexacene | 54.852 | 131.577 | 231.610 | 0.006 |

2,6-Diphenylanthracene | 0.510 | 23.726 | 23.726 | 0.280 |

Rubrene | 44.255 | 2.625 | 2.624 | 0.054 |

TIPS-pentacene | 3.054 | 2.071 | W_{14} = 0.005 | 0.000 |

Crystal . | W_{avg(12,34)} (THz)
. | W_{13} (THz)
. | W_{23} (THz)
. | W_{ortho} (THz)
. |
---|---|---|---|---|

Anthracene | 15.350 | 8.342 | 8.342 | 0.002 |

Tetracene | 1.685 | 2.033 | 71.004 | 0.031 |

Pentacene | 4.035 | 15.419 | 97.988 | 0.019 |

Hexacene | 54.852 | 131.577 | 231.610 | 0.006 |

2,6-Diphenylanthracene | 0.510 | 23.726 | 23.726 | 0.280 |

Rubrene | 44.255 | 2.625 | 2.624 | 0.054 |

TIPS-pentacene | 3.054 | 2.071 | W_{14} = 0.005 | 0.000 |

Within the Goddard extension of the Marcus–Hush transport model, mobilities are predicted along particular conduction channels via weighted sums of dimer couplings.^{9} The unique dimer hopping rates, *W*_{i}, decay according to a cos^{2} function away from the vector between the centers of mass of each monomer, i.e.*,* the dimer–center-of-mass (d-COM) vector. Because the d-COM directions are defined by the original dimer orientations, the unique d-COM vectors are normalized and projected along the Cartesian directions (details in Appendix B). Importantly, the Cartesian projected mobilities are normalized so that only unique conduction channels are combined and all nearest neighbor dimer pairs are included.

The Marcus–Hush–Goddard mobilities (Table V) show similar trends to the BTE predicted mobilities (Table II) but have magnitudes that are significantly smaller. In addition, the relative directional magnitudes of a particular crystal seem to be highly affected by the bulk effects included in BTE. This is especially true in the directions that do not align with the high-mobility channels of the organic crystals, where BTE predicts significantly higher mobilities than Marcus–Hush–Goddard along these directions, more in line with the measured anisotropic mobilities.^{16}

Crystal . | $\mu xM$ (cm^{2}/V s)
. | $\mu yM$ (cm^{2}/V s)
. | $\mu zM$ (cm^{2}/V s)
. |
---|---|---|---|

Anthracene | 0.0056 | 0.6122 | 0.1451 |

Tetracene | 1.8850 | 1.1842 | 0.0000 |

Pentacene | 2.3107 | 1.4996 | 0.0000 |

Hexacene | 3.6215 | 4.9251 | 0.0000 |

2,6-Diphenylanthracene | 0.0012 | 0.6552 | 0.4629 |

Rubrene | 0.0000 | 4.0435 | 0.0297 |

TIPS-pentacene | 0.2871 | 0.1004 | 0.0000 |

Crystal . | $\mu xM$ (cm^{2}/V s)
. | $\mu yM$ (cm^{2}/V s)
. | $\mu zM$ (cm^{2}/V s)
. |
---|---|---|---|

Anthracene | 0.0056 | 0.6122 | 0.1451 |

Tetracene | 1.8850 | 1.1842 | 0.0000 |

Pentacene | 2.3107 | 1.4996 | 0.0000 |

Hexacene | 3.6215 | 4.9251 | 0.0000 |

2,6-Diphenylanthracene | 0.0012 | 0.6552 | 0.4629 |

Rubrene | 0.0000 | 4.0435 | 0.0297 |

TIPS-pentacene | 0.2871 | 0.1004 | 0.0000 |

For all of the materials considered, with the exception of anthracene, the predicted ordering of the Cartesian carrier mobilities match between both models (Tables II and V). The anthracene ordering mismatch can be attributed, at least partially, to both how similar its directional conductivity magnitudes are at the Fermi-level and to the fact that the directional ordering changes just below the Fermi level with *μ*_{z} < *μ*_{x} (Fig. 5, gray dotted line). Of course, the Marcus dimer projection model does not include the band curvature at the Fermi level as it assumes only localized charges and thus does not capture this band-like unique phenomenon. For all of the other crystals, reordering of the largest conductivity does not occur near the Fermi level, if at all.

The similarities between BTE and Marcus–Hush–Goddard models are captured by plotting the largest Fermi-level Cartesian mobility of each crystal on the same axis (Fig. 6). The primary difference between the Marcus and BTE substituted vs unsubstituted results arises from the hole induced reorganization energy (Fig. 7). In the Marcus model, a hole is localized on a single molecule as it probabilistically hops through the lattice via channels of minimal resistance. In many cases, addition of a localized hole can significantly change the energy and geometry of the isolated monomers. This is captured in the Marcus reorganization energies (*λ*) and leads to large reductions in the predicted hopping rates. Two trends are observed in the monomer reorganization energy results: (1) The reorganization energy tends to decrease as the length of the n-acene backbone increases and (2) all of the substituted acenes, when compared with their unsubstituted counterparts of the same acene backbone length, are predicted to have substantially larger reorganization energies. So, as molecules of the crystal become less rigid, either with more substitution (i.e.*,* from di- to tetra-phenyl substitutions in 2,6-diphenylanthracene and rubrene) or with more substituent flexibility (i.e.*,* from phenyls to TIPS), the difference in their reorganization energy is expected to increase. Thus, both structural rigidity and hole delocalization length should increase the mobility of organic single crystals.

The periodicity of the BTE model allows for the delocalization of holes, across the whole unit cell. Thus, constant relaxation time BTE assumes that the reorganization energy is approximately zero because of this charge spreading. Since the charge carrier is completely delocalized, movement of the charge leads to less geometric reorganization as it propagates. This, of course, may be an oversimplification of true charge transport mechanism in organic crystals where phonon-induced transient localization events may introduce some significant reorganization.^{11}

The predicted mobilities from the two models might then be expected to converge for longer acenes, as the Marcus reorganization energy steeply decreases as a function of acene backbone length. Surprisingly, this was not observed, and the gap between the BTE and Marcus mobilities of the unsubstituted acenes increased as a function of backbone lengthening (Fig. 6). This directly shows that a different constant relaxation time cannot align the two theories quantitatively. Again, this means that the implementations of these two theories do not converge at fast scattering rates (small *τ*).

### E. Computational costs

Predicting mobilities with both BTE and Marcus requires multiple DFT energy calculations. Thus, the computational efficiency should be measured in total CPU-hr cost (post optimization) to calculate the Cartesian mobilities with either model (Fig. 8). BTE conductivities require a dense k-point grid SCF calculation (denser than the optimization grid) and a finite temperature (300 K) Cartesian conductivity calculation. Marcus effective electronic couplings are derived from single point energy calculations of all of the unique monomers and dimers (the number depends on the crystalline symmetry). Marcus reorganization energies are estimated by four additional charge localization calculations (2 optimizations + 2 single point energies) (see Appendix A). Figure 8 shows the cost of the best computational methodology for each model. The BTE/PBEsol0-3c/def2-mSVP (orange) cost is significantly lower than Marcus/PBE0-D3/def2-TZVP (gray). However, these are of different basis set sizes, which might overestimate the advantage of the BTE method. So, in addition, the cost of BTE/PBE0-D3/pob-TZVP-rev2 (green) is included. These results clearly show the benefit of using the modified double-*ζ* basis over the triple-*ζ* and the cost difference between the BTE and Marcus models for each crystal. On average, the BTE triple-*ζ* basis costs 2.65 ± 0.38 times that of the modified double-*ζ* basis at a comparable level of theory. Marcus mobilities, on the other hand, cost 4.59 ± 8.42 times the average BTE cost. The large standard deviation in the cost of each method arises primarily from the highly substituted rubrene and TIPS-pentacene crystals. In these crystals, the cost to calculate the reorganization energy is substantially larger than in the other crystals. In fact, it costs 13.3 and 48.5 times the average cost in the case of the other crystals, respectively. This result is expected because configurational flexibility is correlated with computational cost. Therefore, for large and highly flexible systems, the BTE model is even more computationally efficient than the Marcus model.

## IV. CONCLUSION

The BTE model presented here provides an efficient and effective way to compute relative anisotropic conductivities and charge carrier mobilities of high-mobility organic semiconductor materials. The BTE model when coupled with an affordable composite functional (PBEsol0-3c) correctly predicted experimental geometries, bandgaps, and the theoretically consistent increasing mobilities across the n-acenes. In addition, the overall anisotropic mobility trends of the dimer Marcus mobilities were reproduced with the BTE model. The magnitude of the predicted mobilities of the substituted acenes with both models is significantly different due to the larger reorganization energies in Marcus theory. Overall, the BTE model consistently provides quantitative and qualitative ordering of mobilities within and between various acene crystals with less computational cost and numerical manipulation, providing an efficient screening method for large sets of crystals.

## ACKNOWLEDGMENTS

We acknowledge the financial support received from Lehigh University and research computing resources provided by Lehigh University partially supported by NSF CC* Compute program through Grant Nos. OAC-2019035 and TG-CHE190011 allocation from Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1548562.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Zachary J. Knepp**: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Methodology (lead); Validation (equal); Visualization (lead). **Gabriel B. Masso**: Data curation (supporting); Validation (supporting). **Lisa A. Fredin**: Conceptualization (equal); Software (lead); Supervision (lead); Validation (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

The optimized structures are available at https://github.com/fredingroup/organics/tree/main/boltztra as CIF files.

### APPENDIX A: MARCUS REORGANIZATION ENERGIES AND ELECTRONIC COUPLINGS

It is common that intramolecular reorganization energies *λ*_{intra} are computed using a four point method that utilizes the adiabatic potential energy surfaces of a single isolated monomer from the crystal structure.^{9} We have

The states $\varphi opt0$ and $\varphi opt+$ correspond to the optimized Cartesian coordinates of an isolated monomer in the neutral and charged states, respectively; *H*^{0} and *H*^{+} represent the Kohn–Sham Hamiltonians for the neutral and charge electronic states, respectively. In general, according to Eq. (1), charge hopping rates are maximized when the reorganization energy is minimized.

The effective electronic couplings (*V*_{ab} or *V*_{i}) can be computed with density functional theory (DFT) through projecting the isolated monomer orbitals onto the coupled dimer orbitals (DIPROs).^{24,38,55} Note that

Effective, here, means that in addition to the charge transfer integral (*J*_{ab}), one also takes into consideration the orbital overlap (*S*_{ab}) between the monomers and site energy fluctuations (*ɛ*_{a} and *ɛ*_{b}) among the unique monomers within a given dimer, where

In our calculations, we make use of the CATNIP code to compute the quantities mentioned above.^{55}

### APPENDIX B: RENORMALIZATION OF MARCUS COUPLINGS ALONG CARTESIAN DIRECTIONS

For direct comparison between BTE and Marcus mobilities, the unique d-COM vectors are decomposed into their three orthogonal Cartesian components such that the sums of the square projections are normalized. That is,

Here, *θ*_{iα} is the angle between d-COM vector *i* and a Cartesian vector *α*. The product of the cos^{2} projection and the initial d-COM vector magnitude may be interpreted as an adjusted hopping length in a particular Cartesian direction,

In this assumption, the charge carriers are still expected to hop between the monomer centers of mass in each dimer, on average, so the distance traveled in particular Cartesian directions depends on the projection of this hop or the *r*^{2} · cos^{2} components. Then, the hopping probabilities are weighted according to the unique dimer hopping rates,

Before computing the weights, the nonunique conduction channel directions such as *r*_{12} and *r*_{34} were combined in a weighted average of the hopping rates in the d-COM direction,

In doing this, the actual weighting in the projection [Eq. (B5)] is limited to the unique conducting directions. After weighting, the Cartesian projected mobilities are represented by the product of the hopping probability, hopping weight, and corrected hopping distance squared,

The sum over the *i* indices results in an estimate for the Einstein diffusion coefficient *D*_{α} in each Cartesian direction. This model is identical to the 2D Goddard equation with a single cos^{2} function.^{9}

In the diffusion coefficient sums, hopping rates between the nearest neighbor dimer pairs of an isolated herringbone slice ({*W*_{avg(12,34)}, *W*_{13}, *W*_{23}}) and an end-on-end dimer (*W*_{ortho}) are included. The planar and orthogonal dimers were required to accurately predict mobilities in all the projected Cartesian directions as some Cartesian directions are not well aligned with the high conductivity herringbone slice. All drift mobilities [Eq. (B5)] were computed at 300 K to match the temperature set in BTE.

## REFERENCES

_{2}reactivity with bisphenalenyls for simple, tunable, and multicycle colorimetric oxygen-sensing films

^{2}V

^{−1}s

^{−1}: A breakthrough of organic semiconductors for field-effect transistors

^{2}O, CuO, and NiO from hybrid density functional theory