Mechanoelectric Response of Single-Crystal Rubrene from Ab Initio Molecular Dynamics

A robust understanding of the mechanoelectric response of organic semiconductors is crucial for the development of materials for flexible electronics. In particular, the prospect of using external mechanical strain to induce a controlled modulation in the charge mobility of the material is appealing. Here we develop an accurate computational protocol for the prediction of the mechanical strain dependence of charge mobility. Ab initio molecular dynamics simulations with a van der Waals density functional are carried out to quantify the off-diagonal electronic disorder in the system as a function of strain by the explicit calculation of the thermal distributions of electronic coupling matrix elements. The approach is applied to a representative molecular organic semiconductor, single-crystal rubrene. We find that charge mobility along the high-mobility direction a⃗ increases with compressive strain, as one might expect. However, the increase is larger when compressive strain is applied in the perpendicular direction than in the parallel direction with respect to a⃗, in agreement with experimental reports. We show that this seemingly counterintuitive result is a consequence of a significantly greater suppression of electronic coupling fluctuations in the range of 50–150 cm–1, when strain is applied in the perpendicular direction. Thus our study highlights the importance of considering off-diagonal electron–phonon coupling in understanding the mechanoelectric response of organic semiconducting crystals. The computational approach developed here is well suited for the accurate prediction of strain–charge mobility relations and should provide a useful tool for the emerging field of molecular strain engineering.

O ne of the major advantages of organic semiconductors (OSs) compared with inorganic materials is their mechanical flexibility. 1 Whereas flexible displays based on organic light-emitting diodes (OLEDs) are already a multibillion dollar industry, 2 one can imagine a great deal more in terms of applicationswearable sensors, artificial skin, and soft robotics, to name a few. In some of these applications, one requires robustness of the electronic properties to mechanical strain, for example, in flexible displays, whereas in others, a strong response is desired, for example, in sensors. As such, a detailed theoretical understanding of the mechanoelectric response in organic semiconductors is indispensable and very timely.
A fundamental, molecular-scale understanding of the mechanoelectric response of OS is just about to emerge. To date, the best studied system in this respect is crystalline rubrene 3−7 (Figure 1), but even for this system, experimental measurements have given rather controversial results. To facilitate the discussion, it is useful to define the mobility− strain enhancement factor, g ij which measures the fractional change in charge mobility, μ i , along direction i for 1% strain along direction j, ε j . (A negative g value indicates that compression, that is, negative strain, leads to an increase in mobility.) Morf et al. found that in tension, the hole mobility of rubrene slightly decreases, as one might expect, g aa ≃ −4, g ab ≃ −9, 5 and similar results were reported by Matta et al., g aa = −6, g ab = −21. 6 It came perhaps as a surprise that the mobility change is larger when strain is applied in the perpendicular direction than in the parallel direction with respect to charge flow (i.e., |g ab | > |g aa |). Most recently, Choi et al. reported radically different results, in terms of both the magnitude and the dependence on the strain direction, g aa ≃ −70 to −110, g ab ≃ 0. 7 In the latter study, particular care was taken to address some of the unaccounted factors in previous experiments. Four-probe measurements were used, taking care of the possible strain dependence of the contact resistance 8 (whereas previous experiments used a twoprobe setup). In addition to the standard field-effect transistor (FET) mobility measurement, the Hall effect was used to obtain the mobility, which has the advantage of probing intrinsic charge carriers. 9 A few theoretical studies on the mechanoelectric properties of rubrene have been carried out as well, 10−12 broadly agreeing with the experimental studies of Morf et al. and Matta et al. but in contrast with those of Choi et al. Gali et al. used standard band theory, which is known to be problematic for the estimation of charge mobility in molecular organic crystals, 10 as is small polaron hopping theory. 13 Ruggiero 14,15 and treated the off-diagonal electron−phonon coupling and the lattice dynamics in the linear and harmonic approximation, respectively. 11,12,16,17 Moreover, in the latter study, a semiempirical electronic structure method was used for the lattice dynamics. Because the electronic parameters that govern charge mobilities are very sensitive to fine details of the intermolecular structure and dynamics, the above approximations could tip the balance in favor of one set of experimental results 5,6 or the other. 7 Here we address these issues by employing an accurate but computationally more expensive technique for the sampling of intermolecular dynamics governing hole mobility in strained and unstrained rubrene without assuming a harmonic approximation for lattice dynamics or a linear approximation for electron−phonon coupling. Moreover, we employ ab initio molecular dynamics with a van der Waals density functional to sample the room-temperature thermal distribution of electronic couplings between rubrene molecules in the strained and unstrained crystal, and we use them in the framework of TLT to compute charge mobilities. The latter has been shown to give charge mobilities that correlate well with the results of the explicit time propagation of the charge-carrier wave function. 18−21 Our computations support the experimental results of Matta et al., and they afford a molecular-scale explanation of the peculiar strain dependence of the mobility (|g ab | > |g aa |).
We first address the question of which density functional is most appropriate to describe the structure of unstrained and strained rubrene. To this end, we have investigated a number of van der Waals 22−27 and dispersion-corrected 28 density functionals and have compared the optimized lattice constants (at 0 K) to the experimental structure at the lowest temperature reported, 100 K 29 ( Figure 2). Assuming that the remaining effects of the thermal motion and zero point energy are small, the optPBE-vdW functional 22 appears to best reproduce the experimental structure. Hence, we find that the optPBE-vdW functional gives the best overall performance, and this functional is used in subsequent structural optimizations and ab initio molecular dynamics.
The structure of unstrained rubrene was obtained by optimizing a 3 × 1 × 1 supercell (12 molecules, 840 atoms) and applying fixed lattice vectors taken from the experimentally determined structure at 294 K. 32 Finite temperature-related volumetric expansion relative to the optimized 0 K structure and zero-point motion effects on the cell dimensions are thus 325. This method is referred to as sPOD/PBE. The scaling factor was obtained from a best fit to ab initio reference values for the HAB11 database of organic dimers, 31 resulting in a very small mean relative error of only 3.5%. See the Supporting Information (SI) for further details. The strain dependence of electronic coupling is shown in Figure 3. We focus on the electronic coupling along the highmobility direction a⃗ , J a (Figure 3a), and note that similar results apply to coupling along the b ⃗ direction, J b (Figure 3b). We go to larger values of strain in the case of compression (negative strain) because we are primarily interested in exploring how far we can increase the mobility. Strain values of up to ∼1% and compression along both directions a⃗ (blue) and b ⃗ (orange) result in a linear increase in J a . We notice that the coupling increases slightly more strongly when strain is applied along b ⃗ , in line with previous studies. 6,10,11 At larger strain values, this difference becomes very pronounced. Whereas J a keeps increasing for strains along b ⃗ > ∼6%, it saturates to a plateau value for strains along a⃗ > ∼6%. We note that in the experiment, rubrene samples under strains of >∼0.4% are prone to forming cracks and become unstable on macroscopic time scales; however, from a theoretical perspective, it is of interest to explore the response at higher strain values because the material might sustain these when applied for very short durations of time.
The dependence of the coupling enhancement on the strain direction can be rationalized by analyzing the structural response to strain in terms of the center of mass distance, d COM , and the tilt angle, θ, of the rubrene molecules and their effect on the overlap between the highest occupied molecular orbitals (HOMO) that mediates the coupling (Figure 3c). Strain along a⃗ leads to a linear decrease in d COM (Figure 3d  The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter the kind seen above does not occur. On the contrary, the "θtwist" motion ensures progressively constructive interference and increasing orbital overlap as θ decreases, resulting in progressively increasing electronic coupling. Whereas the electronic couplings calculated above for 0 K give first clues with regard to the effect of mechanical strain, it is the thermal distribution of the couplings that ultimately governs the charge mobility in OS. We used DFT(optPBE)-MD to sample the thermal motion of the rubrene crystal at zero strain, −0.8% strain along a⃗ , and −0.8% strain along b ⃗ at room temperature. For each type of electronic coupling, J a and J b , we considered two independent dimers in the supercell; that is, we calculated four electronic coupling time series, J a,1 (t), J a,2 (t), J b,1 (t), and J b,2 (t), for each value of strain. Two of these time series, J a,1 (t) and J b,1 (t), are shown in Figure 4 for each value of strain, and numerical results are summarized in Table  1.
Importantly, we observe a significantly larger decrease in the root-mean-square fluctuations of the couplings, σ a = ⟨(J a − ⟨J a ⟩) 2 ⟩ 1/2 , corresponding to a decrease in off-diagonal disorder, for compression along b ⃗ (13% decrease) compared with compression along a⃗ (6% decrease); see Table 1. To gain further insight into the modes responsible for the dynamics, we calculate the spectral density function of the electronic coupling time series from the cosine transformation of the autocorrelation function. The running integral of the spectral density yields the cumulative disorder, including all frequencies up to ω, σ α (ω), allowing us to quantify the relative contribution of each mode Ä where S α (ω) is the spectral density for time series J α (α = a, b) and β = (k B T) −1 . (See the SI for further details.) Including all frequencies ω → ∞ returns the root-mean-square fluctuation of the time series. Figure 5 displays plots of the spectral densities and corresponding σ(ω) functions for J a (Figure 5a) and J b (Figure 5b) over all values of strain, using the same J a/b time series considered in Figure 4a−d. We observe that in the case of strain along b ⃗ , there is a noticeable reduction in the spectral density amplitude for J a in the frequency range 50− 150 cm −1 , ultimately leading to a smaller overall root-meansquare fluctuation, as evidenced in Table 1.
Using the DFT-MD electronic coupling distributions, we calculate the hole mobility for unstrained and strained rubrene    , where λ = 0.152 eV is the inner-sphere reorganization energy. In method 3, the electronic coupling distributions are renormalized with an effective, λ-dependent band renormalization factor. 33,34 See the SI for further details and Table S1 for a summary of numerical values. Agreement with the currently accepted experimental mobility of unstrained, highly pure singlecrystalline rubrene, ∼15 cm 2 V −1 s −1 , 35 is reasonable, within a factor of 3 for all three methods.
The mobility−strain enhancement factors, as defined in eq 1, are summarized in Table 2. For all three methods, we obtain g aa , g ab < 0, and |g aa | < |g ab |. Importantly, the values are in good agreement with the experimental measurements of ref 6, with deviations of g values of less than a factor of 2. Our key result, |g aa | < |g ab |, is due to the significant suppression of electronic coupling fluctuations (off-diagonal electron−phonon coupling) in the range of 50 to 150 cm −1 as well as a small increase in the mean electronic couplings when strain is applied along b ⃗ . (See the discussion above and Table 1.) The suppression of J a coupling fluctuations is not as effective when strain is applied along a⃗ (see Figure 5a), and we see only a negligible increase in the mean electronic couplings, which explains the observed sensitivity of the mobility enhancement on the strain direction.
The two previous computational studies of Landi et al. and Ruggiero et al. arrived at similar results, in particular, that |g aa | < |g ab |. However, Ruggerio et al. did not find a similar suppression of electronic coupling fluctuations with compressive strain, and, surprisingly, g aa was reported to be slightly positive. Their calculations differed in a number of important aspects from our calculations; in particular, they used the linear electron−phonon and harmonic approximations (not used here), and there are differences with regard to the unit-cell dimensions (DFT minimum vs experiment), the method for the electronic coupling calculations, and the DFT functional used, which can all contribute to the differences observed.
Here we investigated the accuracy of the linear approximation for electron−phonon coupling in some detail; see the SI. Our calculations indicate that this approximation can give errors in mean electronic couplings that are on the same order of magnitude as the (small) effect of 0.8% compressive strain. Hence, this approximation can become problematic in situations where high accuracy is required like in the present application. Unfortunately, a consistent investigation of the accuracy of the harmonic approximation could not be carried out because phonon calculations for supercells as large as the ones used in the current ab initio MD simulations remain impractical.
According to our calculations and the experiments of refs 5 and 6, the mobility enhancement at strain values that can be realized in the experiment without causing plastic deformation (ε < 0.4%) is rather modest (<17%). From a theoretical perspective, it is of interest to explore what strain values would be required to achieve more significant mobility enhancements. Because the simple linear relation eq 1 is no longer valid for large compressions, we take the computed electronic coupling values displayed in Figure 3 for −12% compression and the root-mean-square fluctuations obtained from DFT-MD at  g values presented are obtained by using the overall distributions of electronic couplings (see Table 1) as input to TLT calculations. The error bars are equal to half the difference between the g values calculated using the electronic coupling distributions of one set of dimers (J a,1 , J b,1 ) and the other (J a,2 , J b,2 ). b No diagonal electron−phonon coupling; that is, all site energies are set to zero. c Site energy fluctuations from the Gaussian distribution corresponding to λ = 0.152 eV. d Diagonal electron−phonon coupling accounted for via band renormalization.
The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter −0.8% compression. Because we expect that electronic coupling fluctuations will further decrease with increasing compression, this is most likely an underestimate, and we therefore consider the values obtained as lower limits. We obtain μ a > 120 cm 2 V −1 s −1 and μ a > 48 cm 2 V −1 s −1 , corresponding to a 3.5-fold and 1.5-fold mobility enhancement for strain along b ⃗ and a⃗ , respectively (band-renormalized, TLT method 3). These extremely high mobility values are unfortunately not experimentally accessible because structures under such high compressions will quickly undergo mechanical failure.
In summary, we have used DFT-MD simulations with the optPBE-vdW functional and TLT to calculate the strain dependence of mobility in single-crystalline rubrene. Our results are in good agreement with the experimental studies of refs 5 and 6. In particular, they explain the somewhat counterintuitive observation that the mobility enhancement in rubrene is larger when strain is applied in the perpendicular direction than in the parallel direction with respect to the electron flow. However, our study (and the one of Landi et al.) disagrees with the most recent and very carefully conducted experimental work of Choi et al. It is possible that interfacial sample/substrate effects particular to that specific experiment or remaining structural defects, not accounted for in the computational models, are responsible for the discrepancy. However, further experimental studies will be necessary to give a more conclusive answer.
The intrinsic mechanoelectric response in pure, singlecrystalline rubrene is found to be relatively modest, which is an advantage for flexible electronics applications where the electronic properties should remain robust with respect to mechanical strain. However, this also means that applying gentle external compressive strain as a means to boost the charge mobility is not very effective for this material. Yet, recent theoretical evidence suggests that the mechanoelectric response may strongly depend on the organic molecule under consideration, 12 more specifically, on the nodal shape and the relative orientation of the charge-mediating molecular frontier orbitals in the crystal. Hence it might be possible in the future to develop organic semiconducting materials that show either weak or strong dependence on the external strain, as desired for a given electronic application.
Details of the electronic structure calculations, structural optimizations, molecular dynamics settings, electronic coupling calculations using the POD method, assessment of the linear approximation for non-local electronphonon coupling, equations defining the cumulative disorder, overview of TLT and discussion of the different flavors used to account for local electron− phonon contributions, and a