Thomas F.
Harrelson
a,
Varuni
Dantanarayana
b,
Xiaoyu
Xie
c,
Correy
Koshnick
d,
Dingqi
Nai
a,
Ryan
Fair
e,
Sean A.
Nuñez
f,
Alan K.
Thomas
g,
Tucker L.
Murrey
d,
Michael A.
Hickner
f,
John K.
Grey
g,
John E.
Anthony
h,
Enrique D.
Gomez
e,
Alessandro
Troisi
c,
Roland
Faller
a and
Adam J.
Moulé
*ab
aDepartment of Chemical Engineering, University of California, Davis, CA, USA. E-mail: amoule@ucdavis.edu
bDepartment of Chemistry, University of California, Davis, CA, USA
cDepartment of Chemistry, University of Liverpool, Liverpool, UK
dDepartment of Materials Science and Engineering, University of California, Davis, CA, USA
eDepartment of Chemical Engineering, Pennsylvania State University, University Park, PA, USA
fDepartment of Materials Science and Engineering, Pennsylvania State University, University Park, PA, USA
gDepartment of Chemistry & Chemical Biology, University of New Mexico, Albuquerque, NM, USA
hDepartment of Chemistry, University of Kentucky, Lexington, KY, USA
First published on 23rd October 2018
Recent theories suggest that low frequency dynamic intramolecular and intermolecular motions in organic semiconductors (OSCs) are critical to determining the hole mobility. So far, however, it has not been possible to probe these motions directly experimentally and therefore no unequivocal and quantitative link exists between molecular-scale thermal disorder and macroscale hole mobility in OSCs. Here we use inelastic neutron scattering to probe thermal disorder directly by measuring the phonon spectrum in six different small molecule OSCs, which we accurately reproduce with first principles simulations. We use the simulated phonons to generate a set of electron–phonon coupling parameters. Using these parameters, the theoretical mobility is in excellent agreement with macroscopic measurements. Comparison of mobility between different materials reveals routes to improve mobility by engineering phonon and electron–phonon coupling.
Conceptual insightsThe charge mobility in molecular semiconductors is limited by significant nonlocal electron–phonon coupling caused by intermolecular vibrations, otherwise known as phonons. To engineer new molecular materials with superior functional properties and performance, we need to gain fundamental understanding of how these phonons cause dynamic localization of charge carriers. Here we directly measure the phonon spectrum using high-resolution inelastic neutron scattering (INS) experiments. INS provides a more complete phonon spectrum than analogous optical experiments, such as FTIR or Raman. We demonstrate the accurate computation of nonlocal electron–phonon coupling parameters for six different molecular semiconductors using density functional theory simulations of phonons over the full Brillouin zone, experimentally validated by inelastic neutron scattering experiments. Our coupling parameters were then used to compute charge mobilities which quantitatively agree with literature values, demonstrating our ability to accurately predict functional properties of a variety of molecular semiconductors. Thus, our measured results coupled with computational methodology can be used to screen molecular semiconductors to predict a set of high-performance materials, significantly accelerating the navigation through chemical design space for next-generation materials. |
It is natural to expect that μh can be optimized by exploring the materials space with suitable computational models coupling the atomic scale material properties (describing electron, phonon, and electron–phonon interactions) with a macroscopic theory of charge transport. This has not been realized so far because of challenges on both the computational and the theoretical side of the problem. Boltzmann (band) transport theory is the standard theory to describe charge carriers in inorganic semiconductors. Band theory fails to describe molecular semiconductors with experimental mobilities of the order of few cm2 V−1 s−1 because it predicts a mean free path of just few Angstroms, which is inconsistent with the carrier delocalization assumed by band transport. In contrast, incoherent hopping models fail to be physically reasonable because they predict a hopping rate that is too fast (>1013 s−1). Several recent reviews summarize the problems and the common approach emerging to address them.8,9 The new approach describes the interaction between molecular orbitals localized on different molecules with a transfer integral J, which is modulated by thermal motions. For molecular semiconductors, the standard deviation of the thermal modulation (σ) is of the same order of magnitude as J and takes place with a timescale of ≈1 ps. The thermal fluctuations are a manifestation of non-local electron–phonon coupling but, unlike other simpler cases, the fluctuations of J are too large to be treated as a small perturbation (as in band transport), too slow to be averaged out (e.g. renormalized), and too fast to be treated as static disorder. A variety of numerical approaches have been proposed to deal with this scenario often described as “dynamic disorder”, where the electron–phonon coupling is responsible for a transient localization of the carrier.8,10–13
For any model, it is impossible to predict charge mobility without high quality experimental information about the phonons responsible for the dynamic localization, e.g. mostly low frequency motions with a large intermolecular component. While high energy phonons can be routinely probed via FTIR and/or Raman spectroscopy,14–16 low energy phonons can only be measured using a specialized Raman spectroscopy setup and/or terahertz spectroscopy.17 Even then, the selection rules for Raman limit the measurement to the observation of gamma point phonons, which is of limited use for μh calculations because the symmetry of gamma point phonons minimize the relative displacements between molecules in neighboring cells. Phonons at the Brillouin zone boundaries are most important for the computation of μh. Also terahertz spectroscopy of OSCs is dominated by the dielectric response of free or loosely bound charge carriers present in OSCs, which mask the weaker response of phonons. Diffuse electron scattering has been used to measure the cumulative disorder caused by phonons,18 but this method does not resolve the dynamic disorder into its frequency dependent components and cannot directly inform the charge transport model. This lack of information has so far prevented a quantitative test of dynamic transport theories and the validation of any computed phonon spectra at low frequency. Thus measurement of the full phonon spectrum has been a central blocking point for the development of organic electronics.
To measure the phonons that lead to dynamic disorder, we use inelastic neutron scattering (INS), a technique that probes the atomic dynamics of a material. INS has been used to infer a variety of materials properties including the local structure of polymers,19,20 binding interactions in porous materials,21,22 and protein folding in response to a stimulus. INS is particularly well suited for studying dynamic disorder in OSCs because spectroscopy with neutrons does not have any selection rules (every phonon mode in the Brillouin zone is observable) and neutrons interact with atomic nuclei only (not the free charge carriers).19,23 The low energy vibrations are not convoluted with the dielectric response of the material because there is no neutron-charge carrier interaction. Also the inelastic scattering cross-section is largest for protons, so we observe the weighted density of states of proton motions. Since nearly every carbon motion in an OSC is coupled to the motion of 1H atoms, the INS spectrum contains energy-resolved information of all atomic motions in the OSC crystal. Given an experimentally well-defined crystal structure, the INS instrument at Oak Ridge National Laboratory called VISION characterizes vibrational dynamics with high resolution over a broad energy range (1–600 meV or 8–5000 cm−1) better than other available instruments because the specific instrument design and high incident neutron flux maximizes the signal-to-noise ratio and count rate.24,25
An accurate simulation of the spectrum, e.g. via density functional theory (DFT), provides the most detailed picture of phonon modes in the OSC crystal across all measured energy scales. As van der Waals (vdW) forces strongly impact the crystal structure of molecular semiconductors, and therefore the phonons, we simulate each material with a functional that explicitly includes vdW forces. We use the experimentally verified set of phonons to create a parameter-free model of hole transport in six different OSCs. The model connects structure at the atomic length scale to observable phenomena at the device length scale using information obtained entirely from first principles as shown in Fig. 1, which fits into any computational scheme that that searches through chemical design space to optimize μh.
A numerically efficient method to compute μh is provided by transient localization theory (TLT).8 We chose to use this technique because the theory was specifically designed for modeling high-mobility molecular semiconductors, but any other appropriate theory could be used in conjunction with the computed EP parameters. TLT has been shown to qualitatively reproduce μh for a wide variety of OSCs,26 but has not been applied to materials in which the phonon spectrum is known.
The six materials studied exhibit either high μh, or have slight chemical variations that reduce μh (see crystal structures in Fig. 2 and chemical names in methods section) and can be grouped into two series of homologous compounds. The first group is BTBT, m8-BTBT, and c8-BTBT, with an identical conjugated core and a very similar herring-bone crystal structure in the high mobility plane. The only difference is the inclusion of one (m8-BTBT), or two (c8-BTBT) octyl chains along the longitudinal axis of the BTBT core. Nevertheless, μh of c8-BTBT is much higher than m8-BTBT or BTBT.27 It has been suggested that the octyl side chains reduce the dynamic disorder (σ) of the BTBT repeat unit by preventing lateral vibrations.18 We also study a series of three substituted acenes: TIPS-PN, TES-ADT, and TIPS-ADT to test how changing the conjugated core and side chain structure affects the phonon spectrum. TES-ADT has a nearly identical brickwork crystal structure as TIPS-PN,28,29 but has a higher μh.18,28,29 TES-ADT substituted the terminal benzenes of the core of TIPS-PN with thiophenes, and has smaller side chains than TIPS-PN and TIPS-ADT. A recent article by Illig et al.18 suggested that the charge mobility in TES-ADT is higher due to reduced dynamic motion of the backbone. We also study TIPS-ADT, which assumes a slipped-stack crystal packing and typically displays much lower μh.29,30 Both sets of OSCs are designed to show how differences in crystal structure and side chain attachment affect the phonon spectrum, and therefore the EP coupling.
Most importantly for μh calculations, we accurately simulate the phonon spectrum at energies below kBT at room temperature (<200 cm−1). Fig. 2 shows a comparison between the measured and simulated phonons for each of the six OSCs studied from 10–600 cm−1. The agreement between measurement and simulation is excellent down to the linewidth of the elastic peak (∼10 cm−1) for five of the six samples, particularly below 200 cm−1, indicating that we accurately simulated the interactions between the conjugated cores. Many of the peaks above 200 cm−1 primarily contain motions of side chains, which are expected to display increased structural heterogeneity leading to minor disagreements between the experiments and our simulations.
Our INS results, combined with XRD data in ESI,† Section S1, identify a low temperature polymorphism in c8-BTBT causing an apparent disagreement between the INS and simulated spectra. Since the simulated structure of c8-BTBT agrees with room temperature XRD measurements,35 the simulated phonons correspond to the room temperature motions. The c8-BTBT phonon predictions without experimental INS verification is justified by the quality of agreement for the other five materials.
In addition to the INS data, we measured Raman and FTIR spectra for each of the substituted acene materials (see ESI,† Section S3). Both, FTIR and Raman spectra substantially differ from their INS counterparts, which demonstrates the complementarity between the three measurements. Despite the limitations of Raman and FTIR due to the selection rules, the observed peaks are not heavily weighted by the motions of hydrogens as in INS, and, thus, contain more carbon–carbon, carbon–silicon, carbon–sulfur information than the associated INS spectra. For each substituted acene material, our simulated phonons nicely agree with the FTIR spectra (see ESI,† Section S4), indicating that the simulated phonons accurately describe the full range of motions.
To assess the impact of the phonons on transfer integrals between different pairs of molecules, we computed the non-local electron–phonon coupling, which is a measure of how the transfer integral Jij between molecules i and j is modulated by the displacement QI along phonon mode I. Assuming a linear coupling, the non-local EP coupling parameters gij,I are defined by the relation,
(1) |
Because of the very large number of modes, a convenient way to report the computed EP coupling is by defining a spectral density,
(2) |
We see that each material contains many peaks in Bij(ω) over a very large energy range indicating that many different types of modes contribute to the dynamic disorder. Some materials show the largest peaks in the low energy region of the spectrum (lower than 200 cm−1), indicating that those phonons can be treated classically when considering their impact on charge transport, in agreement with prior studies.10,11,13,36 However, every material contains a non-negligible contribution from high energy modes to total EP coupling, even when considering the phonon occupation numbers at 300 K (see Table 1). The high energy peaks (∼1500 cm−1) in the EP spectra of all materials are the result of carbon–carbon stretch modes in the conjugated core, leading to disruption of the bonding order causing large shifts in the frontier orbital density. The intensity of high energy peaks in the substituted BTBT's appears to be greater than the same peaks in the substituted acenes, in which EP coupling is dominated by acoustic modes.
Material | J (cm−1) | σ 300K (cm−1) | Low freq. (%) | Pair | 1/m* (a.u.−1) | μ theory (cm2 V−1 s−1) | μ exp |
---|---|---|---|---|---|---|---|
c8-BTBT | −728 | 367 | 70 | B and C | 2.21 | 4.86 | 6.0 ± 1.1 |
802 | 210 | 62 | A | ||||
m8-BTBT | −625 | 547 | 76 | B | 2.11 | 1.43 | 2.9 ± 1.1 |
−706 | 535 | 76 | C | ||||
769 | 254 | 67 | A | ||||
BTBT | 1.03 | 39.0 | 96 | B and C | 0.45 | 0.963 | 0.024 ± 0.007 |
765 | 256 | 77 | A | ||||
TES-ADT | −1339 | 439 | 87 | C | 3.40 | 4.22 | 1.0 ± 0.5 |
−462 | 223 | 94 | B | ||||
TIPS-PN | 587 | 311 | 88 | B | 0.71 | 1.34 | 0.65 ± 0.35 |
19.1 | 124 | 74 | C | ||||
TIPS-ADT | −486 | 279 | 97 | A | 0.61 | 0.545 | 0.30 ± 0.25 |
To get an aggregate picture of how phonons impact μh, we find the variance of the transfer integral (σij,T2) at a given temperature using eqn (3).37
(3) |
While σij,T in Table 1 were computed from the phonons sampled over the full Brillouin zone, we also compute σij,T from only the gamma point phonons (see ESI,† Section S4). These gamma point σij,T values differ from the corresponding values in Table 1 by −40% to 25%. Thus, gamma point phonons and measurements that only probe the gamma point (e.g. Raman) are insufficient to describe quantitatively the dynamic disorder in OSCs. To better appreciate the different relevance of the information provided by Raman and INS spectroscopy, we report both experimental spectra along with Bij(ω) in Fig. 3. There are a few peaks in Bij(ω) that are present in the Raman spectra and absent in the corresponding INS spectra. However, the dominant contributions to EP coupling come from the large, broad peak <50 cm−1, which is only captured by INS because INS contains information on phonons over the entire Brillouin zone, whereas Raman scattering only probes gamma point phonons. Each of the six materials displays considerable phonon dispersion at low energy (see ESI,† Section S4). Therefore, INS can be used to probe dynamic disorder in crystalline OSCs, but INS and Raman should be used together to obtain complete information of the low energy phonons.
Fig. 3 Comparison between INS, Raman, and EP coupling spectra (Bij(ω)) for the substituted acene materials (TIPS-PN, TES-ADT, TIPS-ADT). |
Fig. 4 Spectral decomposition of the EP coupling parameter for different pairs of molecules (A, B, or C) for all molecules studied. Sketches of pairs of molecules are shown in the insets of (a–f). |
Our experimentally verified EP couplings were used to parameterize the TLT models of μh for each material. The input for the model are the transfer integrals (Jij) and their standard deviation in the high mobility plane (σij) with a single characteristic fluctuation time, τR (which is extracted from the spectral density as detailed in the ESI,† Section S5). In comparing theoretical and experimental μh one must consider that defects and surface effects coming from different substrates make the measurement of the intrinsic single crystal μh extremely challenging. Also, μh can be highly anisotropic in the high-mobility plane (ESI,† Section S6), meaning that the crystalline orientation with respect to the substrate strongly impacts the measured μh. Each of the materials has been heavily studied in the literature except for BTBT, and the reported μh in Table 1 represents reasonable averages from recent publications that performed transistor measurements to measure μh.28,30,34,38 Since there is no available evidence of intrinsic, band-like transport in BTBT in the literature (e.g. experiments showing a decrease in mobility with temperature), it is possible that the value in Table 1 may be indicative of thermally activated transport and thus not comparable with our theoretical results. In addition, these experimental μhs do not hold for every processing condition/device architecture; they are used as a qualitative guide to determine the accuracy of our simulations.
Since TLT treats all phonons classically, but the role of quantum high frequency modes varies across the materials considered, an improved model should be introduced to differentiate between the effect of low and high frequency phonons. The only phonons able to localize the charge are those with characteristic fluctuation time similar to or slower than the quantum delocalization dynamics of the charges. Faster phonons can eventually re-normalize J but this has only a minor effect on μh.26 To introduce a first phonon quantum correction to transient localization theory, we have recomputed μh assuming that the σij on J is only due to low frequency modes, which agrees with another recent theoretical approach on naphthalene.39
For the BTBT-based materials, our simulated μhs qualitatively match the measured values (see μtheory, and μexp in Table 1), in that the theory appropriately ranks the materials in order of highest μh. The temperature dependence of μh for all materials follows a power law, μh ∼ T−a, where the exponent ranges from 0.835–1.19 (see ESI,† Section S6), which agrees with literature reports of “band-like” mobility in OSCs.40,41 The extremely low μh value for BTBT is likely due to extrinsic defects within the crystal that are very effective in disrupting the quasi one-dimensional transport. Since there is no conclusive evidence for band-like transport for BTBT in the literature, it is likely that the reported BTBT μh in Table 1 is more indicative of thermally activated transport, and thus not applicable to our theoretical predictions. However, the predictions for c8-BTBT and m8-BTBT demonstrate excellent quantitative agreement with experimental results.
The simulated μhs for the substituted acene materials are also captured accurately, again properly ranking the materials. Assuming that each experimental μh is comprised of multiple measurements across multiple domains that are isotropically averaged over the high mobility plane, then the experimental mobility should match the average theoretical μh. However, each average theoretical μh overpredicts the experimental μh, which potentially provides information on the impact of defects on μh. Such information might be valuable as future studies can investigate how different types of defects impact μh, which can guide how to optimize processing conditions for OSC thin films. While the experimental μhs in Table 1 represent reasonable averages found in the literature,26 there are examples of μh measurements for materials with potentially fewer extrinsic defects than average. For TIPS-PN, ref. 38 show that devices can exceed 1 cm2 V−1 s−1 (reaching as high as 1.8 cm2 V−1 s−1). For TIPS-ADT, ref. 30 show devices with μh as high as 0.6 cm2 V−1 s−1. In both cases, our TLT results are in quantitative agreement.
Amongst the substituted acenes studied, TES-ADT demonstrates the poorest agreement between μtheory and μexp. The reason is because we only consider one chemical polymorph, and one crystal polymorph, but we know that this material exhibits multiple chemical/crystalline polymorphs.28,29 Such polymorphs don't disrupt the crystal packing as the ADT core has the same orientation, only the sulfurs are positioned at different positions on the backbone. Inclusion of these polymorphs in mobility simulations would introduce disorder in the equilibrium transfer integrals (acting as extrinsic defects), which would decrease the simulated mobility. We expect that TES-ADT films have a greater density of extrinsic defects than c8-BTBT, which would result in a greater disparity between simulated and experimental μhs.
In addition to TLT, we compute the effective masses of each material, to demonstrate how band theory (without nonlocal EP coupling) predicts μh. Without nonlocal EP coupling to scatter charge carriers, the observed μh is determined by the impurity scattering in the film. Assuming a constant impurity scattering relaxation time for all materials, the values for μh predicted by band theory correlate with the inverse of the effective mass. In Table 1, we show that the inverse effective mass (1/m*) for c8-BTBT is equal to m8-BTBT, 1/m* for TIPS-PN is nearly equivalent to TIPS-ADT, and that 1/m* for TES-ADT is more than 1.5-fold greater than any other material. None of these trends are consistent with experimental values. In all cases, TLT improves on these predictions through the inclusion of nonlocal EP coupling, demonstrating the importance of properly accounting for nonlocal EP coupling in any theory that computes μh.
One of the goals of a detailed description of the system Hamiltonian is the identification of useful paths to improve OSC performance. We can, for example, investigate the chemical origin of the different mobility (and dynamic disorder) in m8- and c8-BTBT to establish principles for designing new materials. Low dynamic disorder can be associated with a smaller gradient of the transfer integral (∇Jij), but this is not the case for these molecules, as shown in Fig. 5 where this quantity is represented graphically for the relevant molecular pair. We conclude (with further details given in ESI,† Section S5) that it is the phonon mode polarization that causes the μh difference between these two materials. One can envision a materials discovery approach where molecular arrangement with the smallest ∇Jij are identified first and phonon calculations are used to verify if dynamic disorder will be sufficiently small.
Fig. 5 Graphical representation of the transfer integral gradient in transport directions A and B for m8-BTBT and c8-BTBT. |
Although our work is most related to charge transport in organic transistor materials, vibrational dynamics are also important to other branches of OSC technology such as exciton separation in OPVs, dynamical processes in triplet emitters, and heat/charge transport in low temperature thermoelectrics. In addition to OSC technologies, our methodology can be applied to diverse fields ranging from topics in catalysis such as reactant binding isomerism/kinetics in MOFs and biological enzymes, to transport mechanisms in energy storage technologies such as fuel cells, CO2 electrolyzers, and artificial photosynthesis. More generally, we have demonstrated a way of integrating phonon measurements and modeling to drive discovery in any field that requires low energy dynamical information.
6,13-Bis(triisopropylsilylethynyl)pentacene (TIPS-PN) was dissolved in boiling acetone, then filtered and re-boiled before covering and allowing solvent to evaporate over several days in a chemical fume hood. The crystals were collected by vacuum filtration, then washed with 0 °C acetone. 5,11-Bis(triethylsilylethynyl)anthradithiophene (TES-ADT) was dissolved in boiling hexanes, then covered, and allowed the solvent to evaporate over several days in a chemical fume hood. The crystals were collected after complete evaporation of the solvent. 5,11-Bis(triisopropylsilylethynyl)anthradithiophene (TIPS-ADT) was prepared and purified as described in ref. 42.
The inelastic neutron scattering spectra were obtained from the VISION spectrometer24,25 at the Spallation Neutron Source at Oak Ridge National Laboratory. The samples were made up of crystalline powders with crystal dimensions of 10s–100s of μm. Each data set is a powder spectrum with no orientational effects. The samples were cooled down to 5 K before collected scattering data to minimize the effect of the Debye–Waller factor on the final results. Measured neutron data for all 6 samples is publicly available in an archive.43
Phonon calculations were performed using plane-wave density functional theory using the Vienna Ab initio Software Package (VASP)44,45 with projector augmented-wave pseudopotentials.46 Unit cell simulations were performed on BTBT to find the appropriate functional, basis size, k-point sampling density, and convergence tolerances (see ESI,† Section S4). The optPBE van der Waals density functional was chosen with a basis size of 520 eV. The appropriate electronic k-point sampling density was 0.5 Å−1. The electronic energy was converged to 10−8 eV, and the interatomic forces to less than 2 × 10−3 eV Å−1 to consider the molecular geometry fully optimized. The unit cell simulations were performed using the MedeA 2.20.2 software suite with VASP 5.4.1 built in. To perform supercell simulations, we prepared the input files using MedeA, then simulated the supercells on the Edison and Cori supercomputers provided by the National Energy Resource Scientific Computing Center (NERSC). Once the supercells were optimized, we used Phonopy47 to calculate the vibrational modes and energies over a specified phonon k-point sampling grid (different for each material) using the finite difference method in which the atoms were perturbed by 0.01 Å. The vibrational modes, and energies were converted to simulated INS spectra using the oCLIMAX tool provided by Oak Ridge National Laboratory.48 We sampled the phonon Brillouin zone densely enough to achieve convergence in the simulated INS spectrum for each material.
Electron–phonon coupling parameters for each phonon mode of every material are computed via numerical differentiation. The variances of all transfer integrals (Jij) for a particular material are computed according to eqn (3) in the text. These variances are used as an input for hole mobility computation using a tight-binding Hamiltonian approach called transient localization theory. A Hamiltonian is constructed for a finite 2D lattice of 5000 localized sites in which the values for Jij are randomly selected from a Gaussian distribution centered on the equilibrium value of Jij with the variance defined above. The Hamiltonian is diagonalized, and the eigenvectors and eigenvalues are used to compute a localization length for the given snapshot of disorder in the system. The Hamiltonian is reconstructed and diagonalized 25 times to achieve convergence in the computed localization lengths. Full details of transient localization theory is available in ESI,† Section S6.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8mh01069b |
This journal is © The Royal Society of Chemistry 2019 |