Multiscale simulations of singlet and triplet exciton dynamics in energetically disordered molecular systems based on many-body Green's functions theory

We present a multiscale model based on many-body Green’s functions theory in the GW approximation and the Bethe–Salpeter equation (GW-BSE) for the simulation of singlet and triplet exciton transport in molecular materials. Dynamics of coupled electron–hole pairs are modeled as a sequence of incoherent tunneling and decay events in a disordered morphology obtained at room temperature from molecular dynamics. The ingredients of the rates associated to the events, i.e. reorganization energies, site energies, lifetimes, and coupling elements, are determined from a combination of GW-BSE and classical polarizable force field techniques. Kinetic Monte Carlo simulations are then employed to evaluate dynamical properties such as the excitonic diffusion tensor and diffusion lengths. Using DCV5T-Me(3,3), a crystalline organic semiconductor, we demonstrate how this multiscale approach provides insight into the fundamental factors driving the transport processes. Comparing the results obtained via different calculation models, we investigate in particular the effects of charge-transfer mediated high exciton coupling and the influence of internal site energy disorder due to conformational variations. We show that a small number of high coupling elements indicative of delocalized exciton states does not impact the overall dynamics perceptively. Molecules with energies in the tail of the excitonic density of states dominate singlet decay, independent of the level of disorder taken into account in the simulation. Overall, our approach yields singlet diffusion lengths on the order of 10 nm as expected for energetically disordered molecular materials.


Introduction
Processes involving electronic excitations play a vital role in many functional molecular systems, both natural [1][2][3] and synthetic [4][5][6]. Such systems are in general characterized by a significant amount of disorder which can be either static, e.g. due to impurities or growth defects, or dynamic when resulting from thermal fluctuations. As a consequence excited state dynamics (like charge transport, exciton diffusion and conversion processes) are sensitive to the interplay between electronic structure of molecular building blocks and morphological details on nano-and mesoscale. For instance for technological applications in organic electronic devices, like organic light-emitting diodes (OLEDs) or organic solar cells (OSC), an in-depth understanding about the fundamental aspects of this interplay and how it gives rise to dynamical properties and eventually device performance is essential [7][8][9].
Computational models of these dynamics aiming at predictions of material characteristics need to take the multiscale nature of the problem explicitly into account. In general, this relies on the choice of an appropriate transport model suitable for the nature of the material and under what conditions it is studied. This choice is inherently linked to the localization behavior of the excitations [10]. For instance, transport in highly-ordered molecular crystals at low temperatures and without defects typically exhibits strong coupling between the molecules leading to delocalized electronic states and a band-like transport that can be described by the Drude model [11] or its extensions including effects of local electron-phonon coupling [12,13]. Most of OLED and OSC applications however operate at ambient conditions, in which thermal disorder is limiting the diffusion of electrons and a semiclassical model with interacting electronic and nuclear degrees of freedom can be employed [14,15]. For disordered systems at ambient conditions, when the electronic coupling is weak, dynamics of charges is typically modeled as a series of incoherent tunneling processes between localized states described, for instance, in the high-temperature limit of non-adiabatic transfer by the well-known Marcus rate [16]  where ÿis the reduced Planck constant, T the temperature, λ ij the reorganization energy, k B the Boltzmann's constant, ΔE ij the free energy difference between initial and final states, and J ij the electronic coupling element. Multiscale models employ this rate-based approach and determine the ingredients to equation (1) from a combination of quantum (QM) and classical (MM) techniques for an explicitly simulated disordered molecular morphology.
Microscopic simulations of charge transport linking electronic and classical degrees of freedom within the regime of incoherent tunneling have seen great advances in recent years [8,17,18], providing insight into, e.g. the role of molecular polarizabilities in the design of binary emission layers in OLEDs [19] or the importance of mesoscale effects on the level alignment and band bending across a donor-acceptor interface in OSC [20]. Compared to the research into separate electron or hole transport, similar models aiming at investigating the dynamics of excitons, or coupled electron-hole pairs, are relatively sparse [21][22][23][24][25][26][27]. Most studies have focused on singlet transport, as electronic couplings can be approximated by dipole-dipole coupling [28], and comparatively small system sizes with little or no disorder. One of the reasons for this is the large computational cost typically associated to the quantum-level description of excitons with good accuracy, which is required for the evaluation of reorganization and site energies, and coupling elements. The use of many-body Green's functions theory in the GW approximation with the Bethe-Salpeter equation [29,30] (GW-BSE) has recently gained increasing attention in these contexts [31][32][33][34]. Hybrid GW-BSE/MM methods have been developed which are capable of predicting excitation energies of chromophores in complex molecular environments at good accuracy [35,36]. It has also been shown how singlet and triplet exciton coupling elements including direct and charge-transfer (CT) mediated mechanisms can be efficiently extracted on equal footing from this approach [37].
In this paper, we extend the GW-BSE description of excited states in complex molecular systems from a purely static picture to a dynamic one in the multiscale framework described above. As a prototypical system, we study the singlet and triplet exciton dynamics in a DCV5T-Me (3,3) crystal [38,39]. Taking an initial structure from experimental crystal orientations, the system is equilibrated at room-temperature using classical molecular dynamics (MD). Then, we determine all ingredients to the Marcus rate from GW-BSE calculations on monomers and dimers, and a combination with classical methods in mixed QM/MM setups. In this way, we take the long-range coupling of singlet excitons, as well as the CT mediated short-range transfer into account on equal footing, analyze the effects of large couplings and thermal fluctuations [40] on the overall dynamics, and scrutinize the contributions of different types of energetic disorder on observed exciton diffusion properties. This paper is organized as follows: Insection 2, we briefly recapitulate steps in the multiscale model, including computational details of our MD calculations, the essentials of GW-BSE, and the rate-based dynamic models. Results obtained for the reorganization energies, site energies, coupling elements, and singlet and triplet diffusion, as well as their analysis and discussion are given insection 3. A brief summary concludes the paper.

Methodology
The basic workflow of the rate-based model of exciton dynamics is depicted infigure 1. Representative morphologies at room temperature are generated from classical MD simulations. Each generated snapshot is then partitioned into localization sites for the excitons, called segments. Based on inter-segment distances, a neighbor list of hopping pairs is created, for which transfer rates need to be determined. A combination of quantum-chemical calculations on the level of GW-BSE with classical atomistic electrostatics is employed to evaluate the rate ingredients (reorganization energy, site energy, and coupling elements) for all pairs. Diffusion tensors and diffusion length are finally obtained from simulating exciton dynamics explicitly by means of a kinetic Monte Carlo (kMC) algorithm and analyzing its trajectories. In the following, we briefly summarize the details of the individual steps in this procedure. Figure 2 shows different views of the crystal structure of DCV5T-Me(3,3) as determined from x-ray data [38]. While the published cell is triclinic, this orthogonal unit cell with side lengths as indicated in panels (a) and (b) contains eight molecules in staggered arrangements in all three directions. Molecules within a π-stack (along the z-direction) are displaced about 8 Å along the backbone (y-direction) with respect to each other. As a consequence of this non-ideal π-stacking, the minimal center-of-mass distance along the stack is 8.8 Å, much larger that the z-distance of 3.7 Å. Based on this structure, an orthorhombic (7×2×14) supercell with side lengths = L 11.06 nm containing 1568 molecules was created as initial structure. To introduce disorder at ambient conditions, MD simulations of 1.5 ns were run in the NpT ensemble at p=1 bar and T=300 K using the Berendsen barostat [41] and stochastic velocity rescaling thermostat [42]. All simulations were performed using an OPLS-based force field with specific parameterization for the class of terminally dicyanovinyl-substituted oligothiophenes [43], and the GROMACS software [44]. After equilibration, the resulting supercell is slightly contracted in x-direction (L x =10.74 nm), while the y-and zdirection remain practically unchanged (L y =10.78 nm and L z =10.39 nm). Six frames in 20 ps intervals starting from time 1.32 ns were taken for the evaluation of the exciton dynamics.

Electronic excitations viaGW-BSE
Electronically excited states for monomers and dimers extracted from those snapshots are calculated within the GW-BSE approach. In comparison to time-dependent density-functional theory (TD-DFT) GW-BSE does not require explicitly tuned functionals to accurately model charge transfer excitations, while retaining the favorable ( )  n 4 scaling with system size. Calculations were performed with the VOTCA-XTP software [36] as briefly summarized below. More exhaustive reviews can be found in [45,46].
Since single-and two-particle excitations are constructed upon a DFT ground state, first the Kohn-Sham (KS) energy levels and wave functions need to be obtained by solving Here, V ext is the external potential (either of bare nuclei or pseudo atoms), V H the Hartree potential, and V xc the exchange-correlation potential.  Employing the GW approximation of many-body Green's functions theory, as introduced by Hedin [29], quasi-particle (QP) states representing independent electron and hole excitations are found by solving the QP equations: Central difference toequation (2) is the occurrence of the energy-dependent self-energy operator ( ) S ¢ E r r , , in place of the exchange-correlation potential. The self-energy operator is evaluated using the one-body Green's function in QP approximation where W is the dynamically screened Coulomb interaction. To determine the latter, first the random-phase approximation (RPA) [30,47,48] is used in calculating the polarization P, and with it the microscopic dielectric function follows from convoluting the polarization with the bare Coulomb interaction v, i.e. ò=1−vP. This, in turn, is inverted and another convolution with the bare Coulomb interaction finally yields W=ò −1 v. Instead of calculating the full frequency dependence explicitly, a generalized plasmon-pole model [49] is used which extends the RPA result for ω=0 (static polarization) and the associated static dielectric function to the dynamic one. To avoid deviations due to the use of the ground state KS wave functions and energies in practical calculations to determine both G and W, we use a double self-consistent evGW scheme in which quasiparticle energies E n QP are updated both in the RPA evaluation of the polarization and in the self-energy ( ) S E n QP . Within the Tamm-Dancoff approximation [46] 3 (TDA) coupled electron-hole excitations, generated e.g. after optical excitation, are expressed as a linear combination of products of occupied (v) and empty (c) QP functions .  , and the bare exchange (K x ) and screened direct (K d ) terms of the electron-hole interaction kernel, respectively. Finally, Ω is the transition energy of the optical excitation, and η=2 for singlet and η=0 for triplet transitions.
For the practical GW-BSE calculations in this paper, single-point KS calculations are performed using the Gaussian09 package [52], the PBE functional, Stuttgart/Dresden effective core potentials [53], and the associated basis sets that are augmented by additional polarization functions [54] of d symmetry. All steps of the actual GW-BSE calculations employ the Gaussian-type orbital based implementation in the VOTCA-XTP package [17,32,36,55]. It uses an auxiliary basis set to express the quantities occurring in the GW self-energy operator and the electron-hole interaction in the BSE using resolution of identity techniques. In [36], this combination of basis sets and functional has been shown to be an efficient minimal basis, yielding excitation energies in close agreement to all-electron calculations based on cc-pVTZ basis at significantly reduced computational cost. Further technical details can be found in [36].
A rate-based model of exciton dynamics requires as a first step the definition of all possible hopping pairs in the system, for which the rates need to be determined. Such a pair list has to account for the different transfer mechanisms of singlet and triplet excitons in a two-region setup. At first a neighbor search is performed, in which all molecules with a closest contact distance of 0.6 nm or less are included. For these pairs both the singlet and triplet excitonic coupling is mediated by direct interactions of the electron-hole wave functions, which decay exponentially with distance analogously to electronic coupling in charge transport [56]. As we have shown in [37] (see, i.e. figure 7(b) therein), this choice ensures that beyond this cutoff, the singlet couplings can be determined in perfect agreement with explicit GW-BSE results using a Coulomb coupling with a classical 3 For typical donor molecules used in OSCs, we showed that the use of the TDA overestimates π-π transition energies by 0.2 eV but yields correct character of the excitations [32].
representation of the transition densities by atomistic point charges (TrESP). For the singlet model, all additional molecules within distance of 3.5 nm to each other are added to the pair list, and their coupling will be evaluated with TrESP.
For each of the in total resulting roughly 372 000 (6000) pairs for singlet (triplet) exciton dynamics, the transfer rate needs to be calculated according toequation (1) from the reorganization energies for excitation l g a 0 and de-excitation l g a 0 , and the site energies of single monomers ΔE a , as well as the excitonic coupling in dimers J ab . Within the Marcus picture, these monomer energies are determined as total energy differences according to represents the total energy of molecule a in either ground (i=0) or excited (i=γ) state at ground (I=0) or excited (I=Γ) state geometry. Note that the above is inherently linked to the Marcus picture and the assumption of a single classically treated promoting mode. In cases with a strong coupling of the molecular vibrations with the electronic states, alternative rate expressions to the Marcus rate should be considered, in which the vibrational modes are also treated quantum-mechanically [17,18,57] In a molecular material, the energies in equations (8)-(10) consists of (QM) contributions of the single molecule and the effects of interactions with the full molecular environment. We evaluate the latter as a correction to the QM energies using a polarizable molecular mechanics (MM) approach, yielding QM+MM energies as ( ) . The MM model uses atomic partial charge representations and Thole's atomic polarizabilities, which need to be parameterized based on QM data, such as the molecules' electrostatic potential and polarizabilities. In principle, this approach requires for frozen soft molecular conformations (given by, e.g. torsional conformers of thiophene units after post-processing to remove the fast degrees of freedom representing the promoting mode in Marcus theory) the constrained optimization of ground and excited state geometries, and individual partial charge and polarizability fits. Especially due to geometry optimizations and molecular polarizability calculations for the excited states in GW-BSE relying on numerical gradients, this is computationally prohibitively expensive. Instead, we assume that since the DCV5T-Me(3,3) molecules in the crystal only fluctuate marginally around the ideal conformation (i.e. the lack of cis-trans rearrangements) the geometric relaxation effects are similar for all conformers. This allows to approximate the molecule-specific reorganization energies inequation (10) by molecule-independent ones, l l » g g a 0 0 and l l » g g a 0 0 . For the internal part of the site energy, it follows that the adiabatic excitation energy can be approximated as i.e. as the vertical excitation energy (absorption) reduced by the generic reorganization energy for excitation. When forming the site energy difference ΔE ab inequation (1), the constant λ 0 γ drops out and the associated disorder is given by variations of the absorption energy alone. Within these approximations, we first optimize the geometries of DCV5T-Me (3,3) in the ground state using DFT as well as in the first singlet and triplet excited states using GW-BSE, respectively. At the optimized geometries, CHELPG [58] derived partial charges {q i } are calculated and the molecular polarizability tensor a QM is determined. Then, atomic polarizabilities are adjusted so that the associated a MM of the Thole model reproduces the QM reference.
In a next step, a vacuum GW-BSE calculation is performed for each molecule extracted from the MD snapshots and the first singlet (γ=S) and triplet (γ=T) excitation energies Ω S , Ω T and electron-hole wave functions are stored. From the parameterized polarizabilities and partial charges, the MM contributions to the site energies ΔE a are calculated using a cutoff based embedding scheme [17], in which within a radius of 3 nm around molecule a polarized interactions are taken into account while for static interactions are active up to 6 nm.
Excitonic couplings for singlets and triplets among the short-range pairs inside the 0.6 nm cutoff mentioned above are calculated using the dimer projection technique GW -BSE-DIPRO as introduced in [37], in which the excitonic coupling is formed as the expectation value of a supramolecular BSE Hamiltonian with electron-hole wave functions for excitations localized on two separated chromophores. Within this approach, accounting for the effects of coupling mediated by intermolecular CT excitations is possible via perturbation theory. As the singlet coupling beyond this first neighbor shell is dominated by Förster like transition density coupling, we represent the molecular GW-BSE transition densities by CHELPG [58] point charges. Coulomb interactions among these transition electrostatic potential (TrESP) charges [59,60] are then used to approximate the longrange singlet coupling classically.
From the electronic couplings, the site energies and the reorganization energies, the Marcus rates equation (1) for exciton transfer are calculated, and simulations of exciton dynamics is performed using kMC. From its trajectory it is possible to extract, e.g. the exciton's diffusion tensor. This requires that the carrier is in a steady-state and has reached the diffusive regime of transport, which can be assumed to be the case for long-lived triplet excitons. Singlets, in contrast, are short-lived since they can decay radiatively back to the ground state. The inverse radiative decay time τ rad , for a singlet state, is given by Einstein's formula for spontaneous emission [61]: where ε is the dielectric constant of the material, α the fine structure constant, Ω S the singlet excitation energy and m Tr the transition dipole moment of the transition. Neglecting non-radiative decay, due to the difficulties of calculating electron-phonon interaction, the radiative decay rate serves as a lower bound for the full decay rate: Viaequation (14), we explicitly include singlet decay in the KMC simulations, and record the distance excitons travel before decaying.

Polarizabilities and reorganization energy
In the following we focus on the properties of the first singlet and triplet excited states of DCV5T-Me(3,3), respectively. These lowest energy excitations of π→π * character are well separated from the rest of the excitation spectrum and in case of the singlet are the optically active transition. After geometry optimizations as outlined insection 2, molecular polarizability tensors for ground and excited states, as well as the reorganization energies are calculated. After diagonalization, the eigen frame of the polarizability tensor is aligned with the molecular axes as indicated infigure 2: the long dimension (backbone) of the DCV5T-Me(3,3) molecule is aligned in y-direction of the coordinate system, while the short dimension of its molecular plane is aligned to the x-direction. The z-direction is normal to the molecular plane. In this frame, the diagonal of the ground state polarizability tensor is given by and T, respectively. This observation is consistent with the nature of the dominant π→π * transition in these excited states. Promotion of an electron to an anti-binding orbital leaves it in a unbound, extended state, which responds more strongly to external fields. Due to the repulsive exchange interaction in the BSE Hamiltonian for singlet excitations (see equation (7)), this effect is stronger for S than for T.
The reorganization energies were evaluated based on the optimized structures using equation equation (10) as sum of term for excitation and de-excitation and result as λ S =0.37 eV for singlets and λ T =0.77 eV for triplets, respectively. Albeit at different levels of theory, similarly high reorganization energies have been reported for other organic compounds, i.e. λ S =0.35 eV for the singlet state of napthalene (SCS-CC2@cc-pVDZ) [62] and λ T =0.75 eV for the triplet state of Alq 3 (ADC(2)@SVP) [24].

Site energies
The site energy distributions obtained from the QM+MM approach outlined insection 2 are approximately Gaussian for both singlets and triplets. Table 1 shows the standard deviations σ of internal quantum (σ QM ) and external classical (σ MM ) contributions, as well as of the total site energies (σ QM+MM ) for singlet and triplet excitons, respectively. The distributions are largely identical across the six individual frames. We find that s s s » + indicating that the two contributions to the site energies are uncorrelated. For the calculation of singlet decay rates viaequation (14), the ground to excited state transition dipole moments m Tr were obtained from individual vacuum monomer GW-BSE calculations. To approximately include the environment effects Ω S is replaced by ΔE QM+MM . As to our knowledge the exact dielectric constant of DCV5T-Me(3,3) has not been measured, we assume ò≈3, which is a typical value for organic materials [63]. Figure 3 shows the distribution of the obtained singlet exciton lifetimes. All six production frames exhibit similar behavior, with the distributions peaking at about τ rad =0.5 ns, which is comparable to reported singlet exciton lifetimes for similar compounds [64].

Coupling elements
Triplet coupling elements calculated using the GW -BSE-DIPRO method arise from the term K d of the electronhole kernel in the BSE Hamiltonian(equation (7)). This term involves the short-ranged direct interaction of monomer wave functions, similar to the coupling of individual electron and hole single particle states. The coupling is very sensitive to molecular arrangements and in particular exhibits an exponential distancedependence. Figure 4 shows the distributions of    To elucidate the nature of the two peak structure in the distributions, we consider infigure 5(a) the same data as infigure 4 now resolved with respect to the center-of-mass distance Δr of the molecules in the pair 4 . The red regions marked by (1), (2), and (3) indicate the highest concentrations of triplet couplings at Δr of approximately 0.9 nm, 1.4 nm, and 2.7 nm, respectively. These distances correspond to neighbors in different unique directions in the crystal structure, as shown infigure 5(c). Structure (1) corresponds to the shifted πstack arrangement within a z-column as already mentioned in section 2.1. Structure (2) is a staggered configuration between two molecules from adjacent columns, e.g. as the bottom two molecules infigure 2(a). Both structures contribute to the highest observed triplet couplings in the right peak offigure 4. The left peak can be associated to tail-to-tail orientations as in structure (3), which only exhibit weak coupling. Figure 6(a) shows the distribution of singlet couplings in all six frames, which appear nearly identical for J 2 <10 −2 eV 2 . For larger singlet coupling strengths we note (on the logarithmic scale) small variations. To assess the relative importance of the short-ranged GW -BSE-DIPRO and the long-ranged TrESP couplings, we resolve the individual contributions for frame 1 infigure 6(b). It is apparent that the singlet transfer integrals evaluated explicitly on GW-BSE level contribute almost exclusively to the region with J 2 >10 −2 eV 2 , also as a consequence of the splitting of the neighbor list. GW -BSE-DIPRO and TrESP are of similar importance only in a small interval 10 −3 eV 2 <J 2 <10 −2 eV 2 . Fromfigure 6(c) (note the linear y-scale) the limited region in which the quantum-mechanical coupling elements contribute is again clearly visible. For lower values of coupling   4 Note that we consider center-of-mass distances up to 2.8 nm. This is not in contradiction with the cutoff distance for the neighbor list of 0.6 nm for the closest-contact distance. elements the transition density coupling completely dominates the total distribution. The analysis of the distance dependence infigure 5(b) shows the expected slow decay due to the Coulomb coupling of the transition densities. The influence of thermal fluctuations can also be clearly seen from the smeared out appearance of the distributions.
Overall, we observe a median coupling of á ñ = J 0.044 eV, which is significantly smaller than the reorganization energy λ S =0.37 eV. However, we also find a few (e.g. 84 in frame 1) singlet couplings are larger than λ S , some of them exceeding even | | > J 1 eV. Even though only 1.39% (0.02%) of the GW -BSE-DIPRO (total) transfer integrals are unusually high, closer inspection of their nature is warranted. Next to the mechanisms of direct exchange and transition density interactions, exciton states can also couple via intermediate bi-molecular CT states in the dimer [37], if the latter lie energetically close to the monomerlocalized states. This mechanism is accounted for in the GW -BSE-DIPRO method by expanding effective monomer states to second order in actual monomer states and subsequent orthogonalization. However, this perturbative approach overestimates CT-mediated couplings if the CT energies are close to resonance with the ones of a monomer. The obtained high coupling elements indicate that the assumption of monomer-localized excitons is not always justified. Instead in some configurations the exciton can be expected to delocalize over the dimer.

KMC simulations
For the simulation of the exciton dynamics within the rate-based model we consider two different cases: in the first only the external contributions to site energies (leading to σ MM in table 1) are taken into account, while the second additionally contains the internal contributions (corresponding to σ QM+MM ). In all cases, the initial site for the KMC trajectory is chosen randomly and independently from each other and the environment. We first consider triplet diffusion and determine the elements of the diffusion tensor D αβ (in the Cartesian directions, i.e. α, β=x, y, z) from long KMC trajectories with simulation time t according to where DR denotes the distance vector between start and end point, and · á ñ an ensemble average. For each frame, 20 KMC simulations were run, with on average 9×10 11 steps.
The diffusion tensor without internal disorder exhibits strong anisotropy with diffusion along the y-axis (D yy =54.26 ×10 −9 cm 2 s −1 ) one order of magnitude larger than along the z (D zz =5.27 ×10 −9 cm 2 s −1 ) and two orders of magnitude larger than diffusion along the x (D xx =0.84 ×10 −9 cm 2 s −1 ) axis. As the energetic disorder is comparatively small (σ MM =44 meV) the triplet diffusion is governed mainly by the excitonic couplings. As we discussed in the previous section (see figure 5(a) and (c)) the strongest coupling between molecules is present in the y−z plane with z being the stacking direction. Due to shift of neighboring molecules along the y-axis, the displacement for each jump is about three times larger along the y-axis than the z-axis, leading to observed anisotropy.
Taking the internal contributions to the site energies into account, reduces the magnitude of the elements of the diffusion tensor by roughly three orders of magnitude. At the same time, the anisoptropy of the diagonal elements is reduced to about 1:10:5. Both observations are a consequence of the increased energetic disorder of σ QM+MM =126 meV, which becomes relatively more important than the more topological connectivity given by the coupling elements.
An effective triplet exciton diffusion length ℓ | | T can be extracted from the trace of the diffusion tensor according to where τ T is the non-radiative triplet lifetime. Since we cannot explicitly calculate it, we make use of the fact that triplet lifetimes are typically six orders of magnitudes higher than the respective singlet lifetimes [64] and use τ T =1 × 10 −3 s in the following. Figure 7(a) shows the obtained diffusion lengths (also resolved in directions according to ℓ t = D i i i T ) for the simulations without (blue) and with internal energetic disorder (red). The results in general reflect both the anisoptropy of diffusion in the crystalline system as well as the significant influence of the internal site energy disorder, reducing the overall triplet diffusion length from 44.9±0.3 nm to 1.45±0.01 nm.
For the estimation of the diffusion length of singlet excitons ℓ | | S we take a different approach without reference to the diffusion tensor. Instead we explicitly take the exciton decay via the inverse radiative lifetimeequation (14) into account in the kMC simulations.
We first again focus on the case with no internal disorder and scrutinize the influence of the small number of unphysically high coupling elements resulting from CT resonances. For a single frame, we consider two limiting scenarios: one in which the high coupling elements are taken at face value, and one in which those with | | l > J are set to zero. KMC simulations with 10 4 singlet insertions are performed for both cases. The full approach yields an average diffusion length of ℓ | | =  62.7 0.6 nm S , while we obtain ℓ | | =  61.9 0.6 nm S with the high couplings removed. The differences between both methods being only 1.2% clearly indicates that the strongly coupled singlet states on a fraction of the total number of hopping pairs has negligible influence on the overall singlet dynamics in DCV5T-Me (3,3).
Based on this we proceed with the full model and run KMC simulations for all six snapshots, now with 2×10 6 singlet insertions each. The direction-resolved and total averages are shown infigure 7(b). Without taking internal disorder into account (blue bars), we obtain an average diffusion length of 55.3±3.1 nm. Due to the much more connected graph underlying the hopping-based dynamics of singlet excitons due to the contribution of the long-range TrESP coupling elements, the anisotropy in the diffusion is with only 1:2:1 along the three main axes much lower than in the triplet case. For the same reasons, singlet diffusion is also much less affected by including internal disorder in the site energies, as can be seen by the red bars infigure 7(b). The average diffusion length is reduced to 16.3±4.0 nm, about 30% of the value without internal disorder (as compared to 3% for triplets), while the anisotropy remains practically unchanged.
From the KMC simulations it is possible to analyze the individual singlet exciton diffusion processes in more detail. In particular we are interested in the dynamic decay process of the singlets and how it is affected by the molecular properties of DCV5T-Me (3,3). To this end we record for each molecule i the number of times N d i a singlet exciton decays there and determine from that an effective decay probability where N d total is the total number of decayed excitons. We then rank the molecules in the simulation cell according to P i d, eff and focus on the 100 molecules with the highest effective decay probability. Infigure 8 we show this decay probability (blue line) as a function of molecule rank for frame 1. For both models without and with internal disorder, we find that most excitons decay on a very small number of sites. To elucidate the reasons for this inhomogeneous decay, we also consider both the relative radiative lifetimet t t = i i rad, rad, rad max per molecule (yellow line) and the site energies ΔE i (green line). The data indicates thatt rad is relatively uniform across the molecules and actually exhibits a small increase, i.e. singlet excitons live longer on the sites they are most likely to decay at. As can also be seen, this general observation is unaffected by the internal disorder model considered in the calculations. In both cases, the fairly homogeneoust rad cannot give rise to the observed decay characteristic. In contrast, the comparison of the P i d, eff with the site energy reveals that the decay predominantly occurs on molecules with energies in the lower tail of the site energy distribution. When singlet excitons transfer to such sites, their total escape ratew w = å i j ij is significantly reduced. While the decay rate is roughly 2× 10 9 s −1 , we find that for the 100 molecules shown infigure 8,w is typically around 5×10 10 s −1 in the model without internal disorder. Only for the molecules with the ten highest effective decay probabilities and lowest site energies, the total escape rate is on the same order of magnitude as trad 1 , and as a consequence the probability of a decay event noticeably increases. In the model with internal disorder taken into account, this behavior is even more pronounced. For the lowest energy molecule which practically dominates the singlet decay, the total escape rate is only 2× 10 8 s −1 , about one order of magnitude lower than the inverse radiative lifetime. This molecule effectively has the characteristics of a deep trap for the excitons.

Summary
To summarize, we have presented a multiscale workflow based on many-body Green's function theory in the GW approximation and the Bethe-Salpeter equation (GW-BSE) to simulate singlet and triplet exciton dynamics in molecular materials as a series of incoherent transfer processes. Using DCV5T-Me(3,3), a crystalline organic semiconductor, we demonstrated how a combination of quantum-mechanical and classical methods can be used to determine the ingredients of excitonic transfer and decay rates, i.e. reorganization energies, site energies, lifetimes, and coupling elements. KMC simulations were then employed to evaluate dynamical properties such as the excitonic diffusion tensor and diffusion lengths. Comparing the results obtained with different calculation models, we showed that a small percentage of high coupling elements indicative of delocalized exciton states does not impact the overall dynamics perceptively. Overall, our approach provides insight into the mechanism of transport including the importance of energy traps for singlet decay and yields singlet diffusion lengths on the order of 10 nm as expected for disordered molecular materials. Results also show that for quantitative predictions a careful treatment of the internal site energy disorder is essential.