Quantum dynamics simulation of intramolecular singlet fission in covalently linked tetracene dimer

In this work we study singlet fission in tetracene para-dimers, covalently linked by a phenyl group. In contrast to most previous works, we account for the full quantum dynamics of the combined excitonic and vibrational system. For our simulations we choose a numerically unbiased representation of the molecule's wave function enabling us to compare with experiments, exhibiting good agreement. Having access to the full wave function allows us to study in detail the post-quench dynamics of the excitons. Here, one of our main findings is the identification of a time scale $t_0 \approx 35 \text{fs}$ dominated by coherent dynamics. It is within this time scale that the larger fraction of the singlet fission yield is generated. We also report on a reduced number of phononic modes that play a crucial role for the energy transfer between excitonic and vibrational system. Notably, the oscillation frequency of these modes coincides with the observed electronic coherence time $t_0$. We extended our investigations by also studying the dependency of the dynamics on the excitonic energy levels that, for instance, can be experimentally tuned by means of the solvent polarity. Here, our findings indicate that the singlet fission yield can be doubled while the electronic coherence time $t_0$ is mainly unaffected.

(*Electronic mail: sam.mardazad@physik.uni-muenchen.de) (Dated: 17 November 2021) In this work we study singlet fission in tetracene para-dimers, covalently linked by a phenyl group. In contrast to most previous studies, we account for the full quantum dynamics of the combined excitonic and vibrational system. For our simulations, we choose a numerically unbiased representation of the molecule's wave function, enabling us to compare with experiments, exhibiting good agreement. Having access to the full wave function allows us to study in detail the post-quench dynamics of the excitons. Here, one of our main findings is the identification of a time scale t 0 ≈ 35 fs dominated by coherent dynamics. It is within this time scale that the larger fraction of the singlet fission yield is generated. We also report on a reduced number of phononic modes that play a crucial role in the energy transfer between excitonic and vibrational systems. Notably, the oscillation frequency of these modes coincides with the observed electronic coherence time t 0 . We extend our investigations by also studying the dependency of the dynamics on the excitonic energy levels that, for instance, can be experimentally tuned by means of the solvent polarity. Here, our findings indicate that the singlet fission yield can be doubled while the electronic coherence time t 0 is mainly unaffected.

I. INTRODUCTION
Singlet fission (SF) is a spin-allowed photophysical process that generates two triplet excitons from one singlet excited state [1][2][3][4] . It can realize multiple exciton generation (MEG) and provide longer-ranged exciton transport inside a semiconductor (typically of the order of µm for triplets compared to ∼ 10 nm for singlets 5,6 ). This process gives rise to an increase in the charge-carrier to radiation ratio, resulting in a net increase in the efficiency for organic semi-conductor based photo-voltaics. SF has therefore shown the potential to surpass the upper bound of the solar cell efficiency of 33 %, set by the Shockley-Queisser limit and, hence, has attracted great attention 7,8 .
Previous theoretical and experimental studies 9-21 investigated the mechanism behind SF and distinguished between coherently driven and thermally activated SF, where only the latter exhibits a strong temperature dependence. The different mechanisms happen on different time scales namely up to 700 fs for coherently driven SF 17 and > 1 ps for thermally activated SF 10 . Unlike widely studied intermolecular singlet fission (xSF) systems of crystalline tetracene and pentacene, covalently linked chromophore dimers allow much easier control of inter-chromophore orientation and interaction for intramolecular singlet fission (iSF) 10,22 . This is achieved by the chemical modification of the bridging group and tuning of the solvent environment. For example, Wang et al. recently reported fast iSF within 10 ps in 1,4-bis(11-phenyltetracen-5-yl)benzene, as is shown in Figure 1, and illustrated its temperature dependence and significant solvent polarity effects 9 .
It is nowadays well-known that a pure electronic description is insufficient to describe these processes correctly 2,3,23-25 , suggesting the necessity of the incorporation of molecular vibrations into theoretical considerations. The main obstacles for realistic calculations are, besides others, a faithful microscopic modeling of the coupling between electronic states and vibrational degrees of freedom, and a full quantum mechanical treatment of their coupled dynamics. Previously, due to the lack of knowledge of the full coupling Hamiltonian, the first problem has been addressed by investigating the effect of successive elimination of vibrational degrees of freedom in order to identify the relevant modes 11,12 . To accomplish a faithful simulation of the quantum dynamics following the light excitation, a variety of methods have been applied, including multi-configurational time-dependent Hartree (MCTDH) [26][27][28] , Redfield theory 29,30 and time-dependent wave packet diffusion (TDWPD) 31 . However, despite numerous investigations, there have been few numerically unbiased simulations of a full microscopic Hamiltonian that is generated from first principles for both, electronic molecular and vibronic degrees of freedom 12,32 . Common timeevolution schemes, e.g. Redfield theory, are only valid in the weak-coupling regime. To overcome such limitations, a decomposition of the vibrational degrees of freedoms into clusters 11,12 , e.g., multi-layer-MCTDH 33,34 , can be introduced, which is superficially related to the tensor-network approach presented in this work.
Recently, Chin, et al. applied tree tensor-network states (TTNSs) 35 to model xSF in 13,13'bis(mesityl)-6,6'-dipentacenyl (DP-Mes) dimers. Exploring an involved numerical procedure, they effectively computed the non-perturbative real-time dynamics of exponentially large vibronic wave functions of real molecules. In this work, we attempt to adopt another novel one-dimensional tensor network method to simulate iSF dynamics in realistic chemical systems. Thereby, our main goal is to introduce a conceptually simple yet powerful and numerically well-controlled procedure, as described in the following. First, we compute the spectral density from ab initio multireference quantum chemical calculations and determine the effect of the molecule's vibrational modes on the electronic overlap integrals to obtain an effective model, capturing the coupling between electronic states and vibrations. Employing the resulting microscopic model, we perform a simulation of the full post photo excitation quantum dynamics using a novel method from tensor networks 36 .
In particular, this method allows for the description of many bosonic degrees of freedom with a large number of internal states by dynamically selecting the relevant bosonic modes. The resulting, advanced numerical simulations enable us to achieve good agreement with experimental data measured for the tetracene para-dimer, covalently linked by a phenyl group. Furthermore, we are able to unambiguously identify the relevant vibrational modes dominating the post-excitation dynamics by studying the energy transfer within the molecule. Our findings of the coherent and incoherent dynamics' time scales are in very good agreement with previous investigations 11,12 . However, here, these results are based on a realistic microscopic modeling of the full molecule's dynamics. Having access to the molecule's full wave function, we study the connection between energy transfer into the vibrational system and phonon-assisted, coherently driven iSF. Our analysis reveals a competition between an initial coherent dynamics generated by the vibrational system and its suppression due to the formation of heavy exciton-phonon quasi particles. Finally, being able to perform realistic simulations, we tuned experimentally relevant parameters such as the solvent polarity to identify the parameter regime generating the largest SF yield.
The scope of this paper is as follows: In section II, we introduce a Hamiltonian for excitonphonon mixtures describing our system adequately. This employs a global coupling of non-local vibrations to excitonic degrees of freedom. In section III, we briefly sketch the new method and exhibit performance advantages. Furthermore, we are going to discuss the necessity for large local Hilbert spaces in simulating these types of models. Accounting for the rapid development of tensor network techniques for large Hilbert space sizes, we also briefly elaborate on why standard, widely used methods may not be appropriate here and discuss our approach in the context of the TTNS formulation introduced by Chin et al. 35 . Finally, we present the simulated iSF dynamics in section IV. These involve comparison to experiment, as well as detecting the phononic modes dominating the energy transfer. We find a total triplet population between 14 % and 25 % depending on the solvent, comparable to the order of magnitude of experimental findings, i.e., 21 % 9 .
Eventually, we investigate the tuning of experimental parameters, i.e. the solvent polarity, on the total triplet yield.

II. MODEL
To model the SF process, we describe the diabatic electronic states coupled to vibrational modes by means of a Frenkel-exciton-Hamiltonian 37,38 , where the lower-case indices i and j correspond to the excitonic states and the upper-case indices I denote the vibrational modes. Here,ĉ i correspond to the usual excitonic creation and annihilation operators, whileb ( †) I create and annihilate phonons. The diagonal elements (i = j) of V i j resemble the energies of the excitonic states, while the off-diagonal elements (i = j) correspond to the couplings among them. ω I is the frequency of the vibration mode I, and g i j,I represents the coupling between the excitonic Hamiltonian and vibrational modes I. We can distinguish between diagonal and off-diagonal couplings here as well. In order to describe our system adequately, we use five diabatic states. They are the two locally excited (LE) states LE 1 and LE 2 , the two charge transfer (CT) states CT 1 and CT 2 , and the triplet pair (TT) state. While the LEs model the electronic excitation in the dimer subsystems, the CTs are energetically high lying states that, if occupied, represent intermediate cationic or anionic molecular states 39 . By initializing the system in one of the LE states (or a coherent superposition), we model a photo-excitation. It was previously shown that the localized triplets in the chromophore units can be interpreted as an electron-hole bound state 4 . However, these triplets in each subsystem are correlated into an overall singlet state 2,40 , which we call TT. The dissociation of the TT state into locally separated triplets is a topic on its own and not part of this investigation 39,41 . The details of the construction of these five diabatic states from state averaged complete active space self consistent field (SA-CASSCF) excited state calculations can be found in appendices A and D. The diabatic exciton matrix elements V i j thus obtained are shown in Table I. Note that the multi-exciton TT state is optically dark 4,42 . Furthermore, the coupling between LE and TT is almost zero, and, thus direct conversion is strongly suppressed. This originates from the fact that the direct process corresponds to a double-electron transfer. The related two-electron integrals, between HOMOs/LUMOs of different groups or molecules, are usually vanishingly small 3,16,43 . However, the couplings between CT and TT or LE are relatively large (40 −60 meV). This suggests a possible indirect superexchange iSF path (LE→CT→TT) [44][45][46][47][48] . It should be mentioned that due to the lack of dynamical correlation in the CASSCF calculations, the diagonal entries of the excitonic matrix V i j may be inaccurate. Therefore, we tune the LE and TT energies to fit the peak position of experimental absorption spectra. Furthermore, the energies of CT are obtained from a linear search within the range of 3.0645 −3.5645 eV using Redfield dynamic simulations 30 . The details about this procedure can be found in appendix D.
To estimate the coupling strength between excitonic states and vibrations we calculate the spectral density according to and the temperature-dependent fluctuations of the excitonic Hamiltonian elements using The details of the methodology for evaluating the exciton-phonon couplings and the graphical illustrations of some essential vibrational modes can be found in appendices B and C. The results are shown in Table I and Figure 2. From the selected spectral density displayed in Figure 2 elements. Under this oscillation, the tetracene group remains unaltered, and therefore, the diagonal elements are less sensitive to this mode. The off-diagonal elements, however, depend on the overlap between orbitals from different tetracene chromophore units. The latter are dominated only to a small part by the orbitals in the bridge because the tetracene groups are far away from each other. Consequently, twisting the bridge will change only the off-diagonal elements significantly.
Furthermore, the lowest frequency mode around 7 cm −1 also contributes to off-diagonal elements because it changes the dihedral angle between the planes of two tetracene units.
In order to maintain numerical feasibility we have to pick the relevant vibrational modes of the 88 atomic molecules. Out of the 258 vibrational modes we neglect 181 that trigger relative fluctuations smaller than 0.1 %. We also discard the modes with frequencies below 10 cm −1 since their oscillation periods are longer than 200 fs. This time scale is set by the time-range of our quan-tum dynamics simulation. Finally, we couple the remaining 76 modes, ranging from frequencies 10.18 cm −1 (0.0013 eV) to 1714.2 cm −1 (0.2125 eV), to the excitonic system.

III. METHODS
Simulating the full quantum dynamics of Equation (1)  We overcome these issues by employing the newly developed method of projected-purification (PP) 36,55 . Within this approach, we restore the U(1)-particle number conservation symmetry by doubling the size of the (phononic) system. Furthermore, we introduce a dynamical procedure to truncate the local phononic Hilbert space dimensions d adaptively by the choice of our truncated weight (δ = 1.0 × 10 −8 ). This enables us to treat large local Hilbert spaces with n ph = 63 very efficiently with the approximation being well-controlled by the maximally allowed sum of discarded Schmidt values.

A. Representing the vibrational modes
Let us briefly summarize the key points of this procedure. We start with a wave function expanded in a local basis labelled by irreducible representations of a global operatorN Here, L corresponds to the total number of single particle orbitals and the n i run from {0, . . . , d − 1}. Exploiting global U(1) symmetries typically provide a significant speedup in numerical sim-ulations 56 , however, the number of phonons is not conserved in Equation (1). We circumvent this problem by introducing an auxiliary environment Hilbert space, which is a copy of the original one, Here, pure indices n j correspond to physical degrees of freedom, while those carrying a barn j involve the auxiliary space. We further introduce the gauge condition This transforms Equation (5) into which happens to be an eigenstate ofN +N, i.e., the total particle number operator expanded to the enlarged Hilbert space, with particle number Ln ph . Figure in terms of the operators breaking the global U(1) symmetry. The remaining operators are tensored with identities in the auxiliary space. With this representation ofĤ, one always stays within the purified sub-manifold 36 .
Adopting this representation, we have obtained a formal global U(1) symmetry for the phonon system that is not particle number conserving, at the cost of introducing a reservoir degree of freedom for each bosonic orbital, as shown in Figure 3b. Instead of treating a few huge matrix blocks, we can now compute with a lot of small ones. We exploit the potentially large number of phys aux phys aux Here, the index n indicates the singular values belonging to the corresponding irreducible representation. Discarding the smallest singular values fulfilling ∑ n,τ (s n τ ) 2 ≤ δ corresponds to eliminating occupation number elements ρ nn not required for an adequate description of the state. Hence, the initial local Hilbert space dimensions d of the vibronic modes are truncated to significantly smaller values, rendering the simulations numerically feasible. During the dynamics, the truncation process of the state enforces the system to keep only the modes required for an adequate approximation of the state (i.e., with sufficient weight).
The restored global U(1) symmetry of the phononic system combined with this truncation scheme results in a highly diminished memory requirement and, more importantly, reduction in CPU time by a factor of 100-1000, due to decomposition and parallelization respectively, as is shown in Figure 3c. We compare the total CPU times during the initial stage of the dynamics, following a photo excitation of a bright exciton with and without exploiting the PP mapping and for different local Hilbert space dimensions d in the case of the trivial phonon representations. Note that in the conventional approach, the memory requirements grow rapidly, causing a breakdown of the simulations without the PP mapping after times of only 8.5 fs or 9.9 fs, respectively. Treating the exciton-phonon coupling in a numerically controlled manner is crucial in order to describe the dynamics of the vibrational modes correctly. For that purpose, we monitor the probability of different phonon occupation numbers by calculating the RDM during the time evolution.

B. Numerical stability and effects of insufficiently large vibrational Hilbert spaces
The RDM with respect to a vibronic mode J is obtained by tracing out all degrees of freedom from the total density operator except for the target mode J, where {I} is the set of all phononic modes. Exemplarily, we investigate the dynamics of a partic- This means that the reported deviations contribute overall errors in the time evolution scaling as therefore, we are able to truncate away unnecessary modes without any additional numerical effort.

C. Relation to other tensor network techniques
Finally, we comment briefly on why two other common tensor network methods appear unsuitable for a full quantum dynamics simulation of the current problem (for a more detailed comparison see 55). We do not discuss MCTDH here since it was extensively explained in previous works 11,26,27,[32][33][34] . is easily adopted to various setups. Based on the previous discussion, we believe that combining the PP mapping with advanced entanglement renormalization could provide a further big leap.

A. Absorption spectra and energy transfer
In order to verify the applicability of our model and methods, we evaluate the absorption spectrum, 66 Here, C(t) is the auto-correlation function, and τ is an appropriately short effective relaxation time. |ϕ(0) is the initial excitonic state, and |χ(0) is the vacuum phononic state. Our system was prepared in the bright initial state configuration |ϕ(0) = (|LE 1 + |LE 2 )/ √ 2. An analogous analysis for the localized and dark initial state can be found in appendix F.
The resulting calculation is displayed in Figure 5a. In order to resemble the relative magnitude of the 0-0 and the 0-1 peaks of the experimental spectrum 9 , the relaxation time is chosen to be τ = 230 fs. We find the theoretical data capturing the main experimental features for the vibra-tionally resolved absorption spectrum. In particular, the position and relative magnitudes of the several dominating vibronic peaks are in excellent agreement. Furthermore, small wavelength features, e.g., a shoulder and a small peak (around ∼ 380 nm and ∼ 420 nm, respectively), are found to be correct. However, their absorption probability is slightly overestimated. We attribute the experimentally observed increased signal widths to the neglected phononic modes, which generate additional scattering processes between diabatic states and vibronic excitations.
In Figure 5b we see the dynamics of the populations of the diabatic states. Note that TT has a small delay time vs CT. This indicates that the indirect mechanism dominates the process. Both the populations of LE and CT show oscillations whose frequency is about 0.12 fs −1 . We find this value corresponding to the energy gap between two eigenstates of the pure excitonic Hamiltonian.
These two eigenstates are accounting for over 96 % of LE and CT when expanding the eigenbasis in terms of the components of the excitonic basis, i.e., they carry the main weight. This demonstrates how far the driving force from excitonic coupling is important for iSF in this system on short transient time scales, in agreement with our findings from section IV B. Note that the coherent oscillations disappear in the scope of the dynamics. This implies that the vibrational system is chosen large enough to achieve relaxation into an incoherent excitonic dynamic, contrasting previous works 11,12 .
Having access to the full wave function, we are able to study the energy transfer between the excitonic and phononic systems. For that purpose, we evaluate the partial energies of the exciton and phonon system, as well as the coupling energy by separating Equation (1). The obtained time dependencies are displayed in Figure 5c. Most importantly, we find only few dominating oscillations in the coupling and vibrational energies at finite values, after an initial transient regime.
As shown in Figure 5d, a Fourier transformation reveals that there are two frequencies dominating the energy conversion between the exciton-phonon coupling and the phononic subsystem. They are located at around 1420 cm −1 and 1620 cm −1 , respectively. We identify these signals with vibrational modes as follows:

B. Solvent and triplet energy dependence
In a subsequent stage, we investigated the dependency of the triplet production on the energies of CTs and TT. Importantly, the CTs on-site energies can be tuned in the experiment by using various solvents with different polarities 68 . For instance, Alvertis et al. 69 recently revealed the significant energy splitting of CT states of the orthogonal tetracene dimer (DT-Mes) in polar solvents and its vital role in the switching between coherent and incoherent iSF. However, because the inter-chromophore couplings are much smaller in our covalently linked tetracene dimer system than in DT-Mes, due to an additional phenylene group between two tetracene units, the local electric field by polar solvents will not result in significant differences between the two CT states in our case. Therefore, the energy splitting between the two CT states is neglected in this work.
As shown in Figure 6 we scanned a whole parameter region shifting the energies. We measure the offsets of the CTs and TT from the values given in Table I by the gaps ∆E(CT − LE) and ∆E(CT − TT), respectively. Most prominently, we find that lowering the energy of the CT states, i.e., increasing the solvent polarity, yields a larger TT population rate at the early stage. Moreover, the gained triplet population can be elevated by increasing the energy of the TT state. Combining these observations, we conclude that one dominating process for the TT yield is the relaxation of CT populations into the energetically closest diabatic states, suggesting the indirect iSF mechanism (LE→CT→TT). This picture is consistent with the observation that the TT production rate is strongly suppressed in the lower right corner of Figure 6a. In addition, we note that the triplet yield obtained from the system parameters (Table I)  slope of the TT production rate changes independently of the specific choice of the on-site potentials. As discussed in appendix G, this can be addressed to a change in the nature of the excitonic dynamics. At times t < t 0 , all excitonic states realize a coherently driven dynamics. However, at times t > t 0 , this mechanism changes into a classically driven relaxation of mixed states of the LEs, CTs, and TT subsystems. Previously, the coherent increase in the TT population has been interpreted in terms of phonon-assisted hoppings employing Fermi's golden rule 11 . However, we also found that the electronic coherence time roughly coincides with the oscillation time of the dominating vibrational mode nos. 184 and 185. In fact, their wavenumber of 1420 cm −1 roughly translates to a period of 25 fs. Noting also that ∼ 2/3, i.e., the dominating part of the TT population is generated during the coherent dynamics, this could create a competition between phonon assisted hoppings and the time scale in which these processes are possible. In appendix G, we investigated the renormalization of the excitonic couplings after rewriting the Hamiltonian using a Lang-Firsov transformation 72 . Importantly, we find that the renormalization of the couplings between the TT state and the remaining excitonic system changes from an enhancement to a suppression after t 0 = (34.725 ± 0.165) fs. It is a remarkable observation that the description of the phonon-assisted hopping in terms of polaron dynamics yields a perfectly matching prediction of FIG. 7. CT-TT hopping enhancement with the schematic sketch of quasi-particle formation. The coefficients λ i j,I = g i j,I /ω I define the ration between the exciton-phonon coupling and the vibration frequencies.
the electronic coherence time. In this quasi-particle picture, the coherent phonon-assisted hopping processes are blocked by the momentum transfer into the vibrational system. In fact, this momentum transfer creates a cloud of vibrations accompanying the excitonic states 61 , and thereby increases the effective exciton mass. Our analysis further shows that within the coherent dynamics, the dominating hopping is between CT and TT states, pointing out the importance of the indirect hopping mechanism. A schematic representation of the different regimes and the corresponding type of quasi-particles governing the dynamics is shown in Figure 7.

V. CONCLUSION AND OUTLOOK
We studied the TT production rate via iSF for a covalently linked tetracene dimer using a largescale full quantum-mechanical treatment of the relaxation dynamics after photo-excitation. Our study extends the theoretical analysis conducted so far 11 by treating couplings between the diabatic states and the vibrational modes obtained from first principles. Incorporating a large number of molecular vibrational modes, we achieve a realistic description of the full molecule's dynamics, following photo-excitation. We obtain excellent agreement with experimentally observed absorption spectra and identify only a small number of vibrational modes, dominating the energy conversion between the different subsystems. Analyzing the renormalization of the exciton couplings in terms of exciton-phonon quasi-particles, we can trace back the origin of the phonon-assisted coherent TT generation to an enhanced delocalization of the created quasi-particles, which after some time transform into heavy, localized ones. A systematic variation of molecular parameters that can be related to the solvent polarity allows us to identify those regions in parameter space exhibiting the largest triplet yield.
As it is of special interest to raise the TT population to produce a high charge carrier yield, our findings can be employed to identify promising manipulations of the molecular structure. Here, we suggest two pathways. First of all, decreasing the energy gaps between CT and LE or TT is crucial for producing the higher TT yield. Our results indicate that both modifications of CT and TT energies mainly affect the overall scale of the triplet yield. Here, during the initial coherent dynamics, a quick increase in the TT is observed which transforms to a classically driven relaxation, happening on time scales t t 0 = 35 fs. We note that the observed time scale for coherent dynamics is in agreement with previous studies, both theoretical and experimental [9][10][11][12][13][14][15][16][17] . However, this time scale is very close to the oscillation period ∼ 25 fs of the vibrational mode nos. 184, 185, and 186 dominating the energy transfer between the excitonic and vibrational systems. Therefore, our findings suggest a second route that, though experimentally more challenging, would be an increase in the electronic coherence time t 0 . Noting that the dominating fraction of the TT yield is generated from the initial, coherent dynamics, a combination of larger coherent times, i.e., a reduction in the frequency of the oscillations of the overall molecule's backbone with an increased solvent polarity could, generate significantly larger TT occupations.
There are still many aspects to be further considered for future simulations of SF dynamics.
First of all, we might take into account more chromophore units. It was suggested recently that the spatial separation of the correlated triplets might be a crucial point in the harvest of a large yield 9 .
Since the presented numerical framework can adaptively reduce the number of vibrational modes, it might be a promising tool to explore the dynamics of larger molecules, too. Additionally, finite temperature effects, which may become relevant for larger molecules 9 , can be incorporated by a thermofield approach, which is a standard technique for tensor network methods. From a methodical point of view, a direct comparison between established time-evolution schemes such as multi-layer MCTDH 33,34 and the presented tensor network representation would be desirable. In particular, a faithful investigation of these different approaches could yield insights which method should be used to optimally describe SF dynamics. Furthermore, we believe that a combination of our representation of the vibrational degrees of freedom with TTNS and entanglement renor-

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Competing interests
The authors declare no competing interests.

Appendix A: Electronic Structure Calculation
We employ CASSCF with a 6 − 31G(d) basis set in order to calculate all the excited states involved in our iSF model due to the multi-excitation nature of TT. The four orbitals near highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) are delocalised. Therefore we use the Pipek-Mezey method 74 on these orbitals to generate four localised ones (from HOMO-1 to LUMO+1). Each of them is confined in single tetracene chromophore.
We construct all the diabatic states from the four localised orbitals and calculate all the excitonic Hamiltonian elements V i j using complete active space configuration interaction (CASCI). All calculations are performed using the OpenMolcas package 75 .
Here we only take into account the linear exciton-phonon coupling terms. We distort the whole molecule (at equilibrium geometry) along selected vibrational modes' displacement for ±0.01 Å.
Then we calculate the excitonic Hamiltonian with the same procedure as mentioned above. With the Hamiltonian at three distinct points in space (±0.01 Å,0 Å) , we can perform linear fitting of each elements V i j and calculate the gradient. Note that the coupling g i j,I in Equation (A1) is connected to dimensionless displacement Q I . One has to multiply the gradient by a factor of h/ (m I ω I ).

Appendix C: Spectral Density Data
Spectral densities and fluctuations for some excitonic Hamiltonian elements which were not included in the main text are displayed in Figure 9.

Appendix D: Redfield search
To obtain a reasonable shifted CT energy, we use Redfield theory 29,30 to simulate the SF dynamics at the temperature of 300 K. In Redfield theory, the dynamics are described using where ρ is the density matrix of the excitonic states, and α, β , γ, κ label eigenstates of the excitonic Hamiltonian. Furthermore the Redfield tensor R αβ γκ is defined as ε is the energy of the corresponding eigenstate and C j = δ ε j (t) δ ε j (0) is the energy fluctuation correlation function of the state j. We then tune the energy of CT within the range of 3.0645 −3.5645 eV, according to the CASSCF results and the previous solvent effect investigations 9 , while fixing the energy of the other excitonic states. In order to perform the Redfield dynamics we choose the bright excited state as was discussed in the main text. Since E CT = 3.3645 eV has the highest TT population we choose this as the TDVP input.

Appendix E: Numerical control parameters
In order to convince ourselves of the validity of the simulations we plot the maximum bond dimension and norm deviation in Figure 10.  TDVP is an inherently local method which is only exact for nearest-neighbour Hamiltonians or infinite bond dimension 51 . Otherwise a so called projection error is introduced, which causes the time evolution to be non-unitary. In order to estimate this error one can look at the deviance from the initial norm as in Figure 10. As we can recognise, the total norm difference is constant throughout the simulation. Furthermore, this is in accordance with the truncation error. However, it is indeed initial state dependent, i.e. the bright exciton initial state has twice the norm error of the localised one while the dark initial state actually lies one order of magnitude higher. This is explained with regard to the bond dimension, i.e. the phononic modes needing high bond dimensions. Since the TDVP projector is exact for an infinite bond dimension, i.e. for an exact state 51 , the stronger the truncation of the state, the stronger the breaking of unitarity.

Appendix F: Initial state dependence
We investigate the dependency of the excitonic population dynamics on the initial state. To do so, we perform the time evolution with the excitation localised in one of the chromophores and with an anti-symmetric excitation of both chromophones, as can be seen in Figure 11. It is to mention that the anti-symmetric initial state is dark, i.e. it is not realisable in experiments from the ground state, due to transition rules. We see that for LE initial state |ψ (0) = |LE 1 the dynamics are similar to the bright initial state in the main text, i.e. there is also a delay from the CT. However the SF yield is around half of the yield of the bright initial state.  Figure 11a corresponds to the localised initial state while Figure 11b is the dark one. Both were carried out with E TT = 2.9493 eV.
As it was argued before, the mechanism behind the coherent SF is the superposition of the LE and TT state 70,71 . This is the reason for the small yield of the localised state. It is remarkable that a delocalisation over two chromophores gives a factor of two in the TT yield. This points into the direction of previous investigations 9 emphasising the importance of spatial delocalisation.
Using the dark initial state we see that a much higher yield than for the system in the main text is reached. After a plateau-like region the value approaches ∼ 10 %. However, there is a maximum after which the population begins to decay back mostly into the LE states.

Appendix G: Excitonic RDM analysis
In the main text we already discussed the appearance of two different scales characterising the time evolution of the TT occupancy. In order to analyse the origin of these time scales we studied the RDM of the excitonic system. In general, the time evolution of the excitonic density operator ρ ex after tracing out the vibrational degrees of freedom can be expanded in its eigenbasiŝ ρ ex (t) = ∑ i, j∈S i|ρ ex (t)| j |i j| = ∑ n λ n (t) |λ n (t) λ n (t)| . (G1) Here, we explicitly exhibited the time-dependence of the eigensystem and defined S = {TT, CT 2 , CT 1 , LE 2 , LE 1 }. In Figure 12 we show the absolute value of the decomposition of H TT (t) = − |TT log λ 3 (t) TT| (G5) Note that due to the normalisation of ρ ex we have ∑ n λ n = 1. In this decomposition of the excitonic density operator (at times t > t 0 ), the time evolution of the CTs and TT states is described by the dynamics of a purely mixed state without quantum mechanical coherences. Furthermore, the dynamics of the LEs w.r.t. the remaining excitonic states is of the same type and coherent dynamic only happens within the LEs. In order to investigate the origin of the found electronic coherence time t 0 we rewrite equation Equation (1) For each exciton-phonon coupling matrix Λ I we computed the eigenvalues which we found to be strictly positive. Additionally, each element off † IfI is between zero and one. These observations motivate an expansion of the hopping elements of the transformed Hamiltonian to first order in the generator coefficients where we explicitly used the fact that we are working in a single-exciton subspace. Here, we also introduced the hopping energy renormalisation ∆Û jl , which is referred to in Figure  We find a change of sign at a simulation time t 0 = 34.725 fs. Note that for t < t 0 we observe a negative sign in the dominating matrix elements corresponding to hopping from the CT into the TT state indicating an enhancement of the transition amplitude. For t > t 0 the situation changes into the opposite and hoppings into the TT state are suppressed. Importantly, we also find a nearly vanishing renormalisation for the direct path LE −→ TT. This emphasises the significant role played by the super-exchange path in the TT production during the coherent dynamics.