Comprehending radicals, diradicals and their bondings in aggregates of imide-fused polycyclic aromatic hydrocarbons

Quantum effects such as ferromagnetism were regarded as rare in organic materials. When reduced to radical states, imide-fused polycyclic aromatic hydrocarbons (IPAHs) have shown room-temperature ferromagnetism in our recent work, to be a potential candidate as ferromagnetic semiconductor. Here, we use variational Davydov ansatz parametrized by density functional theory to investigate the structural and optical properties of IPAHs and their radicals at both molecule and aggregate levels. Our calculation reveals that hydrogen mainly gives rise to radicals and proves the formation of a mid-gap polaronic state, which is further evidenced by UV-vis absorption spectra simulations, in good agreement with experiments. The significant change of dispersion between the π–π stacking structure and planar structure implies the formation of radical–radical bonding (pancake bonding), which is revealed by simulations of NIR absorption signals and serves as the physical basis of long-range ferromagnetic orders. Absorption spectra of perylene diimide (PDI), terrylene diimide (TDI) and their radicals are also predicted.

where  (1) ,  (2) and  (3) represent the vertical energies of HOMO, LUMO and LUMO+1, respectively.  (6) In this article, we find the maximum shift of all atoms moving in three directions in all vibration modes at Cartesian coordinates, and limit the maximum shift in this direction to no more than -0.02, -0.01, 0.01 and 0.02 Å respectively. We set that all atoms in all modes change in equal proportion and get structures moved by four distances. Then, we use the quadratic polynomial to fit the energy differences before and after the moving atoms. 0 is the fitted constant term, which is always 0, as it is the energy difference at equilibrium geometry. are the first-order and secondorder coupling coefficients at the n-th energy level of the q-th vibration mode respectively. Where, , qs r represents the shift of the s-th atom of the q-th vibrational mode. s m is the mass of the s-th atom.

Expression of Davydov Ansatz and Observable Quantities
In this work, we adopt the variational Davydov ansatz to study the excitation dynamics and obtain UV-Vis absorption spectra of IPAHs, which has been shown its advantages in simulating the quantum dynamics in numbers of chemical and physical domain. [25][26][27][28][29][30] Considering three electronic orbitals, the vibrations of the nuclei, and their couplings up to secondorder truncation, the multi-D2 method is taken into account in our calculations, which can be written as  (10) In which  is the variational parameters,  is the time derivative of  . The Lagrangian L is given as ( ) ( )  = − −  22 (11) The absorption spectrum can be obtained by

( )
Wt is called autocorrelation function which is the overlap between initial and final states. Then we can get spectral lines with the Fourier transformation of the autocorrelation function. In practical experiment, lines of the spectra are broaden caused by the resolution of the spectrometer. 30 In order to reproduce this effect, here we multiply

( )
Wt with a phenomenological damping curve for Lorentzian broadening, where  is a model-dependent damping time.

Parameters for Davydov Ansatz Calculation
By DFT calculation, we obtain frontier orbital energies below:   Considering that multiple IPAHs systems need to be calculate, and the vibration modes of NDI with the minimum atomic numbers reach 78, it is necessary to neglect the weak coupling mode to simplify the calculation. We only take the mode with for calculation, and the corresponding frequency and coupling strength are as follows:        Table S8. Frequency (  q ), the first-order electronic-vibrational coupling strength of HOMO (   Table S9. Frequency (  q ), the first-order electronic-vibrational coupling strength of HOMO (   Table S10. Frequency (  q ), the first-order electronic-vibrational coupling strength of HOMO (        In our calculation, the couplings, ( ) 12 uu between orbitals are expressed by the ratio of corresponding transition dipole moments, and the time-dependent normalization coefficients at t = 0, ( )  Then, the orbital energy couplings, u1 (u2) and the time-dependent normalization coefficients at t = 0, Ai (0), Bi (0) and Ci (0) can be expressed below: For NDIHs anion vertical dimer, it is found that at the NIR region, the main transition contribution is S0-S5/S6, which mainly comes from the orbital transition of HOMO to LUMO+5/LUMO+6. In the meantime, the energies of LUMO+5 and LUMO+6 are basically equal. Therefore, only two orbitals, HOMO and LUMO+5, are considered in the Davydov ansatz calculation. In this case, many second-order coupling strengths are 0, and a few large second-order coupling strengths make the calculation time much longer. For improving the calculation efficiency, only the first-order electronic-vibrational coupling strength which satisfies is considered here.

Experimental Spectrum Characterization
The UV-vis absorption spectrum characterization was performed on NDI ionized derivatives with D2O, because OH overtone bands when using H2O would interfere the absorption spectrum characterization and D2O could eliminate the absorption peak of H2O. 31 The NDI ionized derivatives solution in hydrazine hydrate was heated in glovebox in 80°C to remove solvent, afterward the NDI ionized derivatives powder was dissolved in degassed D2O in concentration of 10 -4 mol L -1 in quartz cuvette, which was sealed to avoid contact with O2 and characterized in ambient condition. Although the diradical characteristics of NDIH-diradical cannot be described at the DFT, HSE06/6-311G** level (    For obtaining NDIHs vertical dimer, we set the vertical stacking distance to 3.45 Å according to the common π-π stacking distance, 34 and then set the sliding distances along the Y direction as 0.0 to 7.0 Å to move the monomer of the second layer, respectively. After geometric optimization, these structures can be divided by geometric center distance and energy into three groups (Fig. S5). The second group which involves the structures with an initial slip of 4.0 and 5.0 Å in the Y direction, has the lowest energy and the shortest geometric center distance. When the slip distance exceeds a limit of 6.0 Å, the attraction between the two monomers decreases, and the dimer dissociates, resulting in a sharp increase in the geometric center distance and energy. The geometric center distance of the first group structures is between the second group and the third group, but the energy is the highest. We can find that the structures of the first group and the second group are relatively similar. The difference is that the bending of the second group structures is more intense. It is considered that the strong radical-radical bonding reduces the energy of the second group. Therefore, radical-radical bonding is discussed in detail in our article.

Supporting Figures and Tables
In the following context, we only consider the vertical dimer with the lowest energy, taking the NDIHs dimer with initial slip of 4.0 Å as an example, and call it NDIHs vertical dimer.
For further explore NDIHs dimers, we compare the ESP distributions and orbital images of NDIHs horizontal dimer and NDIHs vertical dimer. It is observed that ESP distribution range of NDIHs vertical dimer is narrower, but the percentage of positive ESP in the total area is larger than that in the horizontal case (Fig. S6). For NDIHs horizontal dimer, the energy of two adjacent energy levels is similar. For example, the energy of HOMO-1 is -4.806 eV, the energy of HOMO is -4.780 eV, the energy of LUMO is -2.160 eV, and the energy of LUMO+1 is -2.133 eV. This makes the electron density distribution of adjacent orbitals very similar (Fig. S7a). The difference is that the wavefunction phase of one of two monomers in adjacent orbitals is opposite. For NDIHs vertical dimer, the isosurface of charge density is also similar for different phases between two adjacent orbitals, such as HOMO-1 and HOMO. The overlaps of electron clouds in the orbitals of LUMO, LUMO+1, LUMO+2 and LUMO+3 further prove the formation of pancake bonding.   Table S20. Interaction energies of NDIHs horizontal dimer and NDIHs vertical dimer, obtained by B3LYP-D3BJ/ 6-311+G** and sSAPT0/jun-cc-pVDZ, respectively. B3LYP-D3BJ/ 6-311+G** (kcal/mol) sSAPT0/jun-cc-pVDZ (kcal/mol)  For NDIHs tetramer, we also set the initial vertical stacking distance to 3.45 Å and the sliding distances between two layers along the Y direction as 0 to 7.0 Å. The energy of 0.0 Å, which represents the initial vertical stacking, is the lowes t (Fig. S8).
The molecular orbital from HOMO-3 to LUMO+3 can be simply divided into two groups (Fig. S9). The first group is HOMO-3 to HOMO. The orbitals in this group are similar to the superposition of monomer NDIHs LUMO orbits. The electronic clouds of HOMO-3 and HOMO-2 are concentrated on the two NDIHs at the bottom, and those of HOMO-1 and HOMO are concentrated at the top. The electron density of HOMO-3 and HOMO-1 overlaps in opposite phases, while that of HOMO-2 and HOMO overlap in the same phase. The second group is from LUMO to LUMO+3, and their distribution is similar to that of the first group with different electron density. The phase of LUMO and LUMO+2 is same, while that of LUMO+1 and LUMO+3 is opposite. ESP distribution of NDIHs tetramer shows the negative region is larger (Fig. S10).   Table S21. Interaction energy of NDIHs tetramer obtained by B3LYP-D3BJ/ 6-311+G** and sSAPT0/jun-cc-pVDZ, respectively. B3LYP-D3BJ/ 6-311+G** (kcal/mol) sSAPT0/jun-cc-pVDZ (kcal/mol) [a] The pair division can be seen in Fig. S8.
For NDI, PDI and TDI, ESP minimum values decrease with the increase of atomic numbers. Starting from the first strip whose central coordinate is positive, the percentage of positive ESP of NDI is 68.15%, that of PDI is 72.42%, and that of TDI is 75.27% (Fig. 1e; Fig. S11).          Table S30. Optimized Cartesian coordinates (Å) of NDIHs tetramer of initial slip 0.0 Å at HSE06/6-311G** level in the gas phase, charge = 0, spin multiplicity = 1.