Impact of Fluoroalkylation on the n-Type Charge Transport of Two Naphthodithiophene Diimide Derivatives

In this work, we investigate two recently synthesized naphthodithiophene diimide (NDTI) derivatives featuring promising n-type charge transport properties. We analyze the charge transport pathways and model charge mobility with the non-adiabatic hopping mechanism using the Marcus-Levich-Jortner rate constant formulation, highlighting the role of fluoroalkylated substitution in α (α-NDTI) and at the imide nitrogen (N-NDTI) position. In contrast with the experimental results, similar charge mobilities are computed for the two derivatives. However, while α-NDTI displays remarkably anisotropic mobilities with an almost one-dimensional directionality, N-NDTI sustains a more isotropic charge percolation pattern. We propose that the strong anisotropic charge transport character of α-NDTI is responsible for the modest measured charge mobility. In addition, when the role of thermally induced transfer integral fluctuations is investigated, the computed electron–phonon couplings for intermolecular sliding modes indicate that dynamic disorder effects are also more detrimental for the charge transport of α-NDTI than N-NDTI. The lower observed mobility of α-NDTI is therefore rationalized in terms of a prominent anisotropic character of the charge percolation pathways, with the additional contribution of dynamic disorder effects.

To the best of our knowledge, there have been limited computational investigations concerning the charge transport properties of NDTI derivatives. In a recent study [31], the effects of electron-withdrawing groups and electron-donating substituents in the α position of NDTI derivatives were computationally investigated to establish design rules for novel semiconductors with both high charge transport properties and environmental stability. In some cases, however, crystal structures were not available, and the study was carried out under the assumption of a one-dimensional (1D) character of the charge transport.
The availability of organic single crystals, because of their general purity, low degree of defects and long-range order, is ideal for analyzing the intrinsic properties of the semiconductor as a function of the molecular arrangement. In this study, two recently synthesized fluoroalkyl-modified NDTI were considered (see Figure 1). One derivative featured fluoroalkyl substituents in the α position of the terminal thiophenes, while the second featured fluoroalkyl substituents on the N position and are hereafter labeled as α-NDTI ( Figure 1a) and N-NDTI (Figure 1b), respectively. The introduction of fluorinated substituents is well known to improve the device stability under ambient conditions, but it also impacts the molecular packing arrangements, thereby leading to different device performances. Indeed, single-crystal OFETs of both fluoroalkyl-modified NDTI showed good electron transport properties, with α-NDTI displaying an average electron mobility of 0.037 cm 2 V −1 s −1 and a highest value of 0.065 cm 2 V −1 s −1 , while N-NDTI scored an average value of 1.27 cm 2 V −1 s −1 (highest value = 1.59 cm 2 V −1 s −1 ) [32]. To model the charge transport properties of the two NDTI derivatives, we employed the non-adiabatic hopping approach recently adopted by us to rationalize the trends in ntype charge transport of various fluorinated and chlorinated PDI derivatives [33,34], as well as the anisotropy of charge transport in fluoroalkylated NDIs [35]. The integrated approach encompasses the quantum chemical (QC) evaluation of the electronic couplings and reorganization energies, followed by the calculation of charge transfer rate constants with different formulations [36]. The latter are then injected in kinetic Monte Carlo (KMC) simulations to model the charge percolations and mobilities. Dynamic disorder effects have been shown to influence the charge mobilities of organic semiconductors by inducing large fluctuations in the transfer integrals driven by slow intermolecular modes [37- To model the charge transport properties of the two NDTI derivatives, we employed the non-adiabatic hopping approach recently adopted by us to rationalize the trends in n-type charge transport of various fluorinated and chlorinated PDI derivatives [33,34], as well as the anisotropy of charge transport in fluoroalkylated NDIs [35]. The integrated approach encompasses the quantum chemical (QC) evaluation of the electronic couplings and reorganization energies, followed by the calculation of charge transfer rate constants with different formulations [36]. The latter are then injected in kinetic Monte Carlo (KMC) simulations to model the charge percolations and mobilities. Dynamic disorder effects have been shown to influence the charge mobilities of organic semiconductors by inducing large fluctuations in the transfer integrals driven by slow intermolecular modes [37][38][39][40][41][42]. In previous investigations, we already considered thermally induced dynamical disorder effects [33][34][35]43], showing that for core unsubstituted PDI derivatives, the remarkable modulation of selected electronic couplings contributes to rationalizing the magnitude and directionality of the measured electron mobilities. Here, we compare the intra-and intermolecular parameters governing charge transport for the two NDTI derivatives, highlighting the impact of charge mobility anisotropy and eventually the role of dynamical disorder effects in affecting the charge mobility.

Results and Discussion
The crystal structures of αand N-NDTI derivatives contain a certain degree of crystallographic disorder. More precisely, there are two possible orientations of thiophene rings inside the crystal, leaving the rest of the molecule identical in the two cases. For both NDTI derivatives, the dominant orientation was labeled ANTI-1, and the other was labeled ANTI-2 ( Figure S1). Both were considered in the calculation of electronic couplings and charge mobilities. To evaluate the effect of different thiophene orientations in the same solid phase, a crystal containing a mixture of the two orientations was also considered for N-NDTI (hereafter labeled mix). After single-molecule geometry optimization, the two N-NDTI structures corresponding to different thiophene orientations converged, as expected, to the same equilibrium geometry. In contrast, two slightly different optimized structures of α-NDTI ANTI-1 and α-NDTI ANTI-2 were obtained because of different conformers of the fluoroalkyl chains.
At their respective equilibrium structures, the lowest unoccupied molecular orbital (LUMO) energies of the two compounds (Table S1 and Figures S2 and S3) were low (−3.9 eV for α-NDTI and −3.8 eV for N-NDTI), in line with previous calculations on NDTI derivatives [24,27,31] and with the requirements necessary to prevent the electron polaron from reducing ambient species [44]. In this regard, α-NDTI should be slightly favored over N-NDTI, since it displays a slightly lower LUMO energy.

Electronic Couplings and Charge Transfer Pathways
The crystal of the α-NDTI display slipped one-dimensional (1D) stacking of its molecules ( Figure 2 and Figures S4 and S5). The identification of all possible near neighbors of a given molecular site and the correspondingly computed charge transfer integrals (see Section 3) showed that only pathway P1, directed along the a axis (Figure 2), carried a significant electronic coupling, whereas all the others displayed negligible values (see Table 1 and Table S3). This is because P1 corresponds to the only molecular dimer in which intermolecular distance and π stacking allow for an efficient overlap between LUMO orbitals. Indeed, for α-NDTI, the combined effect of N and α substituents determined the solid phase packing with a markedly 1D columnar arrangement of the molecules. The electronic coupling associated with path P1, evaluated for ANTI-1 and ANTI-2, displayed slightly varying values in the range of 36-40 meV (see Table 1).  The N-NDTI crystal was characterized by layers of molecules oriented with their conjugated plane almost perpendicular (edge-on) to the b,c plane (see Figures 3 and S6). The molecules in each layer alternate two different site types that determine a twisted cofacial stacking. Notably, the packing between columns of π-stacked molecules, deter-  The N-NDTI crystal was characterized by layers of molecules oriented with their conjugated plane almost perpendicular (edge-on) to the b,c plane (see Figure 3 and Figure S6). The molecules in each layer alternate two different site types that determine a twisted cofacial stacking. Notably, the packing between columns of π-stacked molecules, determined by the reduced number of substituents in this case, afforded smaller distances compared with α-NDTI. Because of the large distance between layers along the a direction, the most relevant paths for charge transport involved molecules in the same crystal layer (parallel to the b,c plane; see Figure 3), while negligible couplings were computed along the a direction. Electronic coupling calculations revealed three main paths that were expected to contribute significantly to charge transport (Table 1 and Figure 3). The P1 path defined a columnar percolation channel along the c axis between pairs of twisted molecules, with a coupling of 35 meV for N-NDTI ANTI-1 and 50 meV for ANTI-2. The relative twisting of the two molecules forming each P1 dimer led to efficient π stacking with a small intermolecular and interplanar distance, resulting in efficient LUMO overlapping and significant electronic coupling.
Path P2 involved molecules belonging to the same site type, (therefore having the same spatial orientation) with a coupling of 10 meV (8 meV for ANTI-2). Because it was directed along the b axis, it would contribute to making charge transport more isotropic. Finally, path P3 defined a zig-zag percolation across two different columns of π-stacked molecules with a minor coupling of 2 meV (1 meV), and it would contribute to seldom jumps to adjacent columns. The larger P1 coupling for ANTI-2 could be rationalized by a more efficient overlap between LUMO orbitals in the thiophene region (see Figure S7). In contrast, the overlap efficiency was slightly reduced for P2 and P3 in ANTI-2.
For the mix crystal of N-NDTI (see Figure S8), the P1 and P2 dimers were composed of one molecule featuring ANTI-1 and one featuring ANTI-2 thiophene orientations. In this case, because of different centers of mass, two different couplings of 46 and 38 meV were obtained for the two components of the P1 path (see Table 1), while the couplings were 9 and 8 meV for P2. In both cases, these were values within the range of those computed for the dimers formed by two identical molecules (ANTI-1 or ANTI-2). pected to contribute significantly to charge transport (Table 1 and Figure 3). The P1 path defined a columnar percolation channel along the c axis between pairs of twisted molecules, with a coupling of 35 meV for N-NDTI ANTI-1 and 50 meV for ANTI-2. The relative twisting of the two molecules forming each P1 dimer led to efficient π stacking with a small intermolecular and interplanar distance, resulting in efficient LUMO overlapping and significant electronic coupling. Path P2 involved molecules belonging to the same site type, (therefore having the same spatial orientation) with a coupling of 10 meV (8 meV for ANTI-2). Because it was directed along the b axis, it would contribute to making charge transport more isotropic. Finally, path P3 defined a zig-zag percolation across two different columns of π-stacked To summarize, fluoroalkyl substitution strongly impacted the solid state organization of the two NDTI derivatives and resulted in charge transfer integrals favoring a 1D charge percolation in α-NDTI and a more isotropic charge transfer network in N-NDTI. Charge transfer rate constants (Equation (4)) are, however, dependent not only on electronic couplings but also on reorganization energies, as is discussed in the next section.

Reorganization Energies
The intramolecular reorganization energies computed using the four-point adiabatic potential (AP) method (see Tables S5 and S6) and from the computed Huang-Rhys (HR) factors S m (Table S7) agreed very closely (0.319 eV for α-NDTI and 0.305 eV for N-NDTI), which implied that the harmonic approximation was acceptable, at least for the higher vibrational frequencies, since their contributions to λ i were dominant. Indeed, as shown in Figure 4, the individual components λ m i (from HR factor S m ; see also Equation (S6)) to λ i revealed, as is generally found for organic semiconductors, a larger contribution from vibrations above 1200 cm −1 , most of which were associated with carbon-carbon bond stretching of the conjugated molecular core, where the largest geometry change occurred (see Figures S9 and S10). Figure 4 also shows that the active frequencies (larger λ m i values) were very similar for the neutral and charged species. Compared with previous investigations on unsubstituted NDTI (λ i = 0.263 eV), [31] our computed λ i values were slightly larger on account of the role of the flexible substituents. Notably, inspection of the computed HR factor S m (see Figures S11 and S12) showed that a large number of low-frequency intramolecular vibrations, mostly associated with the flexible substituents, displayed unexpectedly large S m values that were very likely to be overestimated as a result of the anharmonic character of low-frequency vibrations. More specifically, for α-NDTI and N-NDTI, the highest HR factor was computed for frequencies at 96 cm −1 and 25 cm −1 , respectively, being associated with accordion-like skeletal deformations [45] (see Figure S13).
Because vibrational frequencies below 200 cm −1 can be considered classical degrees of freedom at room temperature, the effective parameters collected in Table 2 and used to evaluate the charge transfer rate constants (see the Marcus-Levich-Jortner (MLJ) formulation in Equation (4)), were determined by retaining only vibrations above 200 cm −1 . For comparison, the effective parameters computed using the entire set of vibrational contributions are collected in Table S7.
frequency intramolecular vibrations, mostly associated with the flexible substituents, displayed unexpectedly large values that were very likely to be overestimated as a result of the anharmonic character of low-frequency vibrations. More specifically, for α-NDTI and N-NDTI, the highest HR factor was computed for frequencies at 96 cm −1 and 25 cm −1 , respectively, being associated with accordion-like skeletal deformations [45] (see Figure  S13). Because vibrational frequencies below 200 cm −1 can be considered classical degrees of freedom at room temperature, the effective parameters collected in Table 2 and used to evaluate the charge transfer rate constants (see the Marcus-Levich-Jortner (MLJ) formulation in Equation (4)), were determined by retaining only vibrations above 200 cm −1 . For comparison, the effective parameters computed using the entire set of vibrational contributions are collected in Table S7.  [46,47]. d Only one value was reported in this case, since both ANTI-1 and ANTI-2 structures, featuring different thiophene orientations in the crystal, converged to the same single-molecule equilibrium geometry.
Finally, we note that the magnitudes of the computed reorganization energies, compared with the largest computed transfer integrals (never exceeding 50 meV), implied a substantial reliability of the non-adiabatic hopping mechanism ( >> ) [39,48].

Charge Transfer Rate Constants and KMC Simulations
The MLJ rate constants , in the absence of an applied electric field and computed for the charge transfer paths discussed in the previous section, are collected in Table 1 and S3. The computed values for path P1 of both NDTI derivatives were large (10 12 -10 13 s −1 )  Table 2. Intramolecular reorganization energy λ i , effective frequency ω eff , and effective HR factor S eff , with contributions from intramolecular classical vibrations λ classic and from the outer sphere λ o to λ o+classic in Equation (4) for α-NDTI and N-NDTI.  [46,47]. d Only one value was reported in this case, since both ANTI-1 and ANTI-2 structures, featuring different thiophene orientations in the crystal, converged to the same single-molecule equilibrium geometry.
Finally, we note that the magnitudes of the computed reorganization energies, compared with the largest computed transfer integrals (never exceeding 50 meV), implied a substantial reliability of the non-adiabatic hopping mechanism (λ i >> V ij ) [39,48].

Charge Transfer Rate Constants and KMC Simulations
The MLJ rate constants k eT , in the absence of an applied electric field and computed for the charge transfer paths discussed in the previous section, are collected in Table 1 and Table  S3. The computed values for path P1 of both NDTI derivatives were large (10 12 -10 13 s −1 ) and similar, as was expected from the λ i and V ij values. Starting from k eT , we generated the KMC (Brownian) trajectories collected in Figure 5 and Figures S14 and S15. and similar, as was expected from the and values. Starting from , we generated the KMC (Brownian) trajectories collected in Figures 5, S14 and S15. Figure S14 confirms the 1D charge transport of α-NDTI, as expected from the computed transfer integrals and as previously documented for other NDTI derivatives [31]. Similarly, the columnar charge transport resulted in a remarkable anisotropy of the computed time-of-flight (TOF) charge mobilities as depicted in Figure 6a. The KMC trajectories in Figures 5 and S15 show that the charge transport in N-NDTI was essentially bidimensional and parallel to the b,c crystallographic plane as a result of  Figure S14 confirms the 1D charge transport of α-NDTI, as expected from the computed transfer integrals and as previously documented for other NDTI derivatives [31]. Similarly, the columnar charge transport resulted in a remarkable anisotropy of the computed time-of-flight (TOF) charge mobilities as depicted in Figure 6a. Overall, the computed charge mobilities from the zero field (Brownian) and TOF KMC simulations (see Table 3) indicated a similar charge transport efficiency for both NDTI derivatives, with computed values to the order of 1 cm 2 V −1 s −1 . Such similar computed mobilities were, however, in contrast with the experimental observation of a superior charge transport efficiency for N-NDTI vs. α-NDTI.  A possible explanation for the lowest observed mobility of α-NDTI may be related to experimental factors, such as contact resistance limitations, additional structural disorder at the single crystal level, and chemical impurities. However, in seeking intrinsic factors, we suggest here that the remarkably anisotropic charge transport of α-NDTI plays a relevant role. The X-ray diffraction (XRD) measurements indeed showed the c-axis of α-NDTI single crystal standing out of the OFET substrate, while the a,b plane was arranged parallel to the substrate. However, the calculations showed that only along the a direction the mobility was large. In other directions, namely in the a,b plane, which was parallel to the OFET substrate (compare μmax ≈ 1 cm 2 V −1 s −1 and μmin ≈ 10 −6 cm 2 V −1 s −1 in Table 3) and in the direction perpendicular to it, the mobility reduced to almost zero. This might justify the discrepancy between the computed and experimental values, with the latter averaging the OFET mobility in the a,b plane (rather than only along one direction).
XRD measurements of the N-NDTI single crystal displayed the a-axis standing on the substrate, while the b,c plane was arranged parallel to the substrate. In contrast to α-NDTI, the KMC simulations showed a similar charge transport efficiency in every direc- The KMC trajectories in Figure 5 and Figure S15 show that the charge transport in N-NDTI was essentially bidimensional and parallel to the b,c crystallographic plane as a result of non-negligible transfer integrals for paths P1, P2, and P3 ( Figure 3). The b,c plane almost coincided with the Y,Z Cartesian plane ( Figure 5), on which the trajectories' projections were therefore the largest. In comparing N-NDTI ANTI-1 and ANTI-2, their charge transport trajectories showed similar displacements along Z (almost coincident with c). For ANTI-1, the extension of the charge displacement in the two directions (Y,Z) was quite similar, although the P1 rate constant (along c ≈ Z) was larger than that of P2 (along b = Y) (Table 1). However, because the intermolecular distance for path P2 was larger than that of P1, the charge hopping through P2 covered a larger distance along b (Y) compared with that covered through P1 along c (≈Z). This interplay between rate constants and intermolecular distances spanned by the charge led to a modest anisotropy overall in the b,c (or Y,Z) plane for N-NDTI ANTI-1 and a more pronounced one for ANTI-2. For the latter, indeed, the rate constant for jump P1 largely dominated over that of P2. Similar differences could be seen in the computed anisotropies of the TOF charge mobilities (see Figure 6b). Focusing on the three crystal structures of N-NDTI investigated, namely ANTI-1, ANTI-2, and the mix, the first (ANTI-1, purple in Figure 6b) showed slightly lower mobilities than the others and the lowest anisotropy. The second (ANTI-2 in green) displayed larger and more anisotropic mobilities, and finally, the last one (mix in cyan) was characterized by an intermediate situation.
Overall, the computed charge mobilities from the zero field (Brownian) and TOF KMC simulations (see Table 3) indicated a similar charge transport efficiency for both NDTI derivatives, with computed values to the order of 1 cm 2 V −1 s −1 . Such similar computed mobilities were, however, in contrast with the experimental observation of a superior charge transport efficiency for N-NDTI vs. α-NDTI.
A possible explanation for the lowest observed mobility of α-NDTI may be related to experimental factors, such as contact resistance limitations, additional structural disorder at the single crystal level, and chemical impurities. However, in seeking intrinsic factors, we suggest here that the remarkably anisotropic charge transport of α-NDTI plays a relevant role. The X-ray diffraction (XRD) measurements indeed showed the c-axis of α-NDTI single crystal standing out of the OFET substrate, while the a,b plane was arranged parallel to the substrate. However, the calculations showed that only along the a direction the mobility was large. In other directions, namely in the a,b plane, which was parallel to the OFET substrate (compare µ max ≈ 1 cm 2 V −1 s −1 and µ min ≈ 10 −6 cm 2 V −1 s −1 in Table 3) and in the direction perpendicular to it, the mobility reduced to almost zero. This might justify the discrepancy between the computed and experimental values, with the latter averaging the OFET mobility in the a,b plane (rather than only along one direction). Table 3. Computed Brownian (zero field) and TOF mobilities for α-NDTI and N-NDTI. XRD measurements of the N-NDTI single crystal displayed the a-axis standing on the substrate, while the b,c plane was arranged parallel to the substrate. In contrast to α-NDTI, the KMC simulations showed a similar charge transport efficiency in every direction of the b,c plane (i.e., the substrate plane) of N-NDTI, and therefore, a favorable mobility is always expected in this case. Furthermore, the KMC trajectories ( Figure S15) also showed a non-negligible charge transport along the a direction, which was relevant for the bottom gate's top contact OFET architecture used in the experiments [49].

Zero Field
A second, minor discrepancy between the computed and observed results concerned the slight underestimation of the largest computed TOF mobility of N-NDTI (0.92 cm 2 V −1 s −1 ) compared with the highest experimental mobility (1.56 cm 2 V −1 s −1 ) [32]. This is somewhat unusual, since predicted mobilities are generally overestimated due to the ideal conditions underlying the models used. Several computational and experimental reasons may account for this discrepancy. Although the conditions for the non-adiabatic hopping were satisfied (for both NDTI derivatives), other factors could have limited the quantitative agreement of the computed charge mobilities, among them being the fact that MLJ formulation relies on a single effective parameter and that our KMC simulations did not include the effects of carrier-carrier interactions, a factor that may be relevant for high charge densities [49]. In this frame, charge mobility is a function of several parameters, among which is the charge density, which was not taken into account here [50,51]. We should further mention that from the experimental side, the measured OFET mobilities could be overestimated, as was suggested in previous works [52,53].
Overall, the above considerations may explain the differences between the observed and computed mobilities. Aside from these factors, dynamic disorder may also play a role in affecting the final computed charge mobility, as is discussed in the next section.

Influence of Dynamic Disorder
Extrinsic sources of disorder can affect carrier transport, as has been shown theoretically [54] and experimentally [55]. In addition, there is a general consensus that dynamic disorder is one of the crucial parameters governing charge transport in organic semiconductors. Generally, thermally induced fluctuations of transfer integrals limit charge transport [37,41,56], but for percolation channels characterized by small electronic couplings, they have been shown to open new hopping pathways, thereby increasing the mobility [34,[57][58][59][60]. In recent works, it has been reported that the strength of the dynamic disorder is highly correlated with the gradient of the electronic couplings [56,61] with respect to the phonon modes, a parameter also known as the non-local electron-phonon coupling (or Peierls coupling) [39,40]. Further, the presence of alkyl chains may influence the dynamic disorder via modulation of electron-phonon couplings. In a recent work on alkylated dinaphthothienothiophene (DNTT) [62], it was shown that transfer integral fluctuations are prevalently dominated by a single sliding mode involving long-axis displacements between pairs of molecules. The enhancement of the thermal disorder by a such sliding mode was shown to dramatically depress the charge mobility. For this reason, the sliding mode was identified as a killer phonon.
In previous works, we already investigated the role of dynamic disorder in PDI and NDI derivatives [33][34][35] by running molecular dynamics (MD) simulations followed by QC evaluation of the transfer integrals. The analysis of fluctuations via a Fourier transform of the autocorrelation function of the transfer integrals [63] provided indications of the most active intermolecular phonons. Long and short molecular axis sliding motions were clearly identified as being largely responsible for the transfer integral fluctuations and, interestingly, they were demonstrated to activate (i.e., a phonon-assisted mechanism) additional charge transport channels in a PDI derivative [34].
Seeking possible evidence of phonon-assisted charge transport, we applied a similar strategy to investigate the role of thermally induced disorder in N-NDTI and computed the transfer integral fluctuations for the two charge pathways displaying limited efficiency (i.e., low V ij (P2 and P3; see Table 1)), therefore being more promising for phonon-assisted charge transport activation. The thermal fluctuations, however, were negligible for the two paths, and therefore we completed the investigation including path P1. However, for the latter, due to its large electronic coupling, the dynamic disorder effects may only be detrimental to charge transport, ultimately lowering the mobility. In Figure 7 (top), we show the computed fluctuations of the P1 transfer integrals. The standard deviation (σ) was only 11 meV, about one third of the average electronic coupling (< V ij >= 35 meV; see Table 1), a value smaller than what is typically found (i.e., generally the same size as the electronic couplings). The phonon frequencies, which were most active in determining the P1 transfer integral fluctuations, are shown in Figure 7 (bottom). Low-frequency intermolecular phonons (below 50 cm −1 ) were dominant, and they have been shown to be associated with sliding modes, as was reported in several previous works [34,62,64,65]. However, because the strength of the non-local electron-phonon coupling is related to σ [40,66], its modest value suggests only minor effects of the sliding modes on the N-NDTI charge mobility. Similar to N-NDTI, it can be argued that the sliding phonon modes would also be the primary lattice vibrations determining the P1 transfer integral fluctuations in α-NDTI. only be detrimental to charge transport, ultimately lowering the mobility. In Figure 7 (top), we show the computed fluctuations of the P1 transfer integrals. The standard deviation () was only 11 meV, about one third of the average electronic coupling (〈 〉 = 35 meV; see Table 1), a value smaller than what is typically found (i.e., generally the same size as the electronic couplings). The phonon frequencies, which were most active in determining the P1 transfer integral fluctuations, are shown in Figure 7 (bottom). Low-frequency intermolecular phonons (below 50 cm −1 ) were dominant, and they have been shown to be associated with sliding modes, as was reported in several previous works [34,62,64,65]. However, because the strength of the non-local electron-phonon coupling is related to σ [40,66], its modest value suggests only minor effects of the sliding modes on the N-NDTI charge mobility. Similar to N-NDTI, it can be argued that the sliding phonon modes would also be the primary lattice vibrations determining the P1 transfer integral fluctuations in α-NDTI. Since electronic coupling fluctuations along P1 are expected to depress charge mobility (due to the large value), it is interesting to compare α-NDTI and N-NDTI, seeking evidence for an additional mechanism accounting for the lower experimental mobility of the former. To this end, the modulations of the P1 electronic couplings were computed Since electronic coupling fluctuations along P1 are expected to depress charge mobility (due to the large V ij value), it is interesting to compare α-NDTI and N-NDTI, seeking evidence for an additional mechanism accounting for the lower experimental mobility of the former. To this end, the modulations of the P1 electronic couplings were computed along the sliding coordinates represented in Figure 8a. The computed coupling dependences are collected in Figure 8b, from which the absolute value of the numerically estimated derivatives (associated to the non-local electron-phonon coupling) was 48 meV/Å for α-NDTI and 33 meV/Å for N-NDTI. The larger electron-phonon coupling for α-NDTI may therefore have contributed to reducing its charge transport efficiency, as was similarly shown in a recent work on alkylated DNTT [62]. Although a more detailed assessment of the entire spectrum of non-local electron-phonon couplings should be carried out for a conclusive response, the preliminary results of such a higher dependency of α-NDTI on thermal fluctuations than N-NDTI may further justify its lower observed charge mobility (Table 3), in addition to the factors discussed in the previous section.

Computational Methods and Models
The bulk charge transport was investigated according to the non-adiabatic hopping mechanism. The charge mobilities were evaluated while either assuming a Brownian motion of the charge carrier-that is, in the limit of a zero field and zero concentration-and under the effect of an electric field (see below). According to the hopping model, the relevant charge transfer event is localized on a molecular pair (dimer) formed by two neighboring molecules. The validity of this model depends on the relative magnitude of the electronic coupling and the internal reorganization energy , with required to be considerably smaller than [39,48,67,68]. We verified (vide infra) that this was the case for both molecules investigated here.
The experimentally available crystal structures of the organic semiconductors were used to identify the possible charge pathways by evaluating the vector displacements of the centers of mass for the molecules surrounding a central reference site in the crystal.
While the electronic couplings were evaluated for the molecular dimers extracted from the crystal structures, calculation of the intramolecular reorganization energies was based on the determination of potential energy surfaces for neutral and charged species. The equilibrium structures were obtained from quantum chemical calculations carried out at the B3LYP/6-31G* level of theory. The nature of the stationary points determined by quantum chemical structure optimization was assessed by evaluating the vibrational frequencies at the optimized geometries. The vibrational frequencies were also employed to estimate the vibrational contributions to the intramolecular reorganization energies through the calculation of the HR parameters (see below).

Computational Methods and Models
The bulk charge transport was investigated according to the non-adiabatic hopping mechanism. The charge mobilities were evaluated while either assuming a Brownian motion of the charge carrier-that is, in the limit of a zero field and zero concentrationand under the effect of an electric field (see below). According to the hopping model, the relevant charge transfer event is localized on a molecular pair (dimer) formed by two neighboring molecules. The validity of this model depends on the relative magnitude of the electronic coupling V ij and the internal reorganization energy λ i , with V ij required to be considerably smaller than λ i [39,48,67,68]. We verified (vide infra) that this was the case for both molecules investigated here.
The experimentally available crystal structures of the organic semiconductors were used to identify the possible charge pathways by evaluating the vector displacements of the centers of mass for the molecules surrounding a central reference site in the crystal.
While the electronic couplings were evaluated for the molecular dimers extracted from the crystal structures, calculation of the intramolecular reorganization energies was based on the determination of potential energy surfaces for neutral and charged species. The equilibrium structures were obtained from quantum chemical calculations carried out at the B3LYP/6-31G* level of theory. The nature of the stationary points determined by quantum chemical structure optimization was assessed by evaluating the vibrational frequencies at the optimized geometries. The vibrational frequencies were also employed to estimate the vibrational contributions to the intramolecular reorganization energies through the calculation of the HR parameters (see below).

Reorganization Energy and Electronic Couplings
The reorganization energy is composed of an intramolecular term λ i and an outer sphere contribution λ o , due to the interaction with the surrounding molecules in the crystal. The former was computed either with the AP method, namely via two point determinations from each potential energy surface (neutral and charged states), or via calculations of the HR factor S m [67,69,70], obtained in turn within the harmonic approximation from the dimensionless displacement parameter B m , which is generally employed to evaluate Franck-Condon (FC) vibronic progressions in electronic and photoelectronic spectra (see the SI for further details) [71,72]. The outer sphere reorganization energy λ o was assumed to be 0.01 eV according to recent determinations [46,47]. The performance of different schemes designed for efficient calculation of the intermolecular transfer integrals and site energies for pairs of molecules were reviewed recently [67,[73][74][75]. In the framework of the dimer approach and one-electron approximation, the electronic coupling (charge transfer integrals) V ij =< φ i Ĥ φ j > , where φ i,j are the highest occupied molecular orbital (HOMO) or LUMO orbitals of the two monomers forming the dimer (for p-type and n-type conduction, respectively), was obtained using a fragment orbital approach. Following previous studies [76][77][78][79], the protocol was based on the determination of the matrix H MOB in the monomer orbital basis (MOB), whose off-diagonal elements were the non-orthogonalized electronic couplings: where ε DI M is the diagonal matrix of the eigenvalues associated to the molecular orbitals of the dimer, C DI M_AOB is the matrix of the eigenvectors of the dimer in the atomic orbital basis (AOB), S MON_AOB is the overlap matrix of the monomers in the AOB, and C MON_AOB is the monomer-localized orbitals matrix. C MON_AOB is, therefore, a block diagonal matrix containing the MO coefficients in the AOB from each monomer, with the off-block diagonals set to zero and the superscript t indicating the transpose.
The computed couplings were then transformed in an orthogonalized basis by performing a Löwdin orthogonalization: where S DI M_MOB is the overlap matrix between monomer orbitals, which was obtained as follows from the MO coefficients of the monomer orbitals and the overlap of the atomic orbitals in the dimer configuration S DI M_AOB : For a dimer, this was conducted on the 2 × 2 H MOB matrix including the HOMO (or LUMO) orbitals of the two monomers [73,80]. A detailed discussion of the approximations involved in the fragment orbital approach was reported in a previous work [81].
The electronic couplings were computed at the same level of theory (B3LYP/6-31G*) adopted for geometry optimization and for the evaluation of the reorganization energies. All QC calculations were carried out with the Gaussian16 suite of programs [82].

Charge Transfer Rate Constants and Kinetic Monte Carlo Simulations
Because the quantum nature of the most active modes governing local electron-phonon coupling cannot be neglected, a suitable formulation of the transfer rate constants k eT associated with each hopping event is provided by the MLJ quantum correction of Marcus' equation [83]: In the MLJ formulation of the charge transfer rate constants, the quantum description of the non-classical degrees of freedom is represented by a single effective mode of frequency ω e f f and the associated HR factor S e f f , determined from the set of computed HR factors (see the SI for additional details). Following previous work [33,34], because frequency vibrations below roughly 150-200 cm −1 at room temperature can be described to a good approximation in classical terms, and because of their possible anharmonicity, their contributions were not included in the evaluation of ω e f f . The exceeding classical contributions were summed to the outer sphere reorganization energy λ o , and the total contribution reads as λ o+classic in Equation (4).
The zero field (Brownian) charge mobilities were determined by computing the diffusion coefficient D with a set of KMC simulations [57,66,[84][85][86][87][88] (see the SI for further details). An approximately linear dependence of the mean square displacement (MSD) < [r(t) − r(0)] > 2 as a function of time t was obtained by averaging over the subsets of 1000 KMC trajectories. The diffusion coefficient D was obtained from the fitted linear dependence of MSD by employing Einstein's equation: The charge mobility was then obtained by Einstein-Smoluchowski's equation: In the presence of an electric field, the TOF charge mobilities were obtained with the following relation by applying an electric field F of magnitude 10 5 V/cm: where d f is the distance traveled by the charge along the F direction and τ is the time required to travel the distance d f , which is assumed to be 50 µm. The mobility was averaged over 100 trajectories.

Simulation of Thermally Induced Dynamic Disorder
To assess the importance of thermal motions to the electronic couplings [64,89], we ran MD simulations combined with the QC evaluation of the charge transfer integrals. Molecular dynamics simulations were run on a 5 × 5 × 5 supercell of the crystal unit cell of N-NDTI. The dynamics of the system was studied with periodic boundary conditions employing the MM3 force field [90] and the Tinker code [91]. Since recent studies have shown that low-frequency intermolecular vibrations can modulate the magnitude of the electronic couplings [33][34][35]40,43,61,[92][93][94][95], we froze all the intramolecular degrees of freedom (rigid body approximation) while allowing intermolecular motions. We ran a 100-ps MD simulation in the NVT ensemble and at T = 300 K using a thermal bath. The integration time step was set to 1 fs, and trajectory snapshots were saved every 30 fs. On a selected range of 12 ps, after equilibration, we evaluated the LUMO transfer integral fluctuations for the most relevant charge pathways of N-NDTI and determined their standard deviation.

Conclusions
We investigated the charge transport properties of α-NDTI and N-NDTI, two recently synthesized NDTI derivatives in which fluoroalkyl chains were introduced to improve the n-type character and device stability in ambient conditions.
The computed electronic couplings showed that fluoroalkyl substitution impacted the solid state organization of the two NDTI derivatives; the active charge transfer channels favored a 1D charge percolation in α-NDTI and induced a more isotropic charge transfer network in N-NDTI. For both, the intra-column (P1) electronic coupling dominated. The computed intramolecular reorganization energies were similar for α-NDTI and N-NDTI and slightly larger compared with the unsubstituted NDTI, underscoring the role of the low-frequency vibrations associated with skeletal deformations.
The similar magnitudes of the computed intra-and intermolecular parameters was reflected in the n-type mobilities, predicted to be to the order of 1 cm 2 V −1 s −1 and comparable for α-NDTI and N-NDTI, in contrast with the experimental observation of a superior charge transport efficiency for N-NDTI.
Disregarding the role of impurities, contact resistance limitations, and other experimental factors, this apparent discrepancy could be reconciled by considering the striking difference in the charge mobility anisotropy. While N-NDTI displayed sizable charge mobilities in every direction parallel to the substrate, α-NDTI exhibited a prominent anisotropic character with mobility that reduced to zero in some directions parallel to the substrate. We suggested that such marked anisotropy of α-NDTI could justify the difference between the computed and experimental values, the latter corresponding to the average OFET mobilities.
Seeking additional sources of charge mobility degradation for α-NDTI, we also compared the role of dynamic disorder on the most effective charge transport path (P1) of both derivatives. Although limited to transfer integral fluctuations induced by an intermolecular sliding mode, the computed electron-phonon couplings confirmed a more detrimental effect for α-NDTI compared with N-NDTI.
The lower observed mobility of α-NDTI was therefore rationalized in terms of the strong anisotropic character of the charge percolation pathways, with the additional contribution of dynamic disorder effects.