Chasing the killer phonon mode for the rational design of low disorder, high mobility molecular semiconductors

Molecular vibrations play a critical role in the charge transport properties of weakly van der Waals bonded organic semiconductors. To understand which specific phonon modes contribute most strongly to the electron - phonon coupling and ensuing thermal energetic disorder in some of the most widely studied high mobility molecular semiconductors, we have combined state-of-the-art quantum mechanical simulations of the vibrational modes and the ensuing electron phonon coupling constants with experimental measurements of the low-frequency vibrations using inelastic neutron scattering and terahertz time-domain spectroscopy. In this way we have been able to identify the long-axis sliding motion as a killer phonon mode, which in some molecules contributes more than 80% to the total thermal disorder. Based on this insight, we propose a way to rationalize mobility trends between different materials and derive important molecular design guidelines for new high mobility molecular semiconductors.


Main Text:
The last decade has witnessed drastic improvements of the electronic properties, environmental and operational stability and processibility of organic semiconductors (OSCs) 1,2 .Designing new materials with high carrier mobilities, µ, remains one of the main research objectives to enable faster operation and lower power consumption of circuits and addressing of advanced liquid crystal and organic light-emitting diode (OLED) displays 1,3 .Yet despite exploring a wide range of material systems, charge carrier mobilities in excess of 10 cm 2 V -1 s -1 have only been achieved in very few molecular semiconductors and highly aligned polymers [4][5][6] .At present, despite significant general advances in the comprehension of transport physics, a truly microscopic understanding is still missing and the search for high mobility materials often remains heuristic.
The complex charge transport mechanism taking place in weakly, van-der-Waals bonded molecular crystals has been a topic of intense debate 3,7 .Clear descriptions of transport exist in the limits of very low (µ << 1 cm 2 /Vs) and very high charge carrier mobilities (µ >> 10 cm 2 /Vs), where transport can be described in terms of incoherent hopping 8 of localized charges between sites and Boltzmann band transport of delocalized Bloch electrons, respectively 9 .However, at present for most state-of-the-art molecular semiconductors mobilities are in an intermediate range, where hopping transport is no longer applicable, but mobilities at room temperature remain too small to satisfy one of the main criteria for Boltzmann transport; the mean free path for carrier scattering with lattice phonons has to be significantly larger than the lattice constant.In recent years, detailed theories of transport for this intermediate regime have been developed that take into account the key characteristics of molecular solids: Due to the soft, van der Waals intermolecular bonding the phonon modes in molecular crystals tend to exhibit large vibration amplitudes on the order of 0.1Å 10,11 .These thermal lattice fluctuations cause temporal variations in the transfer integrals and site energies across the molecular lattice.This so-called dynamic disorder limits the ability of charge carriers to form fully delocalized Bloch electron states and imposes a transient localization on the charges leading to a unique transport regime, in which the carriers exhibit both localized and extended characters [12][13][14][15] : On a timescale shorter than the characteristic timescale for intermolecular vibrations the charges are localized by the thermal disorder in the transfer integrals, but on longer timescales the charge carriers undergo diffusive motion driven by the waves of molecular lattice fluctuations.This transient localization scenario is able to reconcile many of the experimental features of transport, including the concomitant observation of delocalized transport signatures, such as an ideal Hall effect and an increase of the coherence factor with decreasing temperature and increasing pressure [16][17][18][19] , and localized signatures, such as radical-cation like charge-induced optical absorptions 20 or strong electron-phonon coupling effects in their band dispersion 20,21 .
It appears also able to reproduce the observed temperature dependence of the mobility 14 , although experimentally it is difficult to establish in many systems that temperature dependent mobility measurements are not affected by extrinsic factors, such as contact resistance, presence of traps or mechanical stress 22 .
Within the transient localization scenario 15 , the charge transport is governed by few key parameters: To achieve the highest carrier mobilities one should minimize the amount of thermal molecular disorder, and at the same time reduce the sensitivity of carrier motion to such disorder.The second condition can be realized with a proper engineering of the intermolecular overlaps: as recently shown, isotropic band structures with approximately equal transfer integrals in the different crystal directions show the highest resilience to dynamic disorder, as they minimise localization phenomena 15 .Whether this requirement is satisfied in a particular system or not, can be established once the crystal structure is known.However, to achieve the condition of low thermal disorder no effective strategies are currently known apart from the generic requirement for stiff structures in order to avoid large-amplitude vibrational modes 11,19 .Little is known to which extent there is scope to influence the thermal disorder by molecular design.The difficulty arises because of the complex unit cell of molecular crystals, which typically contain on the order of 100 atoms.This means that there are hundreds of phonon modes which can in principle couple to the charge motion and contribute to the thermal disorder.Do all these modes make important contributions or are there specific 'killer' modes that are responsible for generating most of the thermal disorder?
In the former case the complexity would be bewildering and there would be little hope for any systematic molecular design to suppress thermal disorder.However, in the latter case we will need to understand the nature of these 'killer' modes in order to derive molecular design guidelines.A detailed, mode-resolved understanding of the vibrational dynamics and electron-phonon coupling, in particular of low-energy intermolecular vibrations, that are more strongly excited than intramolecular modes, is required for the rational design of the next generation of OSCs.
Low frequency spectroscopy techniques have long been used to investigate phonon modes in molecular solids, for example, in the context of pharmaceutical drug development [23][24][25][26][27] or for identifying polymorphism in OSCs [28][29][30] .Quantum mechanical simulations provide powerful tools to interpret the modes observed in experiments; considerable emphasis has been spent over the last decade to achieve the necessary computational accuracy required to fully capture the very weak forces that are responsible for the low-frequency motions in van der Waals bonded molecular crystals and are more challenging to model than in more strongly bonded systems. 24,25.Here we use this state-of-the-art methodology to compare the vibrational dynamics of several important classes of high mobility molecular semiconductors.To experimentally determine the molecular frequencies we use terahertz time-domain spectroscopy (THz-TDS) and inelastic neutron scattering (INS).We model the spectra by fully-periodic quantum mechanical calculations and use the obtained set of normal phonon modes and their frequencies to determine the vibrational amplitudes and the corresponding mode-resolved electron-phonon coupling constants (see SI section S4 and S5).This methodology enables us to assess the contributions that individual modes make to the total thermal disorder s, which can be understood as a sum over the electron phonon couplings of the individual modes weighted with their amplitudes.Similar approaches have very recently been reported by Tu 31 and Harrelson 32 .Here our focus is on applying this mode-resolved analysis to assess whether there are particular phonon modes that dominate the thermal disorder and understand the impact of any such 'killer' modes on carrier mobilities calculated within transient localization theory.
We focus our work on the classical high mobility OSCs, pentacene and rubrene, as well as a group of thienoacenes, which have emerged as a very promising class of high mobility molecules.Brickwall materials (TIPS-Pentacene, diF-TES-ADT,...) have not been investigated here because of their non-ideal dimensionality of charge transport, lower experimental mobilities and larger molecular motions 3 .Figure 1 highlights the full set of studied materials and their respective molecular packing.BTBT and DNTT, as well as their alkylated derivatives C8-BTBT-C8 and C8-DNTT-C8, all exhibit a herringbone packing similar to that of pentacene, while rubrene, the highest mobility (µ) molecule recorded so far (µ @ 20 cm 2 V -1 s -1 ), exhibits a 1D slipped-cofacial packing 6 .Colored arrows indicate the inequivalent pairs of molecules, which give rise to different transfer integral values within the unit cell.The figure also shows an illustrative example of the characteristics of polycrystalline, bottom-gate, top-contact organic field-effect transistors (OFETs) of DNTT and C8-DNTT-C8.The devices exhibit near ideal characteristics and allow a reliable extraction of mobility values.We consistently observe higher mobilities in the alkylated species (µ = 3.4 cm 2 V -1 s -1 ) compared to the non-alkylated one (µ = 1.4 cm 2 V -1 s -1 ) (see SI Section S2 for detailed information).A similar trend has also been reported in the literature for more optimized single-crystal devices with even higher mobilities 33,34 .This is an example of significant mobility differences between related molecular systems, which are commonly reported in the literature, but are currently not well understood.Our main aim in this paper is to investigate to which extent differences such as this can be explained fully in terms of the isotropy of the band structure, as previously suggested 15 , or whether they also reflect important differences in the mode-resolved electron phonon couplings and resulting thermal disorder.The THz-TDS measurements were performed at 80 K, while INS measurements were performed at around 10K at the ISIS TOSCA beamline (Harwell, UK).The two measurements complement each other: INS is able to cover a broad range of frequencies, routinely from 150 cm -1 to the mid-IR (1500 cm -1 ), and has the advantage that the neutron scattering cross section does not involve any mode-specific selection rules.For the interpretation of THz-TDS measurements optical selection rules need to be considered, but they extend the accessible range to very low frequencies (20 cm -1 in our setup).The spectra (Figure 2) reveal the existence of many low-frequency modes in all the investigated solids; the alkylated C8-BTBT-C8 and C8-DNTT-C8 molecules exhibit an even larger number of modes than the unsubstituted ones due to the increased number of atoms in the crystallographic unit cell.Unlike mid-infrared measurements, where characteristic intramolecular modes can often be attributed to specific functional groups and differences in molecular design or crystal packing lead to slight frequency shifts, low-frequency measurements instead probe the complex potential energy landscape and weak intermolecular interactions of the entire molecular crystal.As such, they require a careful theoretical analysis by solid-state density functional theory (DFT).We have hence calculated the harmonic lattice dynamics of the selected molecular crystals with accurate, all-electron DFT simulations.
Following a very recent benchmark study 35 , we employed the D3-corrected 36,37 Perdew-Burke-Ernzerhof 38 (PBE) functional, which has been shown to outperform other existing approaches in the description of the subtle interplay of covalent and non-covalent interactions in complex molecular crystals (for details see SI section S4).One limitation of our technique is that in the THz-TDS experiments and in the simulations, it was possible to measure and compute phonons only at the Brillouin zone centre; our INS measurements at TOSCA sample over the entire Brillouin zone, although at energies > 200 cm -1 the phonon dispersion is expected to be almost negligible 39 .In principle INS is capable of resolving phonon dispersion, but this would have required large single crystals with a weight on the order of 1g, which unfortunately cannot be grown for these materials.The results of our normal mode calculations (dashed lines in Figure 2) adequately reproduce the experimental THz-TDS and INS spectra, both in terms of vibrational frequencies (with a maximum discrepancy to experiments of 7 cm -1 ) and relative peak intensities.The comparison with the experimental data provides an important validation step for the calculated frequencies and atomic displacements of the normal modes that will be used in the following for the assessment of electron-phonon couplings.Our harmonic lattice dynamics simulations for the unsubstituted molecules do not predict any optically-active vibrational mode at frequencies below the window accessible in THz spectroscopy.However, a number of ultra-soft normal modes (frequencies well below 0.5 THz) are found in alkylated systems.A close inspection of the potential associated with these modes reveals a pronounced anharmonic character, which plausibly originates from the mixing of nearly rigid molecular displacements with the, notably anharmonic, torsions of the side chains.Accurate estimates of the vibrational frequencies for these modes were hence obtained through the exact numerical solution of the Schrödinger equation associated with the anharmonic ab initio potential.Anharmonic modes were found to be substantially harder than what predicted within the harmonic approximation (for further details on anharmonic calculations see SI Section S4).All modes above 1 THz were well reproduced by the harmonic approximation and were analysed accordingly.
Having calculated and validated the full spectrum of vibrational modes, we can now establish comprehensive mode-resolved electron-phonon coupling maps for each one of the different materials under study.This represents a necessary step in any theory of phononlimited electronic transport, which can be either formulated in reciprocal space as customary in Bloch-Boltzmann transport 9,40 or in a real-space picture as we do here.While the two approaches are perfectly equivalent, we opt for the latter as it provides a more insightful These maps consist of the off-diagonal electron-phonon couplings  "# =  "  # ) , reflecting the modulations induced by phonon mode qm of the transfer integral Jx between pairs of molecular sites labelled by the index x (see Figure 1), as well as the intra-molecular , corresponding to modulations of the site energy ei of the i-th molecular site.We shall focus primarily on the off-diagonal couplings to lowfrequency modes, as these are the most detrimental for charge transport (the calculated maps for diagonal couplings are presented in SI section S5).
We first focus on the comparison of DNTT with its alkylated version C8-DNTT-C8.
Figures 3(a,e) compare the calculated amplitudes sq of all the zone centre vibrational modes in the 0-400 cm -1 frequency range.In addition to the aforementioned increase in the number of normal modes upon alkylation, the calculations clearly illustrate the decrease in vibrational amplitude with frequency expected for classical harmonic oscillators at thermodynamic equilibrium ( ./ ∝   / ) ). Figures 3(b,f) compare the calculated values of the off-diagonal couplings bcm.(for x = c these couplings are the largest, equivalent maps for the other molecular pairs are reported in SI section S5).In both materials there are a few specific modes across the spectrum that exhibit particularly large couplings and could in principle have a strong impact on charge transport, while most modes exhibit only weak coupling and are likely to be less relevant.From the quantities presented above we can now evaluate how each of the vibrational modes contributes to the energetic disorder, and therefore track which modes are more detrimental to charge transport.The fluctuations induced by each mode on a given transfer integral (for example, Jc) are actually determined by the product σc = σq βc of the mode amplitude σq and its vibrational coupling βc, plotted in Figure 3(c,g).To assess the total vibrational disorder in a given material, one then has to sum up the contributions of all modes to the different intermolecular transfer integrals, and consider the total fluctuation  = 3 4 / +  6 / +  7 / (Figure 3(d,h), see SI section S5 for details).From the data we can clearly identify a few dominant, low-frequency modes which contribute most strongly to the total disorder, i.e. those corresponding to the most prominent bars and the most prominent steps in the cumulant plot (red line).These are the modes which combine large amplitude and large coupling.In both molecules the most detrimental vibration is a low frequency phonon at 23 cm -1 (DNTT) and 24 cm -1 (C8-DNTT-C8), which involves mainly long-axis relative displacements of neighboring molecules in opposite directions (see Figure 3(i) and (j), respectively).Strikingly, this single "sliding mode" is responsible per se for a conspicuous fraction of the total transfer integral fluctuations, ranging from 45% in DNTT to as much as 75% in C8-DNTT-C8; this justifies referring to it as the 'killer' mode.In the case of DNTT, the remaining energetic disorder can mostly be attributed to few other "in-plane modes" at higher frequency, typically below 150 cm -1 , mostly implying displacements of the core atoms within the herringbone layers (see Figure 3(i)).Modes above 150 cm -1 do not contribute significantly owing to their negligible vibrational amplitude.
Based on our ab initio modeling, the effects of alkylation can therefore be understood as follows: (i) new modes emerge in the relevant frequency window; these however mostly involve the alkyl side chains and not the functional cores, and are thus harmless to charge transport; (ii) in-plane modes are shifted to higher frequencies, hence suppressing their impact on disorder and resulting in energy fluctuations becoming dominated by a single sliding mode; (iii) a detailed analysis, presented in the SI, shows that the presence of side chains reduces the amplitude of out-of-plane sliding motion, in accordance with what was reported by Illig et al. 10 .Surprisingly, while all these effects taken together should lead to a reduction of disorder, we find that (iv) the total disorder increases upon alkylation (red lines in Figure 3(d,h)).This counterintuitive behavior can be understood by studying the transfer integral of individual molecular pairs as a function of the relative intermolecular shift along the molecular long axis (sliding coordinate z), as illustrated in Figure 3(k) and SI Figure S18.
Our finding that the long-axis sliding motion is the 'killer' mode indeed justifies a simple analysis with the relevant quantity determining the thermal disorder being the derivatives  "  ) at the equilibrium position ( = 0).While the curves for DNTT and C8-DNTT-C8 in Figure 3(k) are very similar in shape, the equilibrium position in DNTT is very close to a "sweet spot" where  7 () attains an extremum, implying small values of the derivative (this situation is analogous to that of rubrene 43 , which is known to feature a particularly low value of disorder).Alkylation moves the system out of this sweet spot, increasing the coupling by a factor of ~2, ultimately determining the larger energetic disorder with respect to DNTT.
The protocol devised above can now be extended to the whole set of compounds studied in this work.Figure 4 shows the cumulative plots of the transfer integral fluctuations s in units of  = 3 4 / +  6 / +  7 / , the latter setting the overall scale of the electronic bandwidth 15 .
The different systems can be classified into three distinct classes: (1) Alkylated systems (b,d) feature a single vibration dominating the energetic disorder ---the sliding mode --responsible for the large initial step in the cumulative plots of σ. (2) Non substituted systems presenting a pure herringbone packing, such as DNTT, BTBT and pentacene (a,c,f) present a much reduced disorder from the sliding mode, as they all lie close to a sweet spot; the remaining disorder is distributed among a small number of in-plane vibrations (3) Rubrene (e) singles out because of its peculiar molecular shape, where phenyl rings orthogonal to the backbone act to inhibit any possible sliding motion.This is reflected in the absence of a lowfrequency step-like feature in the cumulative disorder plots.Note that the sliding mode for the molecular pairs a in Figure 1(a), which is a zone-boundary phonon, is not included in our analysis based on Gamma-point vibrations; this leads to an underestimation of  4 .Our analysis here shows that a refined treatment that takes into account zone-centre phonons plus the specific zone-boundary phonons that lead to an asymmetric sliding motion of the a molecular pairs would recover most of the energetic disorder.S5), that determine how a given OSC responds to disorder; the mobility in the most conducting direction can increase by up to a factor of 6 when going from a purely one-dimensional to a fully isotropic structure; (b) the magnitude of the energetic disorder, the effect of which has not been considered in detail yet.Importantly, we show below that the dimensionless parameter   ) as used in Figure 4 above provides an effective way to quantify the energetic disorder of a given material and enables us to rationalize differences in mobility between materials.  ) fully determines the effect of disorder on the transport properties, with the mobility expected to scale as  ∝ (   ) ) >? , with  ≈ 2. According to this simple rule, the variations in disorder reported in Figure 4,   ) = 0.2 − 0.5, entails variations of a factor of 6 in the mobility, totally comparable with the effects of anisotropy.It is clear then that anisotropy and disorder compete at the same level, and both aspects need to be considered in order to capture the experimental trends.
We report in Figure 5 the measured mobilities for the different molecules in the set, and compare them with a full calculation of the mobility using transient localization theory 15 , based on the parameters in SI tables S5 and S7.For the experimental values we have selected the most reliable 44 reported mobilities based on robust four-point-probe measurements and further validated by Hall effect measurements (see SI Table S3), which was possible for all compounds except BTBT, due to its high contact resistance.For the theoretical mobilities we show values calculated along the highest mobility direction, taking into account both offdiagonal and diagonal energetic disorder (methods section 7).The data show a good correlation with the experiments, especially considering the large/known/documented sample-to-sample variability characterising organic semiconductors.Quite expectedly the actual experimental values are lower than predicted, because the calculations disregard extrinsic effects that limit the experimentally achievable mobility values, such as contact resistance limitations, chemical impurities and structural defects.6][47] Also, as mentioned above, our simulations do not take into account the dispersion of phonon frequencies away from the Brillouin zone centre, in particular the effect of the zone boundary phonon that involves a sliding motion of the a molecular pair, this is expected to increase the amount of disorder, hence lowering further the theoretically predicted values as recently reported. 31,32,48show theoretical predictions including both off-diagonal and diagonal electron-phonon coupling.Energetic disorder is calculated accounting for all low-frequency modes up to 300 cm -1 .The grey shaded area is a guide to the eye that highlights the correlation between experimental and theoretical values.
We can now understand how the effects of the structure on the inter-molecular vibrations and thermal disorder are ultimately reflected in the mobility of the materials.Based on the transfer integrals alone, alkylation of DNTT should therefore boost the mobility by a factor of 2.5: in DNTT the combination of  " has the wrong sign, while it becomes ideally isotropic in alkylated DNTT 15 .However, the enhancement of thermal disorder by the long-axis sliding 'killer' mode is so strong that it cancels much of this expected gain; as a result alkylated DNTT is predicted to have only about 30% higher mobility than DNTT, which is consistent with experimental observations on single crystals.An analogous compensation is predicted for the other pair of molecules, BTBT and C8-BTBT-C8.Following the same line of reasoning, the exceptionally high mobility value of rubrene can be attributed to its small energetic disorder due to its reduced sensitivity to the long-axis sliding motion, which compensates for its unfavourable anisotropic transfer integrals.In this material an extra boost comes from two secondary factors: its large lattice parameters and, somewhat counterintuitively, a comparatively small value of the bandwidth scale , which ensures a sufficient thermal population of higher energy and more delocalized states in its quasi 1Dbands 49 .In pentacene the unfavorable combination of  " , the higher disorder compared to rubrene or DNTT as well as unfavourable secondary factors all concur to make it the worst performer in the set.
What is striking from these considerations is that none of the state-of-the-art molecules considered here is ideal in the two aspects that have the largest impact on the charge mobility, i.e. none of them combines low energetic disorder with isotropic transfer integrals.If it were possible to rationally design molecules that meet these requirements simultaneously this could enable molecular semiconductors with very high mobilities, potentially exceeding that of rubrene.In this respect, we have shown here that the alkylated compounds, which were selected as very promising candidates owing to their isotropic electronic properties, are in fact hampered by the presence of very strong disorder and that such disorder is carried almost entirely by a single inter-molecular sliding mode.The beneficial reduction of the amplitude of the sliding fluctuations achieved from a stiffer structure is entirely compensated by a concomitant large increase of their coupling to electrons, enhancing their detrimental impact on transport.The origin of this degradation immediately becomes transparent when crossing the information on the sliding mode dependence of transfer integrals (see in Figure 3k) with that on supramolecular packing (see Table S7).In the unsubstituted DNTT and BTBT crystals, the tilt angle of the molecules with respect to the herringbone plane induces an ideal offset in the relative alignment of adjacent molecules at equilibrium: the transfer integrals are near their extremum, where the coupling of electrons to the sliding mode vanishes.The proximity to this local "sweet spot" is lost upon alkylation, as the tilting angle is reduced.
In terms of rational molecular design guidelines for new high mobility materials the general message emerging from the present analysis is that intermolecular shifts at equilibrium determine not only the set of transfer integrals  " but also their dependence on the long-axis sliding, i.e. they ultimately control both ingredients that govern charge transport in 2D herringbone structures.Since such relative shifts are directly affected by the molecular tilting angle, the latter emerges as a natural parameter to play with in the quest for higher mobility semiconductors.This could be achieved by the engineering of side chains, variation of their length, branching or incorporation of specific binding sites.A second important lesson that follows from our insight is that a proper screening of promising candidates in molecular crystal structure databases can be achieved with relatively simple quantum chemistry methods, without the need of performing expensive solid-state phonon calculations.Indeed, evaluating the set of  " and the corresponding  "  ) already contains the key figures needed to rationalize transport.
1578.6 Å 3 ), were Pawley fitted in DASH 54 and then the most probable space group of P21/c determined by use of ExtSym. 55

Inelastic Neutron Scattering (INS)
High-resolution INS measurements on a powder sample of 2g of BTBT, 1g of C8-BTBT-C8, 1g of DNTT, 1g of C8-DNTT-C8, 1g of rubrene and 1g of pentacene respectively, contained in a flat aluminum cell of cross-sectional area 4×4.8 cm 2 were performed at the invertedgeometry neutron spectrometer TOSCA, spanning energy transfers up to 4000 cm -1 , at the ISIS Pulsed Neutron & Muon Source, Rutherford Appleton Laboratory, United Kingdom 57- 59 .The spectral resolution of TOSCA is similar to that of infrared and Raman techniques, approximately 2 % of the energy in question.The sample was studied at a temperature of ~10 K using a dedicated closed-cycle helium refrigerator, so as to minimize the effect of the Debye Waller factor, thus allowing us to fully distinguish the individual peaks. 60

Organic Field-Effect Transistors (OFETs)
Materials: Octadecyltrimethoxysilane (OTMS) was purchased from Fluorochem and used as received (storage under a N2 atmosphere to prevent hydrolysis).Highly doped n-type Si (100) wafers (resistivity < 0.005 Ω cm, Active Business Company GmbH) were used as the substrates for thin film transistors fabrication.A 300 nm SiO2 layer (capacitance Ci = 11.5 nF cm -2 ) was thermally grown onto the Si substrates as a gate dielectric.
Field-Effect Transistor Fabrication: 30 nm films of DNTT and C8-DNTT-C8 were deposited by thermal evaporation (Angstrom Engineering, Inc.) at a rate of 0.2 Å s -1 and 0.1 Å s -1 respectively, under pressure of 10 -7 Torr and substrate rotation of 25 RPM, onto octadecyltrimethoxysilane (OTMS) treated 300 nm thermally grown SiO2 substrates.The monolayer was deposited according to a method published previously 61 .The substrate temperature was maintained at room temperature during the deposition process.
Approximately 40 nm of gold were thermally evaporated through a shadow mask (L = 350 μm, W = 1000 μm) at a rate of 0.5 Å s -1 to complete the devices.
Device Characterization: Transistors were characterized using an Agilent 4155B Semiconductor Parameter Analyser and standard probe station setup at room temperature in a N2-filled glovebox.Device parameters were extracted using the standard calculation techniques.The devices were stored in the dark in a N2-filled glovebox.
Reliability Factor (r) 44 : is designed to show an overall ability of a reported OFET to "deliver on its promise".This parameter is defined as the ratio of the maximum channel conductivity experimentally achieved in a FET to the maximum channel conductivity expected in the correctly functioning transistor with the claimed carrier mobility, .The parameter r simply shows if the reported µ is a good metrics of the actual electrical performance of a FET.The periodic solid-state density functional theory (DFT) simulations were performed using a developmental version of the CRYSTAL17 software package. 62The Grimme D3corrected 36,37 Perdew-Burke-Ernzerhof (PBE) density functional 38 was coupled with the splitvalence triple-zeta 6-311G(2d,2p) basis set 63 for all atoms.All the simulations were initiated from the experimental single-crystal X-ray diffraction (SCXRD) determined structures, 51,[64][65][66][67] with the exception of C8-DNTT-C8 which was determined for the first time for this study using a combination of DFT simulations and refinement against the experimental powder Xray diffraction (PXRD) data.All structures were fully optimized (atomic positions and lattice vectors) with no constraints other than the space group symmetry of the solid.[70] In most cases, the harmonic approximation was sufficient for reproducing the vibrational potential energy curves for the OSC normal modes, however in the case of the alkylated materials there existed significant anharmonicity (deviation from harmonic behavior) in the lowest frequency modes (<30 cm -1 ).In these modes, the explicit anharmonic vibrational energies were determined by scanning the individual normal mode coordinates and determining the explicit vibrational potential energy curves.These values were then fitted with a sixth-order polynomial, and the one-dimensional anharmonic oscillator Schrodinger equation was used to solve for the vibrational energies.The energy convergence criteria were set to ΔE < 10 -8 and 10 -10 hartree for the optimization and vibrational calculations, respectively.

Calculation of electron-phonon coupling constants and charge carrier mobilities
Local (Holstein) and nonlocal (Peierls) linear electron-phonon coupling calculations 7 in periodic solids were performed following the protocol by Girlando et al. 42 with two important improvements: (i) we rely on a full ab initio description of vibrational modes and their coupling to electrons; (ii) we explicitly account for the anharmonicity that strongly affects low-frequency modes of alkylated systems.DFT calculations (PBE0 functional, def2-SVP basis, ORCA v4.0 code) 71 were performed as a function of the displacement along each mode to compute site energies (HOMO level) and transfer integrals for the symmetry-inequivalent molecular pairs.Transfer integrals are computed as one-electron HOMO-HOMO couplings following the dimer projection approach. 72arge transport simulations within the transient localization framework in the relaxation time approximation 14 were performed as described in a previous paper, 15 for a twodimensional tight-binding model for hole states fully parametrized with ab initio inputs for the specific systems at hand.We have considered lattices of 320,000 molecules as to ensure the convergence of the calculated mobility values at 300K.The numerical propagation of carriers' wavefunctions is efficiently performed with the technique introduced in a previous reference. 73

Figure 1 .
Figure 1.Molecular packing of investigated materials and thin-film OFET

Figure 3 .
Figure 3. Electron-phonon coupling maps and identification of most detrimental modes

Figure 4 .
Figure 4. Comparison of total electron-phonon coupling maps between the investigated

Figure 5 .
Figure 5.Comparison between experimentally measured and theoretically predicted Consideration of the cell volume, molecular volume and space group symmetry suggested Z¢ = 0.5.The simulated annealing component of DASH was used to optimise the orientation and conformation of a half-molecule model of C8-DNTT-C8 against the diffraction data (107 reflections), with the molecular centre of symmetry anchored on the origin of the unit cell.A refineable (1 0 0) correction for preferred orientation was included at the structure solution stage.This yielded a favourable c 2 SA/c 2 Pawley ratio of around 3 for the best solution, with a chemically and crystallographically sensible structure.The validity of the structure was confirmed by a Rietveld refinement against data in the range 4-60° 2q using TOPAS 56 and by fixed-cell energy minimisation of the crystal structure using DFT-D.The energy-minimised structure was then used as the basis of a more accurate starting model of a half molecule, and the structure then re-solved and re-refined in space group P21/a to facilitate comparison with other related structures.The final Rietveld refinement, using a spherical harmonics correction for preferred orientation, yield an Rwp value of 3.6 (c.f.Pawley Rwp of 1.8).

6 .
Density Functional Theory Simulations of vibrational modes