Impact of Charge-Resonance Excitations on CT-Mediated J-Type Aggregation in Singlet and Triplet Exciton States of Perylene Di-Imide Aggregates: A TDDFT Investigation

The modulation of intermolecular interactions upon aggregation induces changes in excited state properties of organic molecules that can be detrimental for some optoelectronic applications but can be exploited for others. The time-dependent density functional theory (TDDFT) is a cost-effective approach to determining the exciton states of molecular aggregates, and it has been shown to provide reliable results when coupled with the appropriate choice of the functional. Here we apply a general procedure to analyze the aggregates’ exciton states derived from TDDFT calculations in terms of diabatic states chosen to coincide with local (LE) and charge-transfer (CT) excitations within a restricted orbital space. We apply the approach to study energy profiles, interstate couplings, and the charge-transfer character of singlet and triplet exciton states of perylene di-imide aggregates (PDI). We focus on the intermolecular displacement along the longitudinal translation coordinate, which mimics different amounts of slip-stacking observed in PDI crystals. The analysis, in terms of symmetry-adapted Frenkel excitations (FE) and charge-resonance (CR) states and their interactions, discloses how the interchange of the H/J character for small longitudinal shifts, previously reported for singlet exciton states, also occurs for triplet excitons.

The role of CT states on photo-induced processes is documented by a large number of investigations. For instance, CT states might facilitate intersystem crossing through spin-orbit coupling [33] and are often involved in the first step of singlet fission [34][35][36][37][38][39][40] through the process of symmetry-breaking charge separation (SBCS). CT states also have a crucial role in exciton dissociation and charge separation in heterojunctions between electron-donating and electron-accepting materials or in homojunctions between crystalline domains with different orientations [41][42][43].
While the modelling and analysis of singlet spin exciton states has received considerable attention [20,22,27,[44][45][46][47][48][49][50], triplet excitons have received comparably less attention and only few investigations have been reported [22,51,52]. A relevant role of triplet states in photoinduced processes is, however, proven and supported by experiments. Enhanced triplet-state generation, following photo-induced charge transfer, was reported by various groups in electron donor-acceptor polymer-blend films in organic photovoltaic devices where nonemissive triplet excitons are responsible for nonradiative-charge recombination [53]. Thus, triplet formation by charge recombination can be a detrimental, unwanted side effect but can also be exploited as a useful triplet-generation method, for instance in photocatalysis and photodynamic therapy [54,55]. Examples are, among others, triplet-state generation in molecular dyads [56] and molecular dimers in which dimerization-induced triplet population has been invoked [57].
A deep insight into the relationship between structural arrangement and triplet excitonstate character is, therefore, essential to identify molecular systems and condensed phase organizations that are most suitable to generate and exploit triplet exciton states.
Perylene-bis(dicarboximide) (PDI) and its derivatives have attracted great interest as chromophores for energy and charge-transport studies, thanks to their propensity to self-organize into ordered assemblies, both in solution and in the solid state via π−πstacking interactions [11,12,58]. Numerous computational investigations on PDI aggregates have been focused on the prediction of exciton states with different quantum-chemical (QC) approaches, including configuration interaction truncated to single excitations (CIS), time-dependent density functional theory (TDDFT) [22,27,46,47,59,60], and highly accurate levels of theory [61][62][63].
TDDFT is a cost-effective approach to evaluating exciton states of molecular aggregates and has been shown to provide reliable results when coupled to a suitable choice of the exchange-correlation functional. To analyze the character of excitonic states predicted by QC calculations, a diabatization procedure can be used to determine the superposition of LE and CT diabatic states in each adiabatic exciton state [21,22,31,32,59].
For aggregates characterized by a symmetric arrangement of chromophores, as those investigated here (Figure 1), when the molecules approach each other, intermolecular interactions mix the LE states to form a symmetry-adapted (SA) superposition of (neutral) LE states, that is, Frenkel excitons (FE). Similarly, the ionic states form delocalized chargeresonance (CR) states [18,64,65] of appropriate symmetry. The symmetry point group of a PDI aggregate, where intermolecular displacements along the longitudinal translation coordinate, z, are considered, (Figure 1) is C 2h . As a result, the most relevant ππ * exciton states, along with the FE and CR diabatic states, all belong to A g and B u symmetry representations. Symmetry can be exploited to easily identify H-and J-aggregation types, since only B u exciton states are accessible as dipole-allowed transitions. Here we adopt a similar strategy to characterize singlet and triplet exciton states of PDI aggregates. To this end, we analyze exciton states in terms of LE and CT contributions, along the lines described in ref. [60], with three major objectives. First, we validate the protocol of ref. [60] by comparing the computed exciton states character (for both singlet and triplet exciton states) with the results obtained from the character-analysis tool, TheoDORE [66]. Second, we show that for triplet exciton states the modulation of CR/FE interactions along the longitudinal translation coordinate determines a switch in the symmetry (A g /B u ) of the lowest triplet exciton state, which corresponds to the unconventional H-/J-character alternation previously documented for singlet excitons [8,9,60]. Finally, we demonstrate that the selection of the functional is less critical for triplet excitons compared to singlet excitons and show that the CAM-B3LYP and ωB97X-D results are very close to each other.

Computational Methods
The protocol employed to analyze the character of exciton states transforms TDDFT amplitudes from the basis of single excitations between the aggregate's orbitals (the delocalized excitation (DE) basis) to the basis of single excitations between molecular site orbitals, the latter defining the diabatic states. It has been previously described [60] and is briefly summarized here. To analyze the exciton character, we express each relevant exciton state in terms of LEs and CTs. To this end, we select the orbital subspace corresponding to relevant ππ * exciton states. For PDI aggregates this includes the HOMO and LUMO of each monomer [21,31,59,60,67] and represents the minimal orbital space (MIOS). The next step involves the determination of the linear combination coefficients, C AGGR_MOB i,j , forming the C AGGR_MOB matrix and describing each aggregate's orbital in the monomer orbital basis (MOB) as [60,68,69]: where the C MON_AOB matrix is a block diagonal matrix containing the molecular orbitals' (MOs) coefficients in the atomic orbital basis (AOB) of each monomer, with off-diagonal blocks set to zero, and S MON_AOB is the overlap matrix of the monomers in the AOB. In general, monomer orbitals belonging to two different molecules are nonorthogonal to each other. Hence, the aggregate's orbitals, C L AGGR_MOB , expressed in terms of orthogonalized monomer orbitals are obtained as: where superscript L indicates Löwdin's orthogonalization [70], and the overlap matrix S AGGR_MOB = C t MON_AOB ·S AGGR_AOB ·C MON_AOB is obtained from the coefficients of the monomer's orbitals, C MON_AOB , and the overlap of the atomic orbitals in the aggregate configuration, S AGGR_AOB .
From the results of the TDDFT calculations on the aggregate, the subset of the n 2 exciton states originated from the MIOS are then selected out of the full set of computed eigenstates. TDDFT amplitudes are expressed based on the DEs, namely excitations between the aggregate's orbitals, and form the columns of the B adia DE matrix. Thus, the following step is required to expand each DE in terms of diabatic LE and CT states (excitations between monomer orbitals). With the aggregate's orbitals expressed in terms of monomer orbitals via the C L AGGR_MOB matrix, each excitation from an occupied i to an empty j aggregate's orbital DE( i → j ) can be expressed as a linear combination of diabatic (LE and CT) excitations (k → l) from an occupied k to an empty l monomer orbital, with the expansion coefficients forming the columns of the unitary matrix U DE →dia given by Exciton states are then readily expressed in the diabatic basis as B adia dia = U DE→dia ·B adia DE , and their character is obtained by summing up the contributions from the CT and LE states.
The corresponding n 2 eigenvalues (excitation energies of the selected excitons) form the diagonal H adia matrix, from which the Hamiltonian in the diabatic LE/CT basis, H dia , can be obtained as [21,30,32,71] Finally, the H dia is rotated in the SA diabatic basis formed by FE and CR excitations to obtain two matrices, H sym dia (one for B u states and one for A g states), whose off-diagonal elements are the interactions between the CR and FE states that ultimately govern the modulation of adiabatic exciton-state energies along the longitudinal translation coordinate. These interaction energies have been shown to correspond to ± combinations of electron-(D e ) and hole-(D h ) transfer integrals [22,31,59,60].
The PDI-monomer structure was the same one used in previous investigations on PDI aggregates, optimized at the BLYP-D/TZV(P) level of theory [72]. The distance between the planes of different monomers was set to 3.4 Å, which is a distance used in previous investigations on dimers of PDI [59,60]. Exciton states were computed for the eclipsed aggregates and for displacements of 0.5 Å up to 8.0 Å along the longitudinal translation coordinate (z) (Figure 1). For singlet exciton states, trimers and tetramers were also considered beside dimers ( Figure S1).
Excitation energies were determined with TDDFT calculations in the Tamm-Dancoff approximation (TDA) [73], using the CAM-B3LYP [74] and ωB97X-D [75] functionals and the 6-31G* basis set. The ωB97X-D functional was previously shown to provide a reliable description of CT character in singlet excitons of PDI dimers [27]. All QC calculations were carried out with the Gaussian16 suite of programs [76]. The CT character of exciton states was determined as described above and using TheoDORE 2.4 [66].

Modulation of Singlet Exciton-State Energies and Characters for PDI Aggregates
Previous investigations of excitation energy profiles and CT character in PDI dimers for interchromophore longitudinal shifts have shown that computed exciton states are generally strong mixtures of the two types of diabatic states. The weight of LE and CT characters, however, strongly depends on the energy difference between FE and CR diabatic states, and on their couplings [27,59,60]. This is the reason why the CT character of the lowest singlet exciton state is critically dependent on the chosen functional. In this regard, it has been shown that the adiabatic energy profiles and CT characters computed at the TD-ωB97X-D/6-31G* level match those obtained at higher level of theory [27,59,60]. Here we performed the character analysis based on TDA-computed exciton states rather than the full TDDFT linear response. While the excitation energy profiles of the four exciton states are quite similar (Figure S2), the FE and CT energy profiles of B u symmetry obtained from TDA-ωB97X-D/6-31G* calculations display a crossing (Figure 2) not seen for TD-ωB97X-D/6-31G* calculations [60], leading ultimately to a CT character of the lowest B u exciton state (at the eclipsed geometry and for small displacements, see Figure 3) greater than 50%, since the weight of the CT character can be traced back to the energy location of FE and CR states [60] (Figure S3). Such large CT character contrasts with experimental results and was previously reported for other long-range corrected functionals and from CIS calculations [27,59]. Obviously, slightly larger intermonomer separations will reduce the contrasting results between TD and TDA-ωB97X-D/6-31G* calculations, which, however, illustrate the crucial role determined by the energy difference between the CR and FE states.  Finally, we note that the CT character determined following the procedure outlined in Section 2 and the results of TheoDORE [66] are in excellent agreement for all the four low-lying states investigated (Figure 3). A similar correspondence is found for the lowest energy exciton states of larger aggregates ( Figure S4).

Modulation of Triplet Exciton-State Energies and Characters for PDI Dimers
Triplet exciton-state analysis was carried out only on TDA calculations since TD triplet wavefunctions contain relevant contributions from de-excitations. The oscillating trend in adiabatic energy profiles (Figure 4a) are easily rationalized by the oscillating trend of the interactions between the CR and FE states ( Figure 4b). As shown previously [22,31,59,60], these interactions are given by the ± combinations of the D e and D h transfer integrals. While couplings between triplet LE states are almost negligible [22], spin triplet CR/FE interactions are of the same magnitude as those computed for singlet exciton states ( Figure S5). In analogy with singlet excitons, these interactions determine not only the CT character of the adiabatic triplet exciton states (Figure 4c), but also the symmetry of the lowest energy state (Figure 4a), with an interchange of A g and B u along the longitudinal translation coordinate. Such an alternation is well established for singlet exciton states and leads to the appearance of CT-mediated J-type spectroscopic features for small longitudinal displacements of one of the two PDI molecules [8,9,58,60] in a range of displacements where Kasha's theory only predicts H-type aggregation. Therefore, CR states create an effective short-range exciton coupling that can induce unconventional J-type spectroscopic features [8,9,46,[58][59][60]69,77,78]. Notably, here we show that the same mechanism is active in the triplet exciton manifold, albeit with less effective CR/FE mixing. Indeed, as previously reported [22], one of the most relevant differences between singlet and triplet exciton states is the larger energy difference between FE and CR states, which triggers a modest CR/FE mixing. This can be appreciated in Figure 4c, showing that the largest CT contribution for the lowest triplet exciton state does not exceed 20%. At the same time, Figure 5a,b compare the results of CR/FE interactions for singlet and triplet A g exciton states of the PDI dimer and enlighten the role of the increased energy difference between FE and CR states in the latter. Nevertheless, the CR/FE mixing is sufficient to modulate the triplet exciton energies, ultimately leading to an alternation between B u and A g states that parallels what is well documented for singlet exciton states.  (Figure 4b) are strong, the exciton-state energies deviate from FE and CR energies. The interaction is less effective for triplet states due to the larger energy difference between FE and CR states.
More specifically, the oscillation of the CR/FE interaction (see the D e ± D h profiles shown in Figure 4b), determines an interchange of the lowest exciton states of A g and B u symmetry in the dimer for longitudinal shifts in the range of ca. 2-3 Å, where the interaction between A g symmetry FE and CR states is almost negligible. At the same time, the interaction between the FE and CR B u states is the maximum and pushes the adiabatic 1B u state energetically below 1A g (see Figure 4a), thereby switching the aggregate character from H-to J-type.
These results suggest that radiative and nonradiative decays from the lowest triplet exciton states of molecular aggregates, intimately related to the nature and symmetry of the lowest energy state, may be modulated by the intermolecular organization, a concept that could be exploited for systems displaying dimerization-induced triplet state populations [57].
We finally compare the CT character analysis obtained from our protocol (solid lines in Figure 4c) with the results of TheoDORE (dashed lines in Figure 4c). The agreement is very good for the two lowest triplet exciton states, while marked differences appear for the remaining two states. The reason is that for several intermolecular configurations the selected ππ * triplet exciton states (whose wavefunction is dominated by excitations within the selected ππ * orbital space) are spread over more than one TDDFT-computed exciton state. As a consequence, TheoDORE analysis picks up only the CT contribution corresponding to the fraction of the selected ππ * exciton state. In contrast, with the procedure described in Section 2, the selected ππ * exciton states are projected out of the entire space of computed exciton states and renormalized. Interestingly, triplet exciton states (belonging to the selected orbital space) appear to mix more strongly than singlet states with other exciton states of the dimer, as suggested by wavefunction analysis and by the fragmentation of the CT contributions in TheoDORE analysis. Such fragmentation contrasts with the almost perfect match discussed for singlet excitons of the PDI dimer ( Figure 2) and suggests that diabatization procedures are likely to be more problematic when analyzing triplet exciton states.
To visualize the nature of exciton states for two representative dimer configurations of the PDI dimer, we have carried out a fragment-based analysis via electron-hole correlation plots using TheoDORE ( Figure 6). The two selected fragments correspond to the two molecules forming the dimer. Exciton states are identified by the nonvanishing elements of the 2 × 2 matrix (the Ω-matrix [66]) represented by different levels of grey. Locally excited contributions appear on the main diagonal (going from lower left to upper right), while CT contributions appear off diagonally. In agreement with the character analysis shown in Figure 4c, Figure 6 shows that for the eclipsed configuration, the character of the 1B u state is purely LE while the character of the lowest energy exciton state of 1A g is partially CT, as indicated by the light grey squares in the electron-hole correlation plot. The situation is reversed at the slip-stacked configuration (2.5 Å longitudinal translation), with the 1B u state displaying some CT contribution.

Influence of Functional on Triplet Exciton Character of the PDI Dimer
While the interplay of CR and FE states is crucial for TDDFT functionals to reproduce higher level results for the calculation of singlet excitons [27,60], the larger energy difference between FE and CR triplet states makes the choice of the functional less critical. This is illustrated by the small differences obtained between TDA-ωB97X-D/6-31G* and TDA-CAM-B3LYP/6-31G* calculations.
The major difference between the two functionals is the lower excitation energies computed at the TDA-CAM-B3LYP level. The difference between the TDA-CAM-B3LYP and TDA-ωB97X-D energies is larger for the two higher energy exciton states, whose CT character is dominant (Figure 7). The lower energy of CT states at the TD-CAM-B3LYP level is not unexpected and was identified as the major source of discrepancy for singlet exciton states computed with the CAM-B3LYP functional compared to ωB97X-D [27]. As a result of the lower energy of CT states with the CAM-B3LYP functional, the CT character of the lowest energy triplet exciton states is slightly larger at the TDA-CAM-B3LYP level, compared to TDA-ωB97X-D, as can be seen in Figure 8. Interestingly, the CR/FE interaction is virtually independent of the chosen functional, with almost identical TD-CAM-B3LYP and TDA-CAM-B3LYP values (Figure 9). While it can be expected that other long-range corrected hybrid variants would provide results similar to those reported here, the use of nonhybrid or hybrid functionals with small contributions of exact (HF) exchange is not recommended because of their tendency to underestimate the energetic position of CT states [46].

Conclusions
We have analyzed the modulation of singlet and triplet exciton states of PDI aggregates computed at the TDDFT level, focusing on the intermolecular displacement along the longitudinal translation coordinate, which mimics different amounts of slip stacking observed in PDI crystals.
The CT character of singlet and triplet exciton states has been determined with a diabatization procedure, and it has been shown that the results agree with other characteranalysis tools, such as the TheoDORE software.
The study has shown that triplet exciton-state energies can be rationalized in terms of the interactions between CR and FE diabatic states, which are found to be of the same strength as those computed for singlet exciton states. Such CR/FE interactions ultimately lead, not only to a mixed LE/CT character of the low-lying triplet exciton states, but also to a CT-mediated J-aggregation mechanism for small longitudinal displacements, similar to what has been previously documented for singlet exciton states.
Finally, we have compared the results of two long-range corrected hybrid functionals and shown that the magnitude of the CR/FE interactions are almost independent on the selected functional. In contrast with singlet exciton states, the larger energy difference between FE and CR triplet states makes the choice of the functional (among long-range corrected hybrid variants) less critical to define the CT contributions in low-lying triplet exciton states.
Overall, these results demonstrate that the same interactions responsible for the CTmediated J-aggregation in singlet exciton states of PDI aggregates are operative in triplet exciton states and show that the role of CR states must be carefully considered as a factor influencing the processes following triplet exciton generation. These findings might be exploited to design chromophores and interchromophore organizations to optimize dimerinduced triplet state population.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/computation10020018/s1, Figure S1:The PDI trimer and tetramer considered in this work; Figure S2:Comparison between adiabatic energy profiles of the four singlet exciton states of PDI dimer along the longitudinal shift calculated at TDA-ωB97X-D/6-31G* (dashed) and TD-ωB97X-D/6-31G* (solid) levels; Figure S3:Comparison between FE and CR singlet states of Bu symmetry calculated at TDA-ωB97X-D (dashed) and TD-ωB97X-D (solid) levels. 6-31G* basis set is used; Figure S4:CT character of adiabatic singlet exciton states of PDI trimer (a,b) and tetramer (c,d) along the longitudinal translation coordinate determined with two different approaches; Figure  S5:Comparison between the modulation of the D e ± D h terms of singlet (dashed) and triplet (solid) exciton states of the dimer along the longitudinal translation coordinate.