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

Singlet fission (SF) is a spin-allowed photophysical process that generates two triplet excitons from one singlet excited state.^1–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¹⁷ and $>$1 ps for thermally activated SF.¹⁷ 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 Fig. 1, and illustrated its temperature dependence and significant solvent polarity effects.⁹ 

It is nowadays well-known that a pure electronic description is insufficient to describe these processes correctly,^2,3,24–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),^27–28 Redfield theory,^29,30 and time-dependent wave packet diffusion (TDWPD).³¹ 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 time-evolution 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, Schröder et al. applied tree tensor-network states (TTNSs)³⁵ 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.³⁶ 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 Sec. II, we introduce a Hamiltonian for exciton–phonon mixtures describing our system adequately. This employs a global coupling of non-local vibrations to excitonic degrees of freedom. In Sec. 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 Schröder et al.³⁵ Finally, we present the simulated iSF dynamics in Sec. 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%.⁹ Eventually, we investigate the tuning of experimental parameters, i.e., the solvent polarity, on the total triplet yield.

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, while $b̂I(†)$ create and annihilate phonons. The diagonal elements (i = j) of V_ij 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_ij,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₁ and LE₂, the two charge transfer (CT) states CT₁ and CT₂, 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.³⁹ 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.⁴ 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 the supplementary material. The diabatic exciton matrix elements V_ij 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).^45–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_ij 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.³⁰ The details about this procedure can be found in the supplementary material.

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 the supplementary material. The results are shown in Table I and Fig. 2. From the selected spectral density displayed in Fig. 2, we can see that for diagonal elements, there are high-frequency modes near 1400 cm⁻¹ dominating the coupling. The remaining spectral densities for other couplings can be found in the supplementary material. Mode nos. 184 (1409.61 cm⁻¹) and no. 185 (1411.18 cm⁻¹) account for these peaks. They change the backbone of the tetracene chromophore unit (see the supplementary material). This will influence the molecular orbital (MO) coefficients and integrals substantially and changes the energy of diabatic states consequently. At the same time, the lower frequency mode at 415 cm⁻¹, which mainly twists the plane of the bridging phenyl group (see the supplementary material), mostly contributes to the coupling in off-diagonal 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⁻¹ 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⁻¹ since their oscillation periods are longer than 200 fs. This time scale is set by the time-range of our quantum dynamics simulation. Finally, we couple the remaining 76 modes, ranging from frequencies 10.18 cm⁻¹ (0.0013 eV) to 1714.2 cm⁻¹ (0.2125 eV), to the excitonic system.

Simulating the full quantum dynamics of Eq. (1) requires a careful and well-converged treatment of the relevant vibrational modes. For that reason, we exploit a matrix-product state (MPS) representation^50–51 of the full quantum many-body wave functions to represent the time-evolution of the photo-excited system. The major challenge is the fact that each vibronic mode, in principle, requires an infinite number of internal degrees of freedom. For practical purposes, we need to truncate the local Hilbert space dimensions to a large, but finite, number d = n_ph + 1. In doing so, care must be taken to choose the local dimensions large enough such that the dynamics is well-behaved. Otherwise, the system will notice the upper limit of the allowed vibrational mode occupations. This truncation needs to be balanced with computational cost scaling as $Od3$ in the local dimension for the two-site time dependent variational principle (TDVP) method, which we used in our calculations to perform the time evolution.^52,53 The numerical calculations were carried out using the SyTen toolkit.⁵⁴ 

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 × 10⁻⁸). 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.

Let us briefly summarize the key points of this procedure. We start with a wave function expanded in a local basis labeled by irreducible representations of a global operator $N̂$⁠,

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 simulations;⁵⁶ however, the number of phonons is not conserved in Eq. (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 bar $n̄j$ involve the auxiliary space. We further introduce the gauge condition

This transforms Eq. (5) into

which happens to be an eigenstate of $N̂+N̄̂$⁠, i.e., the total particle number operator expanded to the enlarged Hilbert space, with particle number Ln_ph. Figure 3(a) illustrates this idea, showing the effect of the gauge constraints at the example of the two lowest excitations of a vibrational mode. It can be shown³⁶ that there exists a bijective mapping between the basis in the original Hilbert space and the basis of the doubled one. Note that the subspace spanned by all basis states $|n1…nL〉|nph−n1…nph−nL〉n1,…,nL$ has the same dimension as the original one. We also map arbitrary global operators to the purified Hilbert space, by substituting creation and annihilation operators with

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.³⁶ 

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 Fig. 3(b). 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 blocks by parallelizing contractions over the symmetry blocks. Furthermore, we make careful use of the real space parallel TDVP algorithm⁵⁷ for some of the solvent dependent calculations.

Importantly, the PP mapping allows us to variationally optimize in the space of retained optimal phononic modes by the tuning of one parameter, namely, the truncated weight during the singular value decomposition (SVD). This can be understood by means of the relation between the retained singular values s_τ during the decimation process of the density matrix renormalization group (DMRG)⁴⁹ and the diagonal elements of the reduced density matrix (RDM) of a single bosonic site ρ,³⁶ 

Here, the index n indicates the singular values belonging to the corresponding irreducible representation. Discarding the smallest singular values fulfilling $∑n,τsτn2≤δ$ 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 Fig. 3(c). 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 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. 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 particular phonon mode (no. 184), which plays a dominant role in the excitonic dynamics (see also Fig. 2 and the supplementary material). In order to clarify the role of the truncated local Hilbert space, we systematically tune the maximally allowed phonon numbers, i.e., n_ph = 7, 15, 31, 63. The result for the vibrational mode no. 184 is shown, in an exemplary fashion, in Fig. 4. It shall be mentioned that intermediate sizes, e.g., n_ph = 23, show the same behavior as the others and are therefore left out for clarity. We evaluate the RDM at three distinct times after photo-excitation, namely, 49, 99, and 148 fs. We observe that the dynamics of the probabilities to excite only a small number of $∼10$ vibrational quanta is correctly described for all maximally allowed phonon numbers n_ph. However, inspecting the tail of the excitation probabilities, significant deviations are visible at the order of δρ_J,nn>10 ∼ 10⁻⁵, as long as n_ph < 63. This translates into a failure of describing the relaxation of the initially excited phonon modes in the scope of the time evolution if the local Hilbert space dimensions are too small. We relate this problem to the contribution of the higher excited vibrational modes to the overall time evolution operator $e−iĤt$⁠. Here, the vibrational energy contributions of $Ĥ$ scale as $n̂J∝∑nn×ρJ,nn$⁠. Inspecting Fig. 4, we find that for n_ph < 63, the error is dominated by excitation probabilities of ρ_J,nn=10−30, i.e., an overall error of

This means that the reported deviations contribute overall errors in the time evolution scaling as $∼10−3$ in the case of mode no. 184.

These observations clearly show that the higher occupation numbers of the vibrational modes can make up for non-negligible contributions to the iSF quantum dynamics, even though the mean occupation and displacement of all vibrations can be quite small. This implies the necessity of employing a large, maximally allowed phonon number n_ph per lattice site. Throughout the simulations, the system dynamically chooses the actually required phonon Hilbert space dimension; therefore, we are able to truncate away unnecessary modes without any additional numerical effort.

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 Ref. 55). We do not discuss MCTDH here since it was extensively explained in previous works.^{11,26,27,33–34}

The pseudo-site (PS) DMRG unfolds the local Hilbert spaces of the phononic part of the system into auxiliary sites and encodes the phonon occupations into a binary representation.^58,59 Thereby, the computational complexity originating from the large local Hilbert space dimensions is reduced drastically at the cost of introducing long-ranged couplings. The PS DMRG has been applied successfully to study the out-of-equilibrium-dynamics of Holstein models.⁵⁹ However, in these systems, the electron–phonon couplings are strictly local. For the tetracene, the situation is quite different as the diabatic states couple highly non-locally to the vibrational degrees of freedom. The arising long-range couplings span the whole system such that combined with the binary encoding of the phononic Hilbert spaces, the calculations with n_ph = 63 would require a total system size of 461 orbitals, which are coupled over the whole range of lattice sites. It has been demonstrated⁵⁵ that already for groundstate calculations, the convergence of DMRG algorithms can be quite challenging for such a large number of lattice sites in the PS representation. In an out-of-equilibrium setup, we expect these problems to become even more severe.

The local basis optimization (LBO) is another very successful method based on rotating the local phononic Hilbert space representations into an optimized basis, allowing for a truncation of the required large local dimensions.^60,61 This effective compression generates a significant speedup in the tensor contractions and has been applied to study the dynamics of coupled electron–phonon systems following a global quench.^62–63 However, in contrast to groundstate calculations, the optimization of the local basis has to be performed at each time step during the time evolution to maintain the optimal description of the vibrational degrees of freedom.⁶⁴ As a consequence, only parts of the effective Hamiltonian can be represented in the optimized basis, the remainder still requires the full local Hilbert space dimension d = n_ph + 1 of the original phononic Hilbert space. To be more precise, the most expensive numerical operation is the application of the effective Hamiltonian, which scales as d³ in the local dimension. Within the LBO scheme, this translates into a scaling d²d_eff, where d_eff is the dimension of the optimized local basis representation. Thus, the numerical costs are reduced at the most by a factor of d_eff/d. Comparing this estimation for the speedup of the LBO method with the speedup presented in Fig. 3(c), we find a significantly larger numerical efficiency exploiting the PP mapping.

Very recently,³⁵ a TTNS representation has been suggested to capture the non-local interactions between molecular electronic states and vibrational modes. As we have shown in Sec. III B, a large number of excitations per vibrational mode inevitably need to be taken into account. Exploring a TTNS representation, the idea is to reduce the numerical costs by grouping strongly correlated regions of the wave function, reducing the required bond dimension m. Assuming a two-site TDVP update (of a rank-3 tensor), the numerically most costly operations scale as $Om3d2w2$ or $Om2d3w3$⁠, depending on the magnitude of the local dimension d and the bond dimension of the Hamiltonian representation w. Using the standard tensor-network representations, we found values of $w∼O100$ when constructing the Hamiltonian,⁶⁵ i.e., a bond dimension $m∼O100$ is an upper limit to maintain the numerical feasibility of the time evolution, which was also used in Ref. 35. However, in a tensor-network representation, the value of m is directly linked to the amount of correlations that can be captured by the corresponding bond. To resolve this conflict, entanglement renormalization by means of machine learning based optimization strategies has been applied to deduce a grouping of vibrational modes, allowing the required reduction in m. In contrast, exploiting the artificial U(1) symmetry introduced by the PP mapping, we find $w∼O10$⁠, enabling us to use m = 6144 without severe runtime losses while still working with a MPS representation, which 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.

In order to verify the applicability of our model and methods, we evaluate the absorption spectrum,⁶⁶ 

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)〉=(|LE1〉+|LE2〉)/2$⁠. An analogous analysis for the localized and dark initial state can be found in the supplementary material.

The resulting calculation is displayed in Fig. 5(a). In order to resemble the relative magnitude of the 0–0 and the 0–1 peaks of the experimental spectrum,⁹ the relaxation time is chosen to be τ = 230 fs. We find the theoretical data capturing the main experimental features for the vibrationally 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 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 Fig. 5(b), 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⁻¹. 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 Sec. 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 Eq. (1). The obtained time dependencies are displayed in Fig. 5(c). 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 Fig. 5(d), 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 and 1620 cm⁻¹, respectively. We identify these signals with vibrational modes as follows:

Note that all of these frequencies correspond to the strong peaks in the spectral densities in Figs. 2(a) and 2(d). The observation that modes 184–186 are most relevant for the energy transfer can be related to the nature of these vibrational modes, i.e., they correspond to collective vibrations of the whole molecule’s backbone (see the supplementary material). We note that previously^11,67 studies of the vibrational contributions to the excitonic dynamics found a similar relation for vibrational modes in the regime of 1420 cm⁻¹. Here, the importance of only few vibrational modes was determined either by successively eliminating oscillation frequency windows¹¹ or by combined experimental and theoretical analyses on a related molecule.⁶⁷ 

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.⁶⁸ For instance, Alvertis et al.⁶⁹ 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 Fig. 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 Fig. 6(a). In addition, we note that the triplet yield obtained from the system parameters (Table I) matches the experimental findings,⁹ i.e., we find a relatively small value of 14%.

Previous studies^11,12,70,71 found a finite time scale for the coherent dynamics. In order to determine this electronic coherence time, in Fig. 6(b), we show the time-evolution of the TT population for selected limiting cases. We find a consistent change in the dynamics at t₀ ∼ 35 fs where the slope of the TT production rate changes independently of the specific choice of the on-site potentials. As discussed in the supplementary material, this can be addressed to a change in the nature of the excitonic dynamics. At times t < t₀, all excitonic states realize a coherently driven dynamics. However, at times t > t₀, 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.¹¹ 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⁻¹ 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 the supplementary material, we investigated the renormalization of the excitonic couplings after rewriting the Hamiltonian using a Lang–Firsov transformation.⁷² 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₀ = (34.725 ± 0.165). It is a remarkable observation that the description of the phonon-assisted hopping in terms of polaron dynamics yields a perfectly matching prediction of 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,⁶¹ thereby increasing 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 Fig. 7.

We studied the TT production rate via iSF for a covalently linked tetracene dimer using a large-scale full quantum-mechanical treatment of the relaxation dynamics after photo-excitation. Our study extends the theoretical analysis conducted so far¹¹ 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₀ = 35 fs. We note that the observed time scale for coherent dynamics is in agreement with previous studies, both theoretical and experimental.^9–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₀. 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.⁹ 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,⁹ 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 renormalization techniques³⁵ could open the path to study the full quantum dynamics of generic, large organic semi-conductors.

See the supplementary material for details on the derivation of the model, as well as details on the numerics during the simulation. Furthermore, it contains a detailed analysis of the electronic RDM in combination with the derivation of the effective renormalized couplings.

We thank Dr. Xiaoyu Xie, Dr. Thomas Köhler, and Professor Uwe Manthe for helpful discussions. We acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy-426 EXC-2111-390814868 (U.S.) and the National Natural Science Foundation of China (Grant No. 22073045) (Y.X., X.Y., and H.M.). Furthermore, we acknowledge support from the Munich Center for Quantum Science and Technology (S.M., M.G., U.S., and S.P.). All calculations made use of the SyTen toolkit,⁵⁴ and some comparison benchmarks were performed using the SymMPS toolkit.⁷³ 

The authors have no conflicts to disclose.

S.M. developed the idea, wrote the code, and performed the calculations. Y.X., X.Y., and H.M. guided with chemical intuition to the chosen molecule. Y.X. and X.Y. did electronic structure calculations and the computation of the vibrational modes. S.M., Y.X., and S.P. interpreted the results and wrote the paper together. S.P. supervised the project and is one of the original authors of the main method. H.M., M.G., and U.S. discussed the results and helped interpret the data.

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