Theoretical study of charge-transport and optical properties of organic crystals: 4,5,9,10-pyrenediimides

The shapes of the HOMOs/LUMOs and the relative positions of adjacent molecules decide the size of the intermolecular electronic couplings.


Introduction
Organic semiconductors are of great interest for their intrinsic scientific challenge and potential applications in organic electronic devices such as organic light-emitting diodes (Lin et al., 2016), plastic solar cells (Gao et al., 2017;Skrypnychuk et al., 2016;Yang et al., 2016), organic lasers (Kuehne & Gather, 2016;Zhang et al., 2016), (bio)chemical sensors Haughey et al., 2016) and organic field-effect transistors (OFETs) Sung et al., 2016;Raghuwanshi et al., 2016;Matsushima et al., 2016;Ford et al., 2016). Compared with p-channel materials, the development of high-performance ambient-stable n-channel materials has largely lagged due to the fact that the transport in n-channel conductors is degraded easily by air, and their low electronic affinity hinders efficient injection of electrons into the empty lowest unoccupied molecular orbital (LUMO). In recent years, numerous attempts were made to overcome these difficulties, and some new n-channel semiconductors have been realized via functionalization of Naphthalene diimides (NDIs) (Bé langer-Chabot et al., 2017;Yuan et al., 2016;Purdum et al., 2016;Kobaisi et al., 2016), diketopyrrolopyrrole (DPP) Yao et al., 2016;Yi et al., 2015), perylene diimides (Yue et al., 2014;Liu et al., 2014;Zhang & Zhao, 2012;Wuerthner & Stolte, 2011) and heteroacenes  with electron-withdrawing substituents or alkyl chains. For example, Yuan et al. synthesized difluoro-and tetrafluorosubstituted Naphthalene diimides, and the OFETs based on these fluorinated NDIs exhibited n-channel field-effect character under ambient conditions with a maximum mobility of 0.1 cm 2 V À1 s À1 (Yuan et al., 2016). Klauk and co-workers found that halogen and cyano substituents are ideally suited to tune the DPP dyes' absorption and redox properties into a favorable regime for applications in organic transistors and solar cells, and the alkyl chain substitution could lead to superiorcontacts between DPP molecules within the layers and closer distances between the DPP layers, both of which are favorable for charge-carrier transport (Stolte et al., 2016). In our previous work, the transport parameters of the cyanated bithiophene-functionalized DPP molecule were simulated in the context of the band model and hopping models, and the theoretical intrinsic electron mobility could reach 2.26 cm 2 V À1 s À1 (Huang et al., 2015). More recently, Zhang et al. designed and synthesized a new type of aromatic diimide, pyrene-diimide (PyDI) small molecules (see Fig. 1), which exhibit excellent charge-transport capability and attractive optical properties (Wu et al., 2017). As a potential family of n-type organic semiconductors however, there are very few reports on their structure-function relationship and electronic structural properties up to now. A theoretical investigation of the relationship between the nature of the different alkyl groups and device performance, and a prediction of their intrinsic electron-transfer mobility are instrumental in guiding the molecular design and development of novel organic semiconducting materials, which is also the focus of this work.
In this work, the charge-transport and optical properties of 4,5,9,10-pyrenediimides (C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI, see Fig. 1) were systematically investigated. Firstly, the reorganization energy associated with charge transport was evaluated using the adiabatic potential energy surface (APS) method, and the influence of the bond-parameter variations on the local electron-vibration coupling were discussed in detail, which evaluates well the effects of different alkyl groups on the reorganization energy. Then, the structural and electronic properties of C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI were investigated, and the angular resolution anisotropic mobility for both electron and hole transport were further evaluated using the newly developed simulation methods. In the end, the steady-state absorption and fluorescence spectra were simulated and a new assignment of fluorescence bands in the experiments is confirmed.

Organization energy, ionization potential and electronic affinity
The reorganization energy , associated with the chargetransport process in organic solid materials, can be evaluated in two ways (Zhao & Liang, 2012). The first is the normalmode (NM) analysis method, which partitions the total relaxation energy into contributions from each vibrational mode: where ÁQ i represents the displacement along normal mode Q i between the equilibrium geometries of the neutral and charged molecules; ! i is the corresponding frequency. The other method is the APS method (the four-point approach), in which can be expressed as follows: Here, E and E AE represent the energies of the neutral and cation/anion species in their lowest energy geometries, respectively; E * and E Ã AE are the energies of the neutral and cation/anion species with the geometries of the cation/anion and neutral species, respectively. Our previous studies showed that the APS method is suitable for both flexible and rigid molecules; in comparison, the NM analysis is more appropriate for rigid molecules due to the large deviation of the lattice vibration from the harmonic oscillator model for flexible molecules (Ma et al., 2017a,b). In this work, we selected the APS method to calculate the reorganization energy.
From the APS of neutral/charged species, the vertical ionization potential (VIP), adiabatic ionization potential (AIP), vertical electronic affinity (VEA) and adiabatic electron affinity (AEA) can be calculated as Full-geometry optimizations of the monomer molecules and the reorganization energy calculations are carried out using the B3LYP functional (Lee et al., 1988)   Molecular structures of C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI. Red: O; blue: N; gray: C; white: H.

Electronic coupling
The intermolecular electronic coupling V ij from site i to site j, which describes the overlap of electronic wavefunctions between the donor and acceptor states, can be written as where S ij , J ij and e i(j) represent the spatial overlap, charge transfer integrals and site energies, respectively. These physical quantities can be calculated as follows: e iðjÞ ¼ hÉ iðjÞ jHjÉ iðjÞ i; ð8Þ Here, H is the Kohn-Sham Hamiltonian of the dimer system and É i(j) represents the monomer HOMOs (for hole transport) or LUMOs (for electron transport) with Lö wdin's symmetric transformation, which can be used as the orthogonal basis set for calculation. The calculations of all electronic couplings in different molecular dimers are performed with the PW91/TZVP of density functional theory (DFT) implemented in the Amsterdam density functional (ADF) program (te Velde et al., 2001).

Anisotropic mobility
The anisotropic mobility is an important intrinsic property of charge transport in organic semiconductors, which depends significantly on the specific surface of organic crystals. Herein, we simulated the angle-resolved charge mobility of 4,5,9,10pyrenediimides by means of solving the master equation, which has been described in detail elsewhere (Yin & Lv, 2008;Yin et al., 2012). Charge-transfer (CT) kinetics through a solid material with many possible residence sites can be described by the master equation: where k ij is the CT rate constant from site i to site j in the crystal considering the correction of the electronic field, p i is the charge occupied density on site i, and 1 À p i is the Coulomb penalty factor, which prevents two or more charges at the same time from occupying the same site. If the CT reaches the so-called steady state, dp i /dt = 0, the p i can be obtained by an efficient iterative procedure given a full set of CT constant k ij values. When an external electronic field E is applied to the crystal, the charge will drift accordingly, and the charge mobility m can be determined from the velocity v as the linear response of the motion to the perturbation: whereÊ E is the unit vector of the applied electric field, R ji is the vector from site i to site j and p tol is the total charge population in the investigated supercell. The calculation here is performed using the periodic boundary condition with a supercell of size 3 Â 3 Â 3, and the external electric field E is set to a relatively small value of 1.0 Â 10 À3 V Å À1 .

Electronic spectra
As the accuracy of a TDDFT calculation is strongly dependent on the chromophore family, it is essential that TDDFT results are validated by comparing the experimental data prior to a detailed interpretation. The accuracies of different functionals such as B3LYP, PBE0 and M06-2x were assessed by comparing the predicted wavelengths and intensities of the lowest energy bands with experimental absorption data. Our tests showed that the simulated absorption spectra based on the geometry structures optimized using M06-2x are in better agreement with the experimental UV-vis spectra, which is inconsistent with previous studies (Huang et al., ,b, 2018Ma & Huang, 2016;Yang et al., 2018Yang et al., , 2017. Therefore, the geometry optimizations of t-C 5 -PyDI at the S 0 state and the S 1 state were implemented using the DFT and TDDFT methods at the M06-2x/TZVP (triple-zeta valence quality with one set of polarization functions) level (Treutler & Ahlrichs, 1995). The self-consistent field convergence thresholds of the energy for both the ground-state and excited-state optimization were used as the default settings (10 À6 ). The excited-state Hessian was obtained by numerical differentiation of analytical gradients using central differences and default displacements of 0.02 Bohr. The geometry optimizations were performed without constraints on bond lengths, angles or dihedral angles. All local minima were confirmed by the absence of an imaginary mode in the vibrational analysis calculations. To evaluate the solvent effect, cyclohexane and chloroform were selected as the solvent in the calculations using the conductor-like screening model (COSMO) method (Klamt & Schü ü rmann, 1993).

Reorganization energy
As one of the key parameters influencing the intrinsic charge-transport rates, the reorganization energies evaluated from the four-point approach are collected in Table 1. Comparison of the reorganization energies of C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI shows that the variation of alkyl chain length has little influence on the reorganization energies associated with intermolecular electron-transfer ( e ) and intermolecular hole transfer ( h ), and the introduction of t-butyl groups reduces the h and e values no more than  Table 1 DFT-B3LYP/6-311G** calculated hole-transfer ( h ) and electrontransfer ( e ) reorganization energies of C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI by the APS approach.  Table 1, the e values of C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI are 0.258, 0.257, 0.251 and 0.251 eV, respectively, which are about 0.07-0.09 eV larger than the corresponding h values. This indicates that the pyrenediimide framework undergoes larger geometry relaxations during the electron-transfer process. For the sake of comprehensive analysis of the relationship between the molecular structure and the reorganization energy, we display the bond length alternation in the C 6 -PyDI molecule at oxidation and reduction, as shown in Fig. 2. For bonds 1 to 8, as shown in Fig. 2(a), the geometric relaxation occurs predominantly on electron transfer; in contrast, the smaller geometric changes in these bonds on hole transfer contribute less to h , which is consistent with the observation of a larger values of e than h . For C-C bonds 9 to 17 in the pyrene core, the bond relaxation on hole transfer is more pronounced than those on electron transfer [see Fig. 2(a)], which means that the changes in these C-C bond lengths contribute more to h than to e . Combined with the calculated h and e values, it can be concluded that the large e value is mainly attributed to the geometric relaxations of five-membered imide rings and C O bonds rather than the pyrene core. Moreover, the changes in C-C bond length in the alkyl chains are also depicted in Fig. 2(b). All bond-length variations of the C-C bonds between the neutral and charged forms are less than 0.003 Å , suggesting that the geometry relaxations of C-C bonds accompanied with hole and electron transport have a very small effect on their reorganization energies. This explains the small difference in the reorganization energies of C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI.

Crystal structure and electronic coupling
In the crystal structures of C 5 -PyDI and C 6 -PyDI, the molecules crystallize in a triclinic system and show a typical one-dimensional -stacking motif, which is usually considered to exhibit stronginteractions and good charge-transport properties. In this motif, only face-to-face intermolecular packing modes can be formed. To facilitate the discussion below, this intermolecular packing mode is defined as the P dimer. For C 5 -PyDI and C 6 -PyDI, the calculated electronic couplings for hole and electron transfer (denoted as V h and V e ) and the mass-centered distances r in the P dimer are summarized in Table 2. It can be seen that intermolecular electronic couplings are relatively small even though the intermolecular distances are very short. For C 5 -PyDI, the V h and V e values are similar, both of which are less than 20 m eV; for C 6 -PyDI, the V e value is 26.3 meV and V h is only 7.2 meV. In order to understand the weak electronic couplings in the one-dimensional -stacking motif, the shapes of the HOMOs and LUMOs in the P dimer of C 6 -PyDI are shown in Fig. 3. We can see that the HOMOs and LUMOs are mainly localized on the pyrene core, and the introduction of alkyl chains has little influence on the frontier molecular orbital charge distributions. In contrast to the perfect face-to-face dimer, there exists an obvious displacement between two neighboring molecules of the P dimer along the molecular axis direction, which leads to a cancelation effect between the bonding and antibonding overlaps and a decrease in the effective coupling projected area. As shown in Fig. 3(b), we can see that the major contribution to the V h value mainly comes from NÁ Á Á interactions, and Á Á Á interactions in the P dimer contribute little to the V h value due to the cancellation effect. In comparison, the distribution character of the LUMO leads to intermolecular OÁ Á Á interactions and partially effective Á Á Á interactions in the P dimer, as show in Fig. 3(d), which explains the fact that the V e value in the P dimer of C 6 -PyDI is about three times larger than the corresponding V h value. Through the above analysis, we can conclude that the weak electronic couplings in the C 6 -PyDI crystal are mainly due to the large Calculated variations in the bond lengths of isolated C 6 -PyDI upon oxidation (black symbols) and reduction (red symbols), the x axis represents the chemical bond (C-C, C-O and C-N) number that is marked using different numbers.

Table 2
Calculated electronic coupling terms V h (hole transfer) and V e (electron transfer) for the different hopping pathways in C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI crystals; r is the intermolecular center-to-center distance. deviation of the P dimer from perfect face-to-face packing, which results in the compensation of bonding and antibonding interactions.
In the crystal structures of t-C 5 -PyDI and t-C 6 -PyDI, the molecules also pack into a lamellar motif; however, every two neighboring molecules form a molecular pair, and the distance between these molecular pairs is much longer than the intermolecular distance within the molecular pair. For example, the distance between t-C 6 -PyDI pairs is 11.629 Å , and the intermolecular distance in the t-C 6 -PyDI pair is only 3.519 Å . For the sake of discussion, the dimer with a small intermolecular distance, i.e. 3.519 Å , is defined as a P1 dimer; and a dimer with a larger distance, i.e. 11.629 Å , is defined as a P2 dimer. Our calculations show that the intermolecular electronic couplings in the P1 dimer are much larger than those in the P2 dimer, which is consistent with the intermolecular distances.
As shown in Table 2, the V h and V e values in the P1 dimer of t-C 5 -PyDI are 91.1 and 152.3 meV, respectively, which are much stronger than those in the P2 dimer (1.5 and 0.2 meV); similarly, the V h and V e values in the P1 dimer of t-C 6 -PyDI are 90.5 and 145.7 meV, respectively; in comparison, both V h and V e values in the P2 dimer are less than 0.01 meV. In this case, the hole/electron transport mobility is primarily determined by the hopping rate between different P1 dimers. The weak electronic couplings in the P2 dimers of t-C 5 -PyDI and t-C 6 -PyDI explain their poor OFET properties.

Anisotropic mobility
The anisotropic hole-transfer and electron-transfer mobility values in the single crystals of C 6 -PyDI and C 5 -PyDI are shown in Fig. 4. It can be seen that their similar crystal structures result in the same angle dependence of mobility; both C 6 -PyDI and C 5 -PyDI show remarkable anisotropic behavior and the highest mobility values appear when the value of È is near 0/180 , that is, along the crystallographic b axis. This mobility distribution as a function of È is consistent with their one-dimensional molecular packing character. The ranges of mobility values of C 5 -PyDI and C 6 -PyDI estimated in the same layer are summarized in Table 3.
We can see that the ranges of the hole-and electron-mobility values in C 5 -PyDI and C 6 -PyDI crystals agree well with the experimental measurements by Wu et al. (2017), which verifies the rationality of our computational method and strategy. Comparison of the predicted hole and electronic mobility values for C 5 -PyDI and C 6 -PyDI indicates that the holes in C 5 -PyDI are intrinsically more mobile than the holes in C 6 -PyDI; while for the electron, the mobility values in the crystal of C 5 -PyDI are obviously lower than those in C 6 -PyDI. Combined with the experimental results reported recently, it can be seen that the holeand electron-transfer mobilities of C 5 -PyDI and C 6 -PyDI in the experiments of Wu et al.
(2017) might have reached their optimum values as n-type or ambipolar materials.
For t-C 5 -PyDI, we simulated the angular resolution anisotropic mobility for both electron and hole transport. As shown in Fig. 4(b), it can be seen that the hole-and electrontransfer mobility values in the ab plane show (a) Illustration of the orientation angle of the transistor channel relative to the crystallographic b axis for C 5 -PyDI and C 6 -PyDI and the orientation angle of the transistor channel relative to the crystallographic a axis for t-C 5 -PyDI. Calculated angle-resolved anisotropic hole (black line) and electron (red line) transport mobility as a function of orientation angle for (b) t-C 5 -PyDI, (c) C 6 -PyDI and (d) C 5 -PyDI. Table 3 Theoretical hole-diffusion mobilities ( h ) and electron-diffusion mobilities ( e ) of C 5 -PyDI, C 6 -PyDI and t-C 5 -PyDI at room temperature (T = 300 K), and some experimental charge-transfer mobility values.

Figure 3
HOMOs (0.02 a.u.) for the P dimer of C 6 -PyDI in (a) side view and in (b) top view, and LUMOs (0.02 a.u.) for the P dimer of C 6 -PyDI in (c) side view and in (d) top view. similar anisotropic behavior: the highest and lowest mobility values were present at È = 0/180 and È = 90/270 , respectively. The maximum hole-transfer mobility value is only 0.004 cm 2 s À1 V À1 , and the maximum electron-transfer mobility value is less than 0.0001 cm 2 s À1 V À1 , which is consistent with the reported highest electron transport mobility (8.72Â10 À5 cm 2 s À1 V À1 ). Although the introduction of the tbutyl group has very slight effects on the reorganization energy associated with charge-carrier transfer, the steric effect induced by its large volume dramatically changes the relative position between adjacent molecules, which has a large interference on the size of V h and V e values. As a result, the OFETs fabricated from t-C 5 -PyDI exhibited lower conductive performance.

Ionization potential and electronic affinity
Aside from mobility, charge injection efficiency is also an important factor that affects the performance of an OFET device, especially for the ambipolar and n-channel OFETs. For OFETs, it requires that the electrode materials have work functions suited for the injection of holes/electrons into the HOMO/LUMO of semiconductor molecules. In molecular orbital theory approaches, the HOMO energy is related to the IP by Koopmanns' theorem and the LUMO energy is used to estimate the electron affinity (ÀE HOMO = IP and ÀE LUMO = EA); however, the ÀE HOMO /ÀE LUMO values are usually inconsistent with IP/EA values in the practical DFT calculations, partly due to the unknown 'exact' exchange-correlation functional. Previous calculations by Zhan et al. showed that the directly calculated vertical IPs are, on the whole, in good agreement with the corresponding experimental IPs (Zhan et al., 2003). Thus, we select IP and EA values of the studied compounds as the evaluation parameters to analyze the charge-injection barrier of the OFET.
For the metal electrode, the key for efficient injection of charge-carrier is that the IP values should be close to or smaller than the work function of the metal electrode; and the EA values should be close to or larger than the work function of the metal electrode. Considering the thermal and oxidative stability of electron-transport materials, the suitable EA values need to be at least 2.8 eV (Newman et al., 2004). The calculated IPs and EAs of C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI are shown in Table 4. We can see that the AIP values of C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI are larger than 7.0 eV. In this case, it needs a strong external electric field to overcome the high hole-injection barriers. It is noteworthy that the introduction of the t-butyl group in C 5 -PyDI and C 6 -PyDI could decrease their AIP values by 0.25 and 0.27 eV, respectively, which is favorable for hole injection. For the predicted EA values, it is observed that (i) all AEA values of C 5 -PyDI, C 6 -PyDI, t-C 5 -PyDI and t-C 6 -PyDI are larger than 1.9 eV, suggesting PyDI can be selected as electron-deficient frameworks for n-type or ambipolar semiconductors; and (ii) addition of the t-butyl group causes a slight decrease in AEA values of C 5 -PyDI and C 6 -PyDI, which is unfavorable to the electron-injection process.

Electronic spectra
Based on the optimized geometries, we simulated the UVvis absorption spectra of t-C 5 -PyDI in the gas phase and in chloroform solvent, as shown in Fig. 5(a). It can be seen that the calculated absorption peak positions in chloroform solvent are similar to the case in the gas phase. Their maximum absorption bands are observed in the near-UV region, which is consistent with the experimental excitation wavelength (270 nm) obtained by Wu et al. (2017). To gain a deeper understanding of the spectra character, the excited states and the frontier molecular orbitals (FMOs) were analyzed in detail, which can provide information about the nature of the excited-state conformations. For t-C 5 -PyDI in chloroform, as shown in Table 5, there exists an S 0 !S 1 transition at 392 nm, an S 0 !S 3 transition at 350 nm and an intense S 0 !S 9 transition at 271 nm. The S 1 state of t-C 5 -PyDI is mainly formed by the transitions from the HOMO to the LUMO; the S 3 and S 9 states are predominantly formed by the transitions from HOMO-1 to the LUMO, and by the transitions from the HOMO to LUMO+2, respectively. Moreover, the fluorescence spectrum of t-C 5 -PyDI is also simulated and the results are shown in Fig. 5(b). It is found that the predicted fluorescence maximum of t-C 5 -PyDI, corresponding to the S 1 !S 0 transition process, is located at 485 nm, which is quite close to the experimental value of 480 nm (Wu et al., 2017).  Table 4 VIPs, AIPs, VEAs and AEAs calculated at the B3LYP/6-311** level (eV).  Figure 5 (a) Calculated absorption spectra of t-C 5 -PyDI in gas and in chloroform.

Table 5
Electronic excitation energies (), oscillator strengths ( f ), corresponding compositions and the configuration interactions (CI) for t-C 5 -PyDI in gas and in cyclohexane.

Conclusions
In this manuscript, we simulated anisotropic charge-transfer mobilities of C 5 -PyDI, C 6 -PyDI and t-C 5 -PyDI, and theoretically predicted the range of their mobility values, which provided reference for the performance optimization of an OFET based on these materials. We systematically analyzed the influences of the alkyl chain on the reorganization energies, crystal packing, electronic couplings and charge injection barrier of PyDI. It is found that the introduction of alkyl chain groups into PyDI has little effect on the reorganization energy, but increases the repulsive interactions of the backbone and thus affects the molecular packing in the crystal. As a result, the electrons in the t-C 5 -PyDI crystal are intrinsically less mobile than the electrons in the C 5 -PyDI crystal. Moreover, the alkyl chain substitution could decrease the electronic affinities and ionization potentials, which decreases the holeinjection barrier and improves the injection efficiency of the holes. We also simulated electronic spectra of t-C 5 -PyDI which reproduced the experimental absorption and fluorescence spectra well.