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.

One 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 ), 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, gwhich measures the fractional change in charge mobility, μ, along direction i for 1% strain along direction 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 ≃ −4, g ≃ −9,[5] and similar results were reported by Matta et al., g = −6, g = −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| > |g|). Most recently, Choi et al. reported radically different results, in terms of both the magnitude and the dependence on the strain direction, g ≃ −70 to −110, g ≃ 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 two-probe 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]

Figure 1

Slice of the conductive a–b plane of crystalline rubrene. The unit-cell and electronic coupling between molecules in the high mobility direction a⃗, J, and the low mobility direction b⃗, J, are indicated. Electronic coupling in the direction orthogonal to the a–b plane is negligibly small.

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 et al. and Landi et al. used the more suitable transient localization theory (TLT)[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| > |g|).

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 ). 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.

Figure 2

Performance of density functionals in predicting the lattice constants of crystalline rubrene. Deviations of (a) the lattice constants and (b) the volume/anisotropy ratio are shown for several van der Waals and dispersion-corrected density functionals (0 K) with respect to experimental values at low temperature (100 K). The optPBE-vdW functional provides the closest match in both plots.

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 included in an empirical manner. Strained structures along a⃗ and b⃗ were obtained by scaling the relevant lattice parameter followed by reoptimization of all atomic positions at fixed strained lattice parameters. Electronic couplings were calculated for dimers taken from the unstrained and strained structures using the projector-operator-based diabatization (POD) technique[30] in combination with the Perdew–Burke–Ernzerhof (PBE) density functional and scaling of the resultant values by a uniform constant of 1.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 . We focus on the electronic coupling along the high-mobility direction a⃗, J (Figure a), and note that similar results apply to coupling along the b⃗ direction, J (Figure b). 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. 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 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.

Figure 3

Strain dependence of electronic coupling in rubrene. (a) J and (b) J at minimum energy geometry for strains between −12 and +3% along a⃗ (blue) and b⃗ (orange). (c) Visualization of the superimposed highest occupied molecular orbitals used for the calculation of J, as obtained from the POD method, for the unstrained case and −12% strain along a⃗ and b⃗. We mark the distance between the center of mass of each molecule (dCOM) and the tilt angle θ. For −12% strain along a⃗, the red “X” marks the onset of overlap between two lobes that contributes destructive interference, resulting in the observed plateau in J versus ε (panel a). (d) Variation in dCOM and (e) θ as a function of strain along a⃗ (blue) and b⃗ (orange).

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, dCOM, 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 c). Strain along a⃗ leads to a linear decrease in dCOM (Figure d), whereas the tilt angle of the molecules, θ, increases only slightly (Figure e). The decrease in dCOM results in the onset of destructive interference between orbital lobes (point marked with “X” in Figure c), counterbalancing any increase in overlap arising elsewhere. The situation is markedly different for strain along b⃗: The tilt angle decreases strongly, whereas dCOM remains virtually unchanged. Destructive interference of 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 and J, we considered two independent dimers in the supercell; that is, we calculated four electronic coupling time series, J(t), J(t), J(t), and J(t), for each value of strain. Two of these time series, J(t) and J(t), are shown in Figure for each value of strain, and numerical results are summarized in Table .

Figure 4

Thermal fluctuations of electronic coupling in unstrained and strained rubrene at 290 K. Time series for (a) J and (c) J with the accumulated average in bold. The unstrained trajectory is shown in dark blue, ε = −0.8% is in cyan, and ε = −0.8% is in red. Accumulated root-mean-square fluctuations for (b) J, σ and (d) J, σ. (e–g) Histograms of the final distributions of J and J. Gaussian fits are shown with dashed lines. The data in panels a–d are for the electronic coupling time series of a single J and J dimer (J, J), whereas the overall distributions, panels e–g, include all calculated electronic couplings (J, J, J, J). A significant decrease in σ (panel b) and a corresponding narrowing of the J distribution (panel g) can be seen for ε = −0.8%.

Table 1

Dependence of Electronic Couplings for Hole Transfer in Rubrene, J and J, on Mechanical Strain, εa

ε	Ja	⟨Ja⟩	σa	Jb	⟨Jb⟩	σb	τTLT
0 (unstrained)	110.2b	107.1 ± 0.9c	28.2 ± 1.3c	–20.0b	–17.0 ± 0.1c	7.7 ± 0.8c	0.068d
–0.8% along a⃗	113.8b	107.3 ± 2.1c	26.6 ± 0.6c	–20.7b	–17.7 ± 0.7c	7.9 ± 0.6c	0.066d
–0.8% along b⃗	115.3b	110.3 ± 0.9c	24.4 ± 0.4c	–21.3b	–17.9 ± 0.1c	7.2 ± 0.4c	0.066d
0, ref [11]	112.4		34.2	–25.9		10.3	0.079
0, ref [12]	149.2		37.8	–20.3		9.2	0.159

All values are in millielectronvolts unless otherwise indicated.

sPOD/PBE electronic coupling for structures obtained by geometry optimization at fixed lattice dimensions using the optPBE-vdW functional. The initial structure and lattice dimensions of unstrained rubrene were taken from the experiment at 294 K.[32]

Thermal average or root-mean-square fluctuations of sPOD/PBE electronic couplings at 290 K, obtained from DFT-MD trajectories generated with the optPBE-vdW density functional. The error bars are equal to half the difference between the values obtained from one set of dimers (Jα,1; α = a, b) and the other (Jα,2; α = a, b).

Time constant (in ps(rad)−1) corresponding to the harmonic mean of the vibrational frequencies, weighted by the average power spectrum of J electronic coupling fluctuations.

Importantly, we observe a significantly larger decrease in the root-mean-square fluctuations of the couplings, σ = ⟨(J – ⟨J⟩)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 . 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 modewhere Sα(ω) is the spectral density for time series Jα (α = a, b) and β = (kBT)−1. (See the SI for further details.) Including all frequencies ω → ∞ returns the root-mean-square fluctuation of the time series. Figure displays plots of the spectral densities and corresponding σ(ω) functions for J (Figure a) and J (Figure b) over all values of strain, using the same J time series considered in Figure a–d. We observe that in the case of strain along b⃗, there is a noticeable reduction in the spectral density amplitude for J in the frequency range 50–150 cm–1, ultimately leading to a smaller overall root-mean-square fluctuation, as evidenced in Table .

Figure 5

Spectral density functions and cumulative disorder, σ(ω) (eq ), for (a) the J time series and (b) the J time series corresponding to the same dimers used in Figure a–d (J(t), J(t)). In panel a, we observe that σ is smaller for strain along b⃗ due to a suppression of the spectral density amplitude in the frequency range 50–150 cm–1.

Using the DFT-MD electronic coupling distributions, we calculate the hole mobility for unstrained and strained rubrene using TLT.[14,15] Three different flavors of TLT are considered to account for the site energy (or diagonal energy) fluctuations. In method 1, all diagonal elements of TLT Hamiltonian are set to zero. In method 2, site energy fluctuations due to intramolecular modes are included by sampling from a Gaussian distribution with width , 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 single-crystalline rubrene, ∼15 cm2 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 , are summarized in Table . For all three methods, we obtain g, g < 0, and |g| < |g|. 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| < |g|, 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 .) The suppression of J coupling fluctuations is not as effective when strain is applied along a⃗ (see Figure a), 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.

Table 2

Strain–Mobility Enhancement Factors Calculated from Equation Using the Transient Localization Theory (TLT) and Results from Other Computational (comp) and Experimental (exp) Studiesa

		gaa	gab	gba	gbb
this workb	comp	–12 ± 2	–41 ± 8	–20 ± 2	–43 ± 7
this workc	comp	–8 ± 8	–25 ± 1	–17 ± 10	–29 ± 1
this workd	comp	–9 ± 1	–32 ± 7	–17 ± 5	–33 ± 5
ref [5]	exp	–4	–9
ref [6]	exp	–6	–21
ref [7]	exp	–70 to −110	0
ref [10]	comp	–9	–12	–15	–18
ref [11]	comp	+4	–12	+4	–13
ref [12]	comp	–8	–16	–10	–16

g values presented are obtained by using the overall distributions of electronic couplings (see Table ) 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, J) and the other (J, J).

No diagonal electron–phonon coupling; that is, all site energies are set to zero.

Site energy fluctuations from the Gaussian distribution corresponding to λ = 0.152 eV.

Diagonal electron–phonon coupling accounted for via band renormalization.

The two previous computational studies of Landi et al. and Ruggiero et al. arrived at similar results, in particular, that |g| < |g|. However, Ruggerio et al. did not find a similar suppression of electronic coupling fluctuations with compressive strain, and, surprisingly, g 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 is no longer valid for large compressions, we take the computed electronic coupling values displayed in Figure for −12% compression and the root-mean-square fluctuations obtained from DFT-MD at −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 μ > 120 cm2 V–1 s–1 and μ > 48 cm2 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, single-crystalline 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.