Exciton transport in molecular organic semiconductors boosted by transient quantum delocalization

Designing molecular materials with very large exciton diffusion lengths would remove some of the intrinsic limitations of present-day organic optoelectronic devices. Yet, the nature of excitons in these materials is still not sufficiently well understood. Here we present Frenkel exciton surface hopping, an efficient method to propagate excitons through truly nano-scale materials by solving the time-dependent Schrödinger equation coupled to nuclear motion. We find a clear correlation between diffusion constant and quantum delocalization of the exciton. In materials featuring some of the highest diffusion lengths to date, e.g. the non-fullerene acceptor Y6, the exciton propagates via a transient delocalization mechanism, reminiscent to what was recently proposed for charge transport. Yet, the extent of delocalization is rather modest, even in Y6, and found to be limited by the relatively large exciton reorganization energy. On this basis we chart out a path for rationally improving exciton transport in organic optoelectronic materials.


Supplementary
a Basis set is fixed to 6-31g(d,p) as commonly used in the literature for similar systems 1-4 . b The geometry for a6T was kept coplanar as done in Ref. 2 since oligothiophenes tend to crystallize in the solid state with coplanar thiophene units 2 . c The calculated oscillator strengths are reported in a.u. d Experimental values are taken for solution and gas-phase molecules from Ref. 5,6 in the case of ANT, Ref. 7,8 for a6T, Ref. 9 for a similar PDI derivative, Ref. 10 for DCVSN5 and Ref. 11 for Y6.
The first important observation is the fact the two long-range corrected functionals give very similar values both for energies and oscillator strengths (f ) of the lowest-energy excitations of the different systems, pointing to a robust description of the transition density of these excitations by different functionals. As the following calculations are aimed at estimating correctly electronic couplings, as explained in the main text, a good evaluation of the transition densities is the most important issue. It is also worth noting that the second lowest (singlet) excited state as found from both functionals is about 0.5-1.0 eV above S1. This means that the Frenkel approximation describing the excited state of the system as a combination of locally excited S1 states is a good first approximation to study transport properties in these OSs.
We also plot the Natural Transition Orbitals (NTOs) 12 corresponding to the S1 excited state in Figure 1 of the main text for all the investigated molecules. The NTOs offer a useful way of visualizing which orbitals give the largest contribution to a given single-particle excitation (the NTOs are obtained here by using the NTOBuilder tool 13 ). For all systems we can see that the NTOs with the largest contributions are indeed essentially the same as the HOMO and LUMO orbitals of the single molecules.

Supplementary Note 2 Internal exciton reorganization energies
We have performed (internal) exciton reorganization energy calculations (Eq. 10 in the main text) with three well-established long-range corrected hybrid functionals: CAM-B3LYP, ωB97X-D and M06-2X. Results are shown in Supplementary Tables 2 and 3 and compared with literature data. It is clear that exciton reorganization energies are quite robust with choice of the functional and they change very little with increasing basis-set size. They also agree well with literature data.  14 , level of theory: SCS-CC/cc-pVTZ b Taken from Ref. 15 , reorganization energy from normal mode analysis, level of theory: ωB97X-D/6-31G(d,p) c Taken from Ref. 16 , level of theory: M06-2X/6-311g(d,p) d Taken from Ref 17 , level of theory: optimally-tuned range separated hybrid with dielectric screening.

Calculation of excitonic couplings
We analyse here various flavours of the excitonic coupling approaches as defined in the main text. The electronic coupling is one of the main ingredients of the Frenkel exciton Hamiltonian employed in this work. We evaluated excitonic couplings using the full Coulomb integral, V Coulomb , in Eq. 11 main text, the TrESP approach, V TrESP , in Eq. 12 main text and the point-dipole approximation (PDA), V PDA , in Eq. 8 below, for all the systems considered here. These values have been compared with the total excitonic coupling, V , obtained using the multi-state fragment excitation difference fragment charge difference approach (MS-FED-FCD) described below. The objective is to assess the extent of the short-range coupling in these closely packed OS solids and the different levels of approximation of the long-range Coulomb interaction. We note in passing that, to check the quality of the TrESP charges found from the ESP fitting procedure (as described in the main text), we have made sure that the sum of the TrESP charges is zero (as it would be the case if the full transition density were integrated over all space) and also that the dipole moments obtained directly from the charge densities, µ = dr I ρ T (r I )r I , are identical with those calculated from the atomic TrESP charges, µ = I q I r I .
For the calculation of MS-FED-FCD couplings (V ) and Coulomb contributions (V Coulomb ) as well as for fitting the TrESP charges, we have used CAM-B3LYP 18 functional for a6T, PTCDI-H and Y6 (as suggested in Ref. 4,19 ) and ωB97X-D 20 for ANT and DCVSN5 (for a more consistent comparison with Ref 1 ). The results are reported in Table 1 of the main text for three of the closest crystal pairs of the investigated OSs.
As reported in the main text, we observe that V Coulomb is generally very close to the total coupling V . The mean relative unsigned error (MRUE) is about 0.7% for the crystal structures (and about 7% along MD).
This small discrepancy can be attributed to the missing short-range part, V Short . The observation that lays the groundwork for efficient approximations of the long-range Coulomb interaction   using TrESP charges (see below).
Notably, we also found that the two long-range corrected functionals give very similar values for the Coulomb couplings of the analysed systems. This is because the transition densities are described similarly by both functionals.

Multi-state FED-FCD
In this section we present two useful diabatization schemes: the fragment charge difference (FCD) 21,22 and the fragment excitation difference (FED) 23 methods, and their combination in the general multistate FED-FCD (MS-FED-FCD) scheme, [24][25][26] used in this work to obtain total excitonic couplings (V in the main text) between maximally localized Frenkel exciton states, that can be used as a benchmark for more approximate methods as discussed in the main text.
We start by considering the case of a dimer composed of two molecules donor (k) and acceptor (l).
The FED scheme can recover the diabatic (localized) basis from delocalized excited states by using an additional operator ∆x, which measures the difference in excitation number between the donor (k) and acceptor (l) molecules. The elements of the ∆x matrix are given in terms of "excitation densities", defined as the sum of attachment (electron) and detachment (hole) densities 27 : where i and j are two adiabatic states and ρ ex ij is the excitation density, defined as, The quantity ∆x has its extrema when the excitation is entirely localized on either the donor (k) or acceptor (l). Without loss of generality, assuming that the adiabatic states i and j are the combination of two diabatic states, |k * l localized on k, and |kl * localized on l, the eigenvectors of the 2 × 2 ∆x matrix represent the transformation from the adiabatic to the diabatic basis (U), and the eigenvalues are either 1 or -1 for |k * l and |kl * 23,24,28 . The diagonal matrix of adiabatic energies (H ad ) can be transformed into the diabatic basis of |k * l and |kl * by U † H ad U = H. The electronic coupling between these states is found as the off-diagonal element of H.
The same strategy can be applied to computing the coupling to charge transfer (CT) states. The FCD method uses an additional operator ∆q to localize the electronic charge: 21 where ρ ij (r) is the transition density between states j and i, if j = i, and the state density if j = i.
Similarly to the FED case, the matrix ∆q can be diagonalized to obtain locally excited states (such as kl * and k * l) with eigenvalue 0, and CT states, such as k − l + or k + l − , with eigenvalues 1 or −1, respectively. This is similar to the FED when the adiabatic states are combination of only two diabatic states. However, in the systems analysed in this work a 2 state diabatization procedure is not sufficient to retrieve excitonic couplings between completely de-mixed and localized (Frenkel) exciton states, which in our case, form the state space for the Frenkel Hamiltonian in Eq. 2 main text.
To overcome this limitation we employed a multi-state FED-FCD procedure. This approach is a generalization of the FED and FCD methods just described 22,24,28 and can make use of multiple adiabatic eigenstates of the donor-acceptor supermolecular system to recover maximally localized and decoupled diabatic states. This strategy was previously successfully employed to recover excitonic couplings in light-harvesting and biological systems 24,25 as well as bridged donor-acceptor moieties 28 . For a more detailed description of this approach we refer to Ref. 22,24 . In brief this algorithm starts from the definitions of the additional operators ∆x and ∆q defined before. For CT states such as k − l + and k + l − ∆x should be zero, while ∆q = ±1 (albeit small numerical uncertainties 22 ). On the contrary, for localized Frenkel-Exciton (XT) states such as kl * and k * l, ∆x equals +1 and −1, respectively, and ∆q = 0. We can define the matrix which then has eigenvalues +1 and -1, respectively for subspaces of CT and XT states. The first step of the multi-FCD-FED scheme is to diagonalize D to separate CT and XT subspaces, similarly to multi-state FCD (previously developed in Ref. 22 ) : The same transformation U 1 is applied to ∆x and ∆q, to obtain ∆x and ∆q . Now, the states can be assigned as either CT or XT depending on the value of D . ∆x and ∆q can be each diagonalized within the XT and CT subspaces respectively, so that: where ∆x is diagonal only in the XT block, whereas ∆q is diagonal only in the CT block. U 2 is block diagonal and rotates the XT and CT subspaces separately. Within the XT subspace, the diagonal values of ∆x can be used to separate kl * (XT1) states from k * l (XT2) states. At the same time, ∆q can be used to divide CT states in k − l + (CT1) and k + l − (CT2). At this point, the four subspaces have been separated.
The diagonal electronic Hamiltonian E can be transformed to this basis as U 2 U 1 EU † 1 U † 2 = H . Finally, we require that the Hamiltonian of each subspace is diagonal, i.e. that the states of the same subspace are not coupled to each other. 22 The unitary matrix U 3 is block diagonal, and each block of U 3 diagonalizes each of the CT1, CT2, XT1, and XT2 blocks of the Hamiltonian. The final Hamiltonian H in Eq. 7 contains the diabatic energies on the diagonal, and the electronic couplings in the off-diagonal blocks.

Point dipole approximation
As it is customary in (dipole-based) Förster theory, at sufficiently large distances between sites (i.e. distance larger than the dimension of the interacting molecules), the transition densities of donor and acceptor can be approximated with the first non-zero term in a multipole expansion, i.e. the transition dipole moment. This approximation is referred to as point dipole approximation (PDA) and the related excitonic coupling, V PDA , can be written as, where r kl is the vector distance between k and l molecules and µ k and µ l their respective transition dipole moments. We found that the PDA approximation breaks down for all molecules except ANT and should not be used (see Supplementary Note 3).

Excitonic couplings The importance of including multiple states in MS-FED-FCD diabatization
It is important to note that in many cases, and also for the systems investigated in this work, a 2 state adiabatic basis, employed in the diabatization procedure to calculate the reference coupling V , is not

Excitonic coupling fluctuations
V TrESP obtained by using the ESP fitting procedure described in the Method section is in very good agreement with V Coulomb for all systems and crystal pairs investigated (see Table 1 in the main text). The notable difference is the fact that TrESP couplings are readily calculated without the need of repeating an electronic structure calculation for each different geometry. This constitutes a great advantage of using TrESP approach in combination with FE-SH, as it permits calculating many thousands of coupling elements (and related nuclear gradients) at each step along the MD, thus allowing the study of large nano-scale systems. Additionally, this approach allows an efficient calculation of both couplings and coupling fluctuations beyond nearest neighbour pairs and to account for long-range interactions within the (Frenkel) Hamiltonian of these systems. These interactions are of vital importance to get a reliable estimate of the diffusion constant as it was recently pointed out in Ref. 30 and Ref. 31  we computed the Coulomb couplings with updated transition densities at each MD step and compared them to the TrESP excitonic couplings with frozen transition charges along 100 snapshots taken from an MD simulation (excitonic couplings for crystal pairs P b , of ANT and a6T and P a , P a1 of PDI, DCVSN5, respectively). To obtain optimal TrESP charges that partially account for more "exotic" configurations during the parametrization, in this case, we carried out the ESP fitting procedure for 50 different structures along MD simulations for each system and then averaged the TrESP charges obtained over all structures (instead of using just a single molecular structure). A sign-tracking procedure was adopted to keep con-

Distance dependence of excitonic couplings
In order to explore the distance dependence and the long-range nature of the excitonic couplings, and to find how well the observations made before for the closest crystal pairs can be generalized to molecules further apart, we performed coupling calculations displacing the the two molecules forming the π-stacked pairs (P b , of ANT and a6T and P a , P a1 of PDI, DCVSN5, respectively) at various distances. The results are reported in Supplementary Figure 5 for some of the systems analysed in this work.
As expected, V , V Coulomb , V PDA and V TrESP , all perform very similarly at large distances, that is when the higher multiple terms of the transition densities can indeed be neglected. We note that the total and DCVSN5 have spatially extended transition densities and, when the distance between the sites falls below the actual spatial extent of the electronic transition density, the PDA breaks down and it should no longer be used (i.e. the actual shape of the transition density needs to be taken into account). This suggests caution in the application of Förster theory when describing exciton dynamics for molecular semiconductors composed of medium-sized and large organic molecules. On the opposite, the V TrESP (green line) gives very good coupling estimates within the full distance regime, proving again to be a useful method. We also note that, for all the systems investigated here, short range effects, V Short , remain negligible even below the shortest crystal pair distance (dashed black vertical line), and the V Coulomb is still a reasonably good approximation of the total V . Yet, short range effects might become important at even shorter distances when the overlap between donor acceptor wavefunctions becomes larger.

Transition dipoles
In Supplementary Figure 6 we show the transition dipoles obtained using TrESP charges. The + and − signs indicate an increase and decrease in bond length going from the ground to the excited system, respectively. The displacements are used to parametrize the force field for the molecules in their lowest energy singlet state. For clarity, displacements symmetric to the ones indicated are not shown. The scaling factor β for force field parametrization according to the TDDFT reorganization energy, as described in the text, is also reported for each system.

Convergence of exciton diffusion constants
Simulation details are given in Methods section of the main text.

Impact of isotropy and dimensionality on exciton transport
Using the transient localization theory framework 37,38 as well as direct non-adiabatic molecular dynamics simulations 35 in the context of charge transport in organic semiconductors, it was recently proven that charges diffuse more effectively in 2D (or 3D) systems featuring isotropic nearest neighbour electronic couplings, rather than in 1D-like materials with a sizeable electronic coupling along a particular transport direction. To asses the impact of dimensionality on the transport properties of excitons we performed exciton transport simulations on reduced 1D models. In particular, we focused our analysis on two extreme cases: ANT and Y6, displaying the lowest and the highest diffusion constants, respectively. For ANT we studied a 1D chain along b direction (see Figure 1 in the main text), while for Y6 a single pillar of molecules along a direction (out of the 4 pillars forming the rod-like nano-crystal adopted in this work) was activated to the transport of the exciton (see Methods). We found that in both cases the MSD and diffusion constant for reduced dimensionality systems are smaller compared to the respective larger models (see Supplementary Figure 10a Table 1 in the main text), thereby making this system already 1D-like and quite strongly anisotropic for exciton transport (see also Figure 4a in the main text). On the other hand, Y6 shows couplings that are of the same order of magnitude between different molecular pairs owing to its 3Dlike structure and favourable packing. Thus, in the latter case, artificially restricting exciton transport to take place along a single pillar of molecules makes the transport slower and less effective. The different coupling anisotropy characterizing the two systems is also reflected in the delocalization of the thermally accessible excitonic band states (see Supplementary Figure 10b

Delocalization of thermally accessible states as a function of interactions range
It has been recently reported that the long-range excitonic interactions can drastically modify the exciton diffusion constant in ordered nanofibers of P3HT depending on the length of the polymeric fibers. 30,31 Prodhan and co-workers 31  To assess the quality of FE-SH simulations we report in Supplementary Table 5  However a smaller reorganization energies of 302 meV was found in the former work, as a consequence of the approximate DFTB approach that relies on local functionals to describe the excited state geometry.
This reorganization energy was employed to model the exciton-phonon coupling leading to a weaker ex-citon relaxation compared to TDDFT used in FE-SH, where the exciton reorganization energy is 561 meV.
Nevertheless, despite these differences and the different underlying dynamics, FE-SH and BC-Ehrenfest methods yield similar diffusion constants (see Supplementary Table 5). The diffusion lengths from FE-SH and BC-Ehrenfest are both slightly underestimated compared to the experiments (see main text for a discussion on this point). a The diffusion length is estimated as described in the main text for all the different methodologies using an exciton lifetime τ = 10 ns, as obtained by experiments 40 . b FE-SH simulations are performed as described in Methods. c Boltzmann corrected Ehrenfest values are taken from Ref. 39 . The authors used reduced dimensionality 1D models along a and b crystallographic directions. d Perturbative rate theory result are taken from Ref. 1 . A rate expression very similar to the Marcus-Levich-Jortner rate for electron transfer was used by the authors by including an effective quantum mode coupled to the excitation energy transfer process. e Experimental estimated for the diffusion lengths are taken from Ref. 41,42 . f Excitonic couplings (averaged along MD trajectories) are given for P b in the crystal. g Reorganization energy from TDDFTB h Reorganization energy obtained using normal mode analysis based on Huang-Rhys factors evaluation and the same level of theory as for FE-SH.

Supplementary
In Supplementary Table 5 we also compare both numerical non-adiabatic propagation schemes with the perturbative rate theory used in Ref. 1 by Aragó et. al as well. The latter approach is reasonably well justified in this system by the high local exciton-phonon coupling and the mostly incoherent hopping mechanism (observed by inspecting non-perturbative FE-SH trajectories, see e.g. Figure 3 main text).
The rate expression used in Ref. 1 is very similar to the Marcus-Levich-Jortner rate for electron transfer 28,43 . It introduces an effective quantum high-frequency mode that is coupled to the exciton transfer process in order to effectively account for quantum-mechanical vibrations. In fact, in the common Marcus expression all vibrational modes in the system are treated as classical harmonic oscillators. While such an approximation is generally good for treating low frequency motions of a general solvent, it is not completely accurate for intra-molecular modes of the molecules, since for these modes ω > k B T . The approach employed by Aragó et al. gives a diffusion constant in the same order of magnitude as FE-SH and BC-Ehrenfest approaches and slightly longer diffusion length compared to both non-adiabatic propagation schemes. This can be explained with the fact that high frequency vibrations play a non-negligible role for the transport and, for systems with a reasonably high activation barrier (e.g. ANT), a classical treatment of modes might not be entirely justified (even at room temperature). This fact provides an additional explanation for the lower diffusivity predicted by the FE-SH and BC-Ehrenfest approaches, which treat all the modes with classical force fields, compared to both experiment and Marcus-Levich-Jortner like expression. On the other hand, the excitonic couplings computed in Ref. 1 were found by using a diabatization scheme that, to the best of our understanding, includes only 2 adiabatic states 1 .
Therefore the coupling between Frenkel exciton states is possibly an effective coupling still contaminated by mixing with charge-transfer states. For this reason, these couplings are slightly larger compared with both TrESP and TD-DFTB couplings (which by constructions do not include such a mixing with charge transfer states).
Overall, considering the differences in simulation set-ups, parameters entering the various Hamiltonians or analytic expressions and the uncertainty commonly characterizing experimental measurements, we show in Supplementary a6T The exciton diffusion length of a thin-film made of a6T molecules vacuum-deposited on quartz was reported to be about 60 nm by the quenching of the photoluminescence 44 . This value is close to the experimental diffusion length reported for anthracene, although for the latter system the exciton lifetime is an order of magnitude longer (10 ns) than that in a6T (1.8 ns) 45 and compensate for the smaller diffusion constant in ANT compared to a6T (see Table 2

DCVSN5
To the best of our knowledge, DCVSN5 experimental measurements for the diffusion length are not available. However this system has been simulated in Ref. 1

Y6
Also in the case of Y6, the comparison between the experimental estimates of diffusion constant and diffusion length and FE-SH in Table 2 of the main text should be taken just as indicative. In fact, Y6 forms a thin-film in real devices. In order to measure the exciton diffusion length, authors in Ref. 16 measure the exciton-exciton annihilation rate γ and calculate the exciton diffusion constant D according to γ = 4πDR 0 . The above equation can be derived through solving the 3D diffusion equation. In the deduction, the 3D dimensionality is embodied in the factor of 4π. In general the diffusion length is dependent on the system dimension and it can be defined L = √ nDτ . Where n = 2, 4, 6 for 1D, 2D and 3D diffusion respectively. Experimentally however, in Ref. 16 , L is defined as L = √ Dτ . Therefore in Table 2 of the main text we provided values from both definitions of the diffusion length applied to a 3D model.

Supplementary Note 7 Comparison exciton vs charge transport
In Supplementary Table 6 we compare transport parameters characterizing hole and exciton transport for some of the organic crystals investigated in this work. We report hole transfer interactions for the same crystal pairs in Figure 1 (main text) and hole reorganization energies and compared these with excitonic interactions (V ) found using MS-FED-FCD and related reorganization energies.
The hole transfer couplings have been computed using two different approaches. On the one hand, we used the well-established fragment orbital density functional theory (FODFT) approach Ref. 34,35,47 . This method is related to the interaction between the diabatic state wavefunctions constructed from orbitals of the isolated donor and acceptor fragments (similarly to what happens in constraint density function theory (CDFT), but unlike in CDFT, the interaction between donor and acceptor fragments is neglected).
FODFT calculations are carried out with the CPMD program package using the PBE exchange correlation functional as previously done in Ref 34,35,47 . We also used a projection operator-based diabatization (POD) approach as Ref. 48  pointed out before 35 that the coupling fluctuations actually enhance the delocalization of the states and help the transport in the hopping regime, while they become detrimental in the band-like regimes as they act as scattering sites. This observation give an additional explanation for the fact that exciton transport in the investigated systems (where there is a reasonably high barrier) is slower than the corresponding charge transport (see Figure 5 in the main text).
The most important difference between hole and exciton is the magnitude of charge reorganization energies as compared to the corresponding exciton reorganization energies. We note that, for comparison, hole reorganization energies were computed using B3LYP functional, which is the most common and widely used functional for electron/hole reorganization energies 3,14,50 because of its good agreement with UPS spectroscopic measurements. 51,52 Taking as an example the anthracene molecule (see Supplementary Figure 14), the hole reorganization energy is ≈ 140 meV while the relaxation of the exciton is about four times larger. This is due to the large charge redistribution occurring upon excitation. In particular, just by looking at the NTOs, it is evident that almost every bonding interaction in the HOMO becomes an anti-bonding interaction in the LUMO and vice versa. This causes each bond to become stronger (in case of a bonding interaction) or weaker (for an anti-bonding one) changing the length of the actual bond and causing a large structural modification upon excitation (see also Supplementary Figure 14). These sizeable reorganization energies, that in FE-SH are used to re-parametrize the classical FF as explained in the main text, are responsible for the strong exciton localization in these systems. Bond The + and − signs indicate an increase and decrease in bond length going from the ground (neutral) to the excited (charged) system for exciton (charge), respectively. The displacements are used to parametrize the force fields in FE-SH and FOB-SH respectively. Note that in this figure, bond displacements for ANT hole are taken from Ref. 34 and refer to B3LYP/6-311G(d) basis set. The NTOs related to the excitation of ANT are also shown.
Having previously analysed both charge and exciton transport parameters, we show in Supplementary   Figure 15 a representative surface hopping trajectory for both transport cases within the herringbone plane of the anthracene crystal. The most striking difference is related to the transport mechanism of charge and energy carriers. Indeed, while the excitonic wavefunction is completely localized on a single molecule and it is mostly hopping from one molecule to a neighbouring one, the corresponding charge carrier wavefunction is much more delocalized ( IPR ≈ 5) and it is characterized by extended diffusive jumps as described in Ref. 35 for other high mobility organic crystal. The thermally accessible valence band states are more delocalized in ANT than the corresponding excitonic states, and this allows the charge carrier wavefunction to extend to several molecules (more than 15 in the example trajectory of Supplementary Figure 15) and to cover a longer spatial displacement compared to the localized exciton.
This is a consequence of the more efficient transient delocalization. As we already showed in Figure 3 in the main text, a larger wavefunction delocalization correlates with a faster transport and a larger diffusion constant. In the particular case of ANT, the diffusion constant is much larger in the case of hole transport than to the corresponding exciton diffusion constant (D b =0.0905 cm 2 s −1 vs 0.0033 cm 2 s −1 , respectively for similar system sizes). Figure 15: Charge vs Exciton transport in ANT Hopping mechanism of exciton (a) and transient delocalization mechanism of electron hole transport (b) in ANT. See caption of Figure 3 main text for an explanation of the plotted quantities.