Multiscale Simulations of Singlet and Triplet Exciton Dynamics in Disordered Molecular Systems based on Many-Body Green’s Functions Theory

We present a multiscale model based on Many-Body Green’s Function Theory in the GW approximation and the Bethe-Salpeter Equation (GW-BSE) for the simulation singlet and triplet exciton transport in molecular materials. Dynamics of coupled electron-hole pairs is 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 were 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 with 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 disordered molecular materials. We present a multiscale model based on Many-Body Green’s Function Theory in the GW approximation and the Bethe-Salpeter Equation ( GW -BSE) for the simulation singlet and triplet exciton transport in molecular materials. Dynamics of coupled electron-hole pairs is 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 ﬁeld techniques. Kinetic Monte Carlo simulations were then employed to evaluate dynamical properties such as the excitonic diﬀusion tensor and diﬀusion 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 with diﬀerent calculation models, we investigate in particular the eﬀects of charge-transfer mediated high exciton coupling and the inﬂuence 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 diﬀusion lengths on the order of 10 nm as expected for disordered molecular materials.


I. INTRODUCTION
Processes involving electronic excitations play a vital role in many functional molecular systems, both natural 1-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 this context, microscopic simulations of charge transport linking electronic and classical degrees of freedom have seen great advances in recent years, 8,10,11 providing insight into, e.g., the role of molecular polarizabilities in the design of binary emission layers in OLEDs 12 or the importance of mesoscale effects on the level alignment and band bending across a donor-accepor interface in OSC. 13 For disordered systems at ambient conditions, dynamics of charges is typically modeled as a series of incoherent tunneling processes described, for instance, in the high-temperature limit of non-adiabatic transfer by the well-known Marcus rate 14 whereh 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 eq. 1 from a combination of quantum (QM) and classical (MM) techniques for an explicitly simulated disordered molecular morphology. Compared to the study the transport of separated electrons or holes, similar models aiming at investigating the dynamics of excitons, or coupled electron-hole pairs, are relatively sparse. [15][16][17][18] Most studies have focused on singlet transport, as electronic couplings can be approximated by dipole-dipole coupling, 19 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 20,21 (GW -BSE) has recently gained increasing attention in these contexts. [22][23][24][25] Hybrid GW -BSE/MM methods have been developed which are capable of predicting excitation energies chromophores in complex molecular environment at good accuracy. 26,27 It has also been shown how singlet and triplet exciton coupling elements including direct and charge-transfer mediated mechanisms can be efficiently extracted on equal footing from this beyond-DFT approach. 28 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 prototpyical system, we study the singlet and triplet exciton dynamics in a DCV5T-Me(3,3) crystal. 29,30 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 can be extracted 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 charge-transfer mediated short-range transfer into account on equal footing, analyze the effects of large couplings and thermal fluctuations 31 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 sec. II, 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 sec. III. A brief summary concludes the paper.

II. METHODOLOGY
The basic workflow of the rate-based model of exciton dynamics is depicted in fig. 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 quantumchemical 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. 29 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 zdirection) are displaced about 8Å along the backbone (ydirection) with respect to each other. As a consequence of this non-ideal π-stacking, the minimal center-of-mass distance along the stack 8.8Å, much larger that the zdistance of 3.7Å. Based on this structure, an orthorhombic (7×2×14) supercell with side lengths L 0 x = 11.06 nm, L 0 y = 10.80 nm, and L 0 z = 10.37 nm containing 1568 molecules was created as inital structure. To introduce disorder at ambient conditions, molecular dynamics simulations of 1.5 ns were run in the N pT ensemble at p = 1 bar and T = 300 K using the Berendsen barostat 32 and stochastic velocity rescaling thermostat. 33 All simulations were performed using an OPLS-based force field with specific parameterization for the class of terminally dicyanovinyl-substituted oligothiophenes, 34 and the GROMACS software. 35 After equilibration, the resulting supercell is slightly contracted in x-direction (L x = 10.74 nm), while the y-and z-direction 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.

B. Electronic Excitations via GW -BSE
Electronically excited states for monomers and dimers extracted from those snapshots are calculated within the GW -BSE approach using the VOTCA-XTP software 27 as briefly summarized below. More exhaustive reviews can be found in Refs. 36,37. 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, 20 quasi-particle (QP) states representing independent electron and hole excitations are found by solving of the  29 Each cell contains eight molecules in staggered arrangements. An orthorhombic (7 × 2 × 14) supercell with side lengths L 0 x = 11.06 nm, L 0 y = 10.80 nm, and L 0 z = 10.37 nm is used as initial structure for the MD simulations.
quasi-particle equations: Central difference to eq. 2 is the occurrence of the energydependent self-energy operator Σ(r, r ′ , E) in place of the exchange-correlation potential. The self-energy operator is evaluated using the the one-body Green's function in quasi-particle approximation as where W the dynamically screened Coulomb interaction. To determine the latter, first the random-phase approximation (RPA) 21,38,39 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 calculating the full frequency dependence explicitly, a generalized plasmon-pole model 40 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 Kohn-Sham 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 QP n are updated both in the RPA evaluation of the polarization and in the self-energy Σ(E QP n ). Within the Tamm-Dancoff approximation 37,41 (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) quasi-particle functions Here, A γ vc is the electron-hole amplitude of excited state γ which is obtained as solutions to the Bethe-Salpeter equation 36,42,43 The Hamiltonian comprises contributions from free interlevel transition energies, D = E QP virt − E QP occ , 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 Kohn-Sham calculations are performed using the Gaussian09 package, 44 Stuttgart/Dresden effective core potentials, 45 and the associated basis sets that are augmented by additional polarization functions 46 of d symmetry. All steps of the actual GW -BSE calculations employ the Gaussian-type orbital based implementation in the VOTCA-XTP package. 10,23,27,47 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 Ref. 27, this combination of basis sets and functional has been shown to be an efficient minimal basis, yielding excitation in close agreement to all-electron calculations based on cc-pVTZ basis at significantly reduced computational cost. Further technical details can be found in Ref. 27.

C. Rate-based model
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. 48 In a second step, all additional molecules within distance of 3.5 nm to each other are added to account for long-range Coulomb coupling of singlet excitons.
For each of the in total resulting roughly 372000 (6000) pairs for singlet (triplet) exciton dynamics, the transfer rate needs to be calculated according to eq. 1 from the reorganization energies for excitation λ 0γ a and de-excitation λ γ0 a , 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 where E i a (I) 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. In a molecular material, these (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 E i,QM+MM The MM model uses atomic partial charge representations and Thole's atomic polarizabilities, which need to be parametrized based on QM data, such as the molecules' electrostatic potential and polarizabilities. In principle, this approach requires for frozen soft molecular conformation (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 polarizabilities 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 eq. 8 by molecule-independent ones, λ 0γ a ≈ λ 0γ and λ γ0 a ≈ λ γ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 eq. 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 49 derived partial charges {q i } charges are calculated and the molecular polarizability tensor α QM is determined. Then, atomic polarizibilities are adjusted so that the associated α MM of the Thole model reproduces the QM reference.
In a next step, a GW -BSE calculation is performed for each molecule extracted from the MD snapshots and the excitation energies of the first singlet (γ = S) and triplet (γ = T ) excitation energies Ω S , Ω T and electronhole wave functions are stored. From the parametrized polarizabilities and partial charges, the MM contribution to the site energies ∆E a are calculated using a cut-off based embedding scheme, 10 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 Ref. 28, 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 49 point charges. Coulomb interactions among these Transition electrostatic potential (TrESP) charges 50,51 are then used to approximate the long-range singlet coupling classically.
From the electronic couplings J AB , the site energies and the reorganization energies, the Marcus rates ω A→B eq. 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 steadystate 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: 52 where ε is the dielectric constant of the material, α the fine structure constant, Ω S the singlet excitation energy and µ 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 eq. 10, we explicitly include singlet decay in the KMC simulations, and record the distance excitons travel before decaying.

A. Polariziabilities 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 sec. II, 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 fig. 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-directions is normal to the molecular plane. In this frame, the diagonal of the ground state polarizability tensor is given by diag(α 0 ) = α G = (75, 293, 34)Å 3 . The excited states show approximately the same polarizability in the x-and z-directions as in the ground state, with α S = (75, 1022, 27)Å 3 and α T = (70, 646, 33)Å 3 . Along the backbone, in contrast, the polarizability is about 3.5 and 2.2 times larger for S 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 (cf. eq. 7), this effect is stonger for S than for T .
The reorganization energies were evaluated based on the optimized structures using equation eq. 8 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) 53 and λ T = 0.75 eV for the triplet state of Alq 3 (ADC(2)@SVP). 16

B. Site energies
The site energy distributions obtained from the QM+MM approach outlined in sec. II are approximately Gaussian for both singlets and triplets. Table I   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 σ 2 QM+MM ≈ σ 2 QM + σ 2 MM indicating that the two contributions to the site energies are uncorrelated.
For the calculation of singlet decay rates via eq. 10, the ground to excited state transition dipole moments µ Tr were obtained from individual 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. 54 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. 55

C. Coupling elements
Triplet coupling elements calculated using the GW -BSE-DIPRO method arise from the term K d of the electron-hole kernel in the BSE Hamiltonian (eq. 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 distance-dependence. Figure 4 shows the distributions of log 10 [(J/eV) 2 ] for triplet excitons in the six frames. Overall, the distributions have similar features, with a main peak corresponding to values of J 2 ≈ 10 −5 eV 2 and a second one about five orders of magnitude lower. Visible differences among the frames are clearly effects of thermal fluctuations on the atomic positions of the molecules. The maximum triplet coupling J max T = 0.12 eV is smaller than the reorganization energy λ T = 0.77 eV. This corroborates the assumption of hopping-type triplet transport in DCV5T-Me (3,3).
To elucidate the nature of the two peak structure in the distributions, we consider in fig. 5(a) the same data as in fig. 4 now resolved with respect to the center-of-mass distance ∆r of the molecules in the pair. 56 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 fig. 5(c). Structure (1) corresponds to the shifted π-stack arrangement within a z-column as already mentioned in sec. II A. Structure (2) is a staggered configuration between two molecules from adjacent columns, e.g., as the bottom two molecules in fig. 2(a). Both structures contribute to the highest observed triplet couplings in the right peak of fig. 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 fig. 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 fig. 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 elements the transition density coupling completely dominates the total distribution. The analysis of the distance dependence in fig. 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 charge-transfer (CT) states in the dimer, 28 if the latter lie energetically close to the monomer-localized 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 momoner-localized excitons is not always justified. Instead in some configurations the exciton can be expected to delocalize over the dimer.

D. Kinetic Monte-Carlo 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 tab. I) are taken into account, while the second additional contains the internal contributions (corresponding to σ QM+MM ). 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 ∆R 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) one order of magnitude larger than along the z (D zz = 5.27 ·10 −9 cm 2 /s) and two orders of magnitude larger than diffusion along the x (D xx = 0.84 ·10 −9 cm 2 /s) 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 fig. 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 longer than the respective singlet lifetimes 55 and use τ T = 1 · 10 −3 s in the following. Figure 7(a) shows the obtained diffusion lengths (also resolved in directions according to ℓ i = √ D ii τ 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 eq. 10 into account in the kinetic Monte-Carlo 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 |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 |ℓ S | = 62.7 ± 0.6 nm, while we obtain |ℓ S | = 61.9 ± 0.6 nm 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 fig. 7 connected graph underlying the hopping-based dynamics of singlet excitons due to the contribution of the longrange 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 fig. 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 i d a singlet exciton decays there and determine from that an effective decay probability where N total d is the total number of decayed excitons. We then rank the molecules in the simulation cell according to P eff d,i and focus on the 100 molecules with the highest effective decay probability. In fig. 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 lifetimeτ rad,i = τ rad,i /τ max rad per molecule (yellow line) and the site energies ∆E i (green line). The data indicates thatτ 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 homogeneous τ rad cannot give rise to the observed decay characteristic. In contrast, the comparison of the P eff d,i 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 rateω 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 fig. 8,ω 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 τ −1 rad , 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.

IV. 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. Kinetic Monte-Carlo 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.