Transiently delocalised hybrid quantum states are gateways for efficient exciton dissociation at organic donor-acceptor interfaces

Domain: physics
Published: 2025-12-22
Editorial outcome: minor_revision

Abstract
The field of organic photovoltaics has witnessed a renaissance in recent years owing to the development of non-fullerene acceptor materials reaching record power conversion efficiencies of >20%. New computational models are needed to rationalise the regimes of photophysics reached in these materials.

Here we report on a novel implementation of eXcitonic state-based Surface Hopping, a powerful non-adiabatic molecular dynamics method for the simulation of photo-induced charge generation in truly nanoscale donor-acceptor interfaces. We observe a transition from an inefficient cold to an efficient hot exciton dissociation mechanism as the electronic coupling between the molecules or the dielectric constant of the material is increased. The hot pathway is observed to occur by Frenkel excitons converting into transiently delocalised hybrid exciton-charge transfer states that subsequently form charge separated states. This avoids the formation of kinetically trapped interfacial charge transfer states that are prone to non-radiative recombination.

Introduction
Organic solar cells (OSCs) have attracted considerable research attention in recent decades, largely due to their relatively low processing costs and mechanical flexibility allowing applications that cannot be accessed using inorganic materials1. The introduction of a two-component OSC or heterojunction, where donor and acceptor phase are separated by an interface, constituted a seminal leap towards the commercial application of OSCs2, whilst the introduction of non-fullerene acceptors (NFAs)3, e.g., Y64 and its derivatives, further advanced the limits of the power conversion efficiency of OSCs, with recorded efficiencies recently exceeding 20% for the first time4. Recent progress in overcoming the stability issues associated with heterojunctions5 has also been made via the synthesis and characterisation of single-component OSCs6,7.

Theory and computation are lagging behind those recent experimental developments leaving a large gap in our fundamental understanding of charge generation in OSCs8. Whilst the initial formation of bound electron-hole pairs or Frenkel excitons upon light absorption in low dielectric organic semiconductor (OS) materials is generally agreed upon, their subsequent conversion to free charge carriers is less clear. Two principal mechanisms have been postulated in the literature: the so-called cold exciton dissociation mechanism, where an exciton dissociates at the donor-acceptor interface into a low-energy (cold) charge transfer (CT) state (in the following denoted interfacial CT (iCT) state) that needs to overcome the large Coulomb barrier via thermal motion to form free charges9–12. In the second mechanism, termed hot exciton dissociation2,13–18, the exciton instead dissociates into electronically excited (hot), non-interfacial CT (niCT) states, which are poised to transition into charge separated (CS) states and free charge carriers.

The use of density functional theory (DFT) has been instrumental in describing opto-electronic properties of organic molecules and donor-acceptor complexes19–22. Moreover, assuming that the dynamics of photo-excited OSs can be treated as a sequence of transitions between localised molecular electronic states with well-defined rate constants (e.g. Marcus or Redfield rates), Kinetic Monte Carlo (KMC) simulations have found some success in reproducing internal quantum efficiencies (IQEs) of certain systems23–25. However, ordered OSs, especially the ones that would permit large IQEs, are in a difficult regime where the assumption of discrete transitions between localised molecular electronic states may no longer be a valid approximation.

Indeed, recent work has shown that electronic state delocalisation26–37, the symmetry of electronic band states (or, equivalently, the sign-pattern of electronic couplings)38,39, and the feedback from the electronic to the nuclear subsystem40–42 are essential for a realistic description of charge carrier26–29,38,39 and exciton transport dynamics30,31,35,40,42 in ordered OSs - these are quantum effects that are missing in typical KMC simulations. This is not expected to be any different for charge generation in the latest generation of high efficiency OSC materials, as they tend to form highly crystalline morphologies4. Hence, the development of simulation methods that include these vital effects is highly desirable.

In this respect, a study by Scholes considered the energetics of delocalised eigenstates of a static site-based Hamiltonian, and concluded that the energetic alignment of the exciton and niCT states promotes free charge generation, and avoids the occupation of the kinetically trapped iCT state43. Further support for this finding requires the prediction of dynamics, specifically the propagation of time-dependent populations of the electronic states involved in charge generation, which has since been attempted with a variety of methods with varying degrees of accuracy. Gélinas et al. modelled their experimental evidence of ultrafast charge separation in a fullerene acceptor heterojunction by propagating the electronic state populations of a similar site-based Hamiltonian with a truncated state space, and concluded that charge generation in such a blend is assisted by the transient occupation of high-energy niCT states with a delocalised anion14. A notable development was made by Vukmirović and Janković44, who designed an approach to describe the time-dependent population dynamics of delocalised electronic states, starting from a Hamiltonian constructed from the full manifold of localised singlet exciton/CT-states in a model heterojunction. A further significant advancement, based on a similar formulation of the electronic structure, was recently achieved by Kassal and coworkers, who simulated the population dynamics of delocalised quantum states with an analogue to KMC termed delocalised Kinetic Monte Carlo (dKMC)32–34.

Non-adiabatic dynamics approaches are a powerful alternative to the above methods and are currently among the most accurate approaches for the simulation of the electronic quantum dynamics of photo-excited organic materials. Multi-configurational Time-dependent Hartree (MCTDH)45–47 is generally viewed as the gold standard, but this method is computationally expensive even in its multilayer version. Here we report on the development and application of eXcitonic state-based surface hopping (X-SH), a computationally less expensive quantum-classical non-adiabatic dynamics method enabling realistic atomistic simulation of photo-induced charge generation in truly nanoscale donor-acceptor interfaces (10 − 100nm) on experimentally relevant timescales (>10ps). Solving the electronic Schrödinger equation coupled to the finite temperature nuclear motion, the method naturally accounts for delocalisation of excitons and charge carriers, diagonal and off-diagonal thermal electronic disorder (as well as static disorder if present) and the feedback from electronic to nuclear dynamics. Our method comprises a major improvement over a preliminary version of X-SH that did not account for off-diagonal thermal electronic disorder48 deemed essential for exciton30,31,35,40,42 and charge transport26–29,38,39 and thus for photo-initiated charge generation.

We apply our simulation tool to systematically investigate the effect of important material parameters on the mechanism and efficiency of exciton dissociation—in particular, the electronic coupling between the molecules in donor and acceptor materials (determining their band width and giving rise to charge carrier delocalisation), and the electron-hole interaction strength (related to their dielectric constant and exciton binding energy). To this end we consider a model heterojunction comprised of a crystalline donor, α-hexathiophene (a6T), and a crystalline NFA, perylene-diimide (PDI) whose derivatives have been previously characterised as components of heterojunction OSCs49,50. Our choice to consider an idealised crystalline system reflects the fact that the latest organic photovoltaic devices with record efficiencies are made of crystalline materials4. Electronic coupling and electron-hole interaction strength are initially chosen to correspond to the physical values as obtained by DFT parametrisation. These parameters are then varied within physically realistic bounds for OSs to explore different charge generation regimes while keeping the atomistic structure of the molecules and the packing structure of the interface unchanged. In this way, the effect each parameter has on exciton dissociation can be identified in isolation. We lastly examine the role played by structural disorder by selecting parameters that led to efficient charge generation in the crystalline a6T:PDI interface, and re-simulating the exciton dissociation process with an interfacial defect introduced to disrupt the donor/acceptor stacking pattern.

We find that for the a6T:PDI interface, charge generation proceeds via the cold exciton dissociation mechanism: fast relaxation of the initial band-like electronic excitation to a localised Frenkel exciton, diffusion of the Frenkel exciton to the interface via hopping transport, followed by formation of a vibronically fully relaxed iCT state that slowly dissociates thermally to free carriers. Intriguingly, as we increase the electronic coupling between the molecules we increasingly observe a second, more efficient hot exciton dissociation channel. Here, Frenkel excitons convert to transiently delocalised hybrid exciton-charge transfer states (hybrid XT-CT states) that directly transition into niCT and CS states. Remarkably, the hot exciton dissociation process is observed to occur at distances of up to several nanometers away from the interface owing to the delocalised nature of hybrid XT-CT states. This way, the formation of trapped iCT states that are prone to non-radiative recombination is avoided. Similar observations are made when the dielectric constant of the donor and acceptor is increased in place of electronic coupling. Additionally, for a parameter combination where the hot pathway dominates, the insertion of an interfacial defect into the a6T:PDI heterojunction significantly reduces the rate of exciton dissociation.

Results
Model electronic Hamiltonian for a6T:PDI heterojunction
Quantum chemical methods (e.g., TDDFT) generate electronically excited adiabatic states but they are too time-consuming for non-adiabatic dynamics of large molecular systems. A more feasible strategy is to fragmentise the system into molecules or molecular fragments and to construct the electronic Hamiltonian of the full system in the space of locally excited or quasi-diabatic electronic states, which we refer to as local states in the following. This is illustrated here for the a6T:PDI heterojunction that is schematically shown in Fig. 1a–c. It contains a 1-dimensional electronically active region of 20 a6T donor and 20 PDI acceptor molecules embedded within an extended 3D crystalline environment that ensures the structural integrity of the electronically active system (see Methods for details of structure generation). We consider a minimal electronic state space for this system consisting of (i) the lowest singlet excited states |ϕk XT where the excitation is localised on molecule k while all other molecules are in their ground state, (ii) the charge transfer (CT) states |ϕCT kl, where molecule k carries a positive charge (hole) and molecule l a negative charge (excess electron) while all other molecules are in their ground state. The corresponding electronic Hamiltonian in this locally excited state basis is given by Eq. (2) and shown schematically in Fig. 1d. The diagonal elements of HbXT and HbCT are referred to as site energies (Eqs.(9) and (11)) whilst the off-diagonal elements are the excitonic couplings between XT states (Eq.(12)) and the electronic couplings between the CT states, respectively (Eq.(16)), and HbXT-CT contains the couplings between XT and CT states. While in our present application we only consider singlet excited states, it is relatively straightforward to extend the electronic state space of X-SH to higher spin and higher energetic states if required.

We further assert that charge generation is dominated by excitons formed in the acceptor phase. This is hinted by our previous finding that exciton transport is faster in PDI than in a6T40. Hence, we disregard exciton formation in the a6T phase. The excitonic state of PDI molecule i, |ϕ i XT, is described by the lowest S1 singlet excitation which is dominated by the natural transition orbitals that closely resemble the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) (98.6%, see Fig. 1c in ref. 40). Excitations above the first-excited singlet state are relatively high in energy and are disregarded. The frontier orbitals relevant to Frenkel exciton formation in PDI, and for electron and hole transport in PDI (LUMO) and a6T (HOMO), respectively, are shown in Fig. 1e. In regard to the CT states, TDDFT calculations predict that the interfacial CT state where the excess electron and hole are localised on the interfacial a6T:PDI pair is 390meV below the Frenkel exciton state that is localised on the interfacial PDI molecule (see Methods). The electronic energy of the CT states further away from the interface is based on the interfacial CT state energy plus a screened Coulomb electron-hole interaction term (dielectric constant ϵ r=3.5, a typical value in organic semiconductors51) and an external electric field contribution typically found in OSCs (105Vcm−1)20,52. This results in an exciton binding energy (defined in Eq. (24)) of 260meV, which is typical of NFA materials53. Excitonic and electronic couplings between the respective states are calculated efficiently using the transition charge (TrESP54) and analytic overlap methods (AOM55,56) parameterised to (TD)DFT calculations, respectively. This enables the efficient update of all electronic Hamiltonian matrix elements (site energies and couplings) at each nuclear timestep of the X-SH simulations, thereby accounting for diagonal and off-diagonal thermal electronic disorder (related to the electron-phonon coupling). Relevant parameters defining the electronic Hamiltonian and the diagonal and off-diagonal thermal electronic disorder are given in Table 1, with further details given in the Methods section.

Density of electronic states
We first examine the electronic structure of the a6T:PDI junction at 0K by diagonalising the electronic Hamiltonian at the minimum energy configuration of the junction. This generates the adiabatic electronic eigenstates of the system as a linear superposition of the local XT and CT states, as shown in Fig. 2a (0K distributions in full colours). To classify the eigenstates we first divide the local CT states into iCT, niCT and CS states. In the iCT state the electron-hole pair resides solely at the donor and acceptor molecules forming the interface, in the niCT state between the interface and the maximum of the Coulomb barrier (that results from the applied external field) and in the CS state beyond the Coulomb barrier maximum. The eigenstates are then denoted XT (CT) eigenstates if they project to >95% on local XT (CT) states and as mixed hybrid XT-CT states otherwise. CT eigenstates are then further divided into iCT (CS) eigenstates if they project to >80% on local iCT (CS) states and as niCT eigenstates otherwise.

We find that the centre of the XT band (black) remains at ≈400meV above the lowest energy iCT states (red) with a bandwidth that is 4H k X l T, where H k X l T is the excitonic coupling between PDI molecules in the direction orthogonal to the interface. The niCT and CS bands are at higher energies than the XT band, about 650meV above iCT, with a bandwidth greater than 4½(HHT km + HET ln)=2, where HHT km and HET ln are the electronic couplings for hole and electron transfer in a6T and PDI in the direction orthogonal to the interface, respectively. Evidently, the Coulomb interaction (shown in orange in Fig. 2) results in CT bands that are wider than what one would expect from electronic coupling interactions alone. Notice that the centre of the CS band is below that of the niCT band due to the applied electric field (Fig. 2a solid distributions). Notice also that the density of hybrid XT-CT states is small due to the energetic separation of XT and CT bands resulting in little mixing. MD simulations in the electronic ground state show that all bands become considerably broadened at 300K (Fig. 2a shaded distributions).

Cold exciton dissociation
The time-dependent relaxation of a photo-induced excitation in the electronic state space of the a6T:PDI junction characterised above is simulated using X-SH. We refer to Methods for a description of the implemented X-SH method and for simulation details. Photo-excitation is modelled by a Gaussian laser pulse generating an initial distribution of excitations. The time evolution of the electronic wavefunction Ψ(t) (Eq.(6)) along a representative X-SH trajectory is shown in Fig. 3. After initial excitation (IPR, Eq. (21)) ≈ 5), the relaxation process can be broadly resolved into 4 stages.

First, rapid exciton relaxation to the bottom of the XT band leading to exciton localisation on one to two molecules (IPR=1–2) within 50fs (Fig. 3b and c). The dissipated excitation energy of about 290meV is due to fast downward electronic transitions within the XT band (Fig. 3d and e) and the geometric relaxation related to an internal reorganisation energy of a single PDI molecule in the S1 state of 195meV.

Second, diffusion of the exciton to the a6T:PDI interface on the ps timescale assisted by a transient delocalisation mechanism (Fig. 3f and g). Here, wavefunction contraction and expansion events (Fig. 3h) induced by short-lived thermal excitations within the excitonic band (Fig. 3i and j) eventually bring the exciton into the interfacial region.

Third, exciton dissociation to the iCT state. When thermal fluctuations bring the interfacial excitonic wavefunction (Fig. 3k) in close energetical alignment with the iCT state, the two states start to mix due to the finite electronic coupling between local iCT and local interfacial XT states. This results in the formation of a localised, low-energy hybrid XT-CT state at the interface (Fig. 3l–o) that rapidly converts to the iCT state (Fig. 3p) within 20–30fs. The iCT state remains solely at the bottom of the CT band due to the strong Coulomb interaction, which localises the charges and allows the interfacial a6T:PDI molecules to relax to their equilibrium geometries in the cationic/anionic charge states.

Finally, charge separation or recombination of the iCT state. We observe frequent oscillations of charge density between iCT and niCT states (Fig. 3q) but the electron-hole pair is ultimately unable to overcome the large Coulomb barrier (700meV) between the vibronically fully relaxed iCT state and the CS states on the 10ps timescale of present X-SH simulations and remains at the bottom of the CT band (Fig. 3r–t). The subsequent dynamics, that is charge separation or recombination of the iCT state, involves longer time scales. This dynamics is not included in the present simulations as these longer timescales cannot be presently accessed. The mechanism described above is referred to as cold exciton dissociation because all excitons relax to the lowest excited (cold) state of the system, i.e. iCT state, before they acquire enough thermal energy to separate to free charges or recombine to the ground state.

To describe the relaxation dynamics more quantitatively, we project Ψ(t) onto the local XT, iCT, niCT and CS states and average over all trajectories to obtain the corresponding time-dependent populations P XT, P iCT, P niCT and P CS, (given by Eqs.(17–20)). We find that the decay for XT fits well a bi-exponential model (Fig. 4a, solid black and dashed cyan lines, R2 = 0.9991),

P XT = Ci exp(−ki t) + Cb exp(−kb t), (1)

where ki and kb can be interpreted as rate constants for the dissociation of excitons that are photo-generated in the vicinity of the interface (subscript i) and the bulk region of PDI (subscript b), respectively. The interface-generated excitons can readily dissociate to iCT states without having to diffuse across the PDI phase giving rise to a fast rate constant ki = 3.03ps−1. The bulk-generated excitons decay about a factor of 30 slower, kb = 8.65×10−2ps−1, indicating that diffusion to the interface is rate-limiting.

The contribution of interface and bulk-generated excitons in our simulation, Ci = 0.15 and Cb = 0.82, reflect the 20 times larger number of PDI molecules in the bulk phase compared to the interface. In general, the contributions will depend on the interface surface to volume ratio of the heterojunction. We note that the resolution of the exciton decay profile into interfacial and diffusion-limited components has also been observed in experiment14,57. The decay of the XT population is mirrored by an increase in the CT population where, in accord with the above example trajectory, the vast majority of the CT population on the present simulation timescale is composed of iCT states (Fig. 4a, solid red lines). Only traces of niCT (solid purple) and virtually no CS states (solid blue) are observed.

Hot exciton dissociation
The above results imply that in order to accelerate free charge generation, the formation of iCT trap states should be avoided and excitons should preferentially relax to niCT and on to CS states. In the a6T:PDI heterojunction studied above, this does not happen because the niCT and CS bands are too high in energy compared to the XT band (see Fig. 2a). In the following we consider two strategies to achieve a better energetic overlap between XT and niCT/CS band states: (i) increase of the electronic coupling in the donor and acceptor materials and (ii) increase of the dielectric constant. The newly assigned parameters defining the electronic Hamiltonian after taking approaches (i) and (ii) are shown in Table 2. All other parameters remain the same as in the simulations of the original system above. We find that doubling the electronic couplings (Table 2, simulation 2) and thus the width of niCT and CS bands create tail states that are close in energy with low-energy states in the XT band (Fig. 2b). Alternatively, a small increase in dielectric constant from 3.5 to 5 (Table 2, simulation 4) shifts the niCT and CS bands down in energy to align with the XT band (Fig. 2c). Adopting the common expression for exciton binding energy53 (Eq. (24)) and approximating the formation energy of free charges from the iCT state by a screened Coulomb potential, we find that the increase in dielectric constant corresponds to a lowering of the exciton binding energy from 260 to 10meV (however, we note that this is an approximation as the electron-hole interaction energy at short range may not exactly follow a screened Coulomb interaction). In both cases, with an increase of the electronic coupling or dielectric constant, the resulting close energetic alignment also increases the density of hybrid XT-CT states.

We have carried out X-SH simulations for all parameter regimes denoted simulations 2–7 in Table 2. The population dynamics corresponding to a doubled electronic coupling magnitude (Table 2, simulation 2) is shown in Fig. 4b. We find that the exciton decay remains bi-exponential (R2 = 0.9991) with ki and kb now marginally larger than in the previous simulations with unmodified electronic couplings. Importantly, the CS state population is now 8.0% at 10ps, up from 0.3% for unmodified couplings. If both the electronic couplings are doubled and the dielectric constant is increased from 3.5 to 5 (Table 2, simulation 5), the exciton dissociation rate increases further, generating the CS state as the majority species after 10ps, 33.1%, Fig. 4c. The CS yield is increased further if electronic couplings are tripled (Table 2, simulation 6). However, this high coupling regime is not easily accessed with conventional organic semiconductors.

The cause of the remarkable increase in CS yield is now a pertinent question. Analysing the X-SH trajectories, we find that the improved energetic alignment of XT and CT states brought about by larger electronic couplings and/or higher dielectric constant gives rise to an additional and more efficient hot exciton dissociation mechanism. This is exemplified by the evolution of the wavefunction Ψ(t) along a representative X-SH trajectory for doubled electronic coupling and dielectric constant increased to 5 (simulation 5, Fig. 5). Whilst the initial stages of exciton relaxation are the same as in the cold exciton dissociation mechanism (Fig. 5a), we now observe the formation of a state not seen before, whereby a transiently delocalised XT state (Fig. 5b) forms a delocalised hybrid XT-CT state (Fig. 5c). This state is a superposition of delocalised XT, iCT and niCT states. In such states, the excitons and charges are higher in energy and delocalised to a greater extent (see Fig. 5d–f and Supplementary Fig. 6) and can be centred further from the donor/acceptor interface (Supplementary Fig. 8). It is then converted into an niCT state (g), which is able to transition to an energetically accessible CS state (h). As in Fig. 3, the right-hand side panels plot the evolution of relevant observables over the time frames of the events shown on the left. This includes: exciton IPR (d); active eigenstate potential energy relative to the lowest-energy state in the relevant band (e, j); active eigenstate position relative to the lowest-energy state in the band (f, k); electron and hole IPR (i). The vertical dashed lines correspond to observable averages over all time steps and all trajectories. The electronic wavefunction Ψ(t) is classified as a CS state after 10ps of X-SH dynamics.

A remarkable feature of the hot exciton dissociation mechanism is that a significant proportion of excitons no longer dissociate at the interface, as traditional models of excitonic solar cells would assume2, but several lattice spacings/nanometres away from it, owing to the delocalised nature of hybrid XT-CT states in this regime (see Supplementary Figs. 4 and 10, using the center of the wavefunction as a distance measure). Notice that the iCT states are still (on average) the energetically lowest excited states and they are still formed (Fig. 4c). Thus, in this regime both the cold and the hot dissociation mechanisms occur simultaneously (Supplementary Fig. 5). Yet, the fraction of excitons dissociating via the efficient hot mechanism steadily increases with increasing electronic coupling and/or dielectric constant.

Hopping Model vs eXcitonic state-based Surface Hopping
Charge generation in OSCs is often modelled with KMC or, equivalently, chemical master equations using excitation energy transfer and electron transfer rates for exciton or charge hopping between different molecules. In the following we investigate how the population dynamics obtained from a chemical Master Equation compares with the ones obtained from X-SH simulations. Hopping rates between sites were calculated with semiclassical Marcus theory58, with the relevant electronic parameters (reorganisation energies, electronic or excitonic coupling and driving force) chosen to be the same as in X-SH simulations. Full details of the calculations are given in Supplementary Discussion 1. The results for the original system where only the cold exciton dissociation mechanism operates (simulation 1 in Table 2) is shown in Fig. 4a (dashed lines). Although the exciton dissociation rate is somewhat underestimated compared to X-SH, the master equation approach faithfully models the population dynamics of CT states in a, but the large discrepancy in niCT and CS populations in regimes shown in b, c, as the Master equation approach does not account for the formation of hybrid XT-CT states and charge delocalisation.

Effect of interfacial defect
So far, we have focussed on the dynamics of free charge photo-generation in a perfectly ordered heterojunction, as a crystalline morphology of the donor/acceptor aggregates has proven beneficial to charge generation efficiency4 and increased the role played by delocalisation26–37. Here, we introduce a prototypical interfacial defect into the PDI acceptor phase (see Supplementary Fig. 15), to examine its potential effect on exciton dissociation and the dynamics of the CT-states. We focus on the dynamics of a system where free charges are being mostly generated by hot exciton dissociation (simulation 5 in Table 2), as we expect the hot pathway to be most affected by the introduction of a defect. The change in population dynamics upon the introduction of the defect is displayed in Fig. 4d. Exciton dissociation has slowed down markedly, as the disruption in a6T-PDI stacking at the interface has reduced the coupling between the XT and CT states and reduced the density of hybrid XT-CT states, which are necessary for exciton dissociation (see Fig. 3). This results in a significantly lower occupation of CT-states throughout the dynamics. Moreover, a smaller fraction of excitons now dissociate via the hot pathway, whilst a more significant fraction dissociate via the cold pathway that is prone to non-radiative recombination. The interfacial defect stifles the hot pathway by disrupting electron transfer between PDI molecules near the interface. This reduces the density of hybrid XT-CT states with charges delocalised further away from the interface, which are instrumental to hot exciton dissociation (see Fig. 5).

Comparison with other theories
The standard approach for the simulation of charge generation in organic interfaces, KMC employing Marcus rates, has seen success in comparison to experiment23,25, but this approach can only be applied to electronically weakly coupled systems, e.g. structurally disordered materials, where excitons and charges are localised on molecules or molecular units. It cannot be justifiably applied to ordered, electronically well connected systems where delocalised electronic states are likely to dominate the relaxation dynamics (as shown above). Whilst we have tackled this problem with X-SH non-adiabatic dynamics simulation, Kassal and coworkers took a different path and developed a compelling extension to KMC coined dKMC32–34. Starting with a fully quantum mechanical Hamiltonian description in terms of a local electronic basis and harmonic nuclear modes, they utilised the polaron transformation to encode any strong system-bath coupling into the electronic subsystem, and propagate the resulting electronic density matrix via the Redfield59 equation. By integrating this piecewise and mapping it to KMC, they propagate the evolution of excitation density through transitions between the eigenstates of the Hamiltonian of the full system, thereby accounting for delocalisation.

In their recent application of dKMC, Balzer et al. made qualitative observations that are similar to those reported in our current work34. They observed that excitons may dissociate fairly remotely from the interface (thereby directly accessing niCT states) as well as an increased incidence of hybridised states with increasing coupling strength. Their earlier work found that even a marginal increase in carrier delocalisation accelerated charge generation from CT states33, which they attributed to an increase in overlap between eigenstates involved in charge separation. The qualitative similarities are striking and reassuring given that dKMC and X-SH are based on very different theoretical foundations. An advantage of X-SH is its full atomistic nature allowing for a most realistic simulation of molecule-specific donor-acceptor interfaces that includes all vibrational degrees of freedom explicitly without assuming harmonic interaction potentials and without the need of guessing model parameters, of which there are many in such simulations. Importantly, X-SH encompasses a wide range of exciton transport and dissociation regimes (in contrast to Marcus or Redfield theory) which makes this simulation tool widely applicable, especially to the latest generation of NFAs. The drawback of X-SH is the higher computational cost when compared to, e.g., dKMC, but as we have shown here X-SH simulations on the 10–100ps timescale of truly nanoscale systems are feasible.

Besides dKMC, fully quantum-dynamical methods have been employed to study charge separation in simplified donor-acceptor interface models. Previous studies have shown that the quantum nature of high-frequency vibrational modes can, in some cases, influence the wavepacket dynamics, as evidenced by coherent beatings in the populations60. These coherent oscillations can affect relaxation processes and charge separation on ultrafast timescales. From our previous comparisons with fully quantum-dynamical simulations (using MCTDH for the same electronic Hamiltonian), we found that such initial dynamics are captured with reasonable accuracy by X-SH, even though the high-frequency modes are treated classically. An important benefit of our approach over a fully quantum-mechanical treatment, beyond the significant computational efficiency, is that our scheme correctly recovers the long-time detailed balance, a fundamental requirement for processes occurring on extended timescales where essentially all vibrational modes are in thermal equilibrium and vibronic coherences are significantly damped48.

Comparison with experiment
Our simulations reproduce salient experimental observations of the manner of exciton dissociation. This includes the resolution of exciton dissociation into rapid interfacial and diffusion-controlled regimes14,57, and the initial relaxation of the exciton largely outpacing charge separation61. A valuable insight into the charge generation timescales predicted by X-SH can be gleaned through a comparison to the work by Bakulin et al.13 They utilised a secondary, IR-wavelength pump pulse following a standard excitation pulse to excite charge carriers located in bound CT states at the heterojunction interface, measuring the subsequent increase in photocurrent. Bakulin et al.13 conducted this measurement on two NFA systems (among other fullerene acceptors), which showed that exciton dissociation occurs on the timescale of several ps, as in several other works16,17,57,62,63. They explained the increase in photocurrent through the transient occupancy of higher-energy delocalised CT states that preceded total charge separation, which occurred on the timescale of 10–50ps. Qian et al.10 also measured the rate of exciton dissociation in a number of low-offset NFA blends, using transient absorption spectroscopy where the resulting signal of the probe pulse was decomposed into excitons and distinct charge carriers. Significant exciton dissociation occurred in most blends within 10ps. Thus, X-SH can qualitatively reproduce the aforementioned experimental timescales of exciton dissociation and charge separation, with an electronic Hamiltonian solely parameterised with respect to DFT, that explicitly incorporates the effect of fluctuating couplings between electronic states.

As mentioned earlier, we only considered singlet states in this work, disregarding any potential formation of low-lying triplet states. Gillett et al. provided a comprehensive experimental assessment of the role of triplet CT states in low-offset and NFA-based OSCs64. They can be close in energy to their singlet counterparts because the exchange energy decreases rapidly with increasing electron-hole separation65. Triplet CT states may form via a geminate spin-mixing mechanism from singlet CT states on nanosecond timescales; however, Gillett et al. found no experimental evidence for this geminate pathway in the NFA blends they examined. Instead, they identified a dominant non-geminate mechanism, in which the recombination of free charges generates triplet CT excitons that can undergo fast back charge transfer to Frenkel exciton triplets. However, the non-geminate recombination pathway requiring free charge carriers happens on much longer timescales than the current X-SH simulations can access. Yet, depending on the material under consideration, triplet formation may occur faster, in which case it provides additional entropic driving force for free charge generation66,67.

Design rules
Finally, in light of the results obtained from current X-SH simulation, we discuss favourable materials characteristics for organic photovoltaics that may assist experimentalists in selecting suitable combinations of donor and acceptor materials. 1. Small exciton binding energy, Eb ≈ 0meV, brought about by, e.g., an increase in dielectric constant to >5. This has also been predicted through experiment68, but with an emphasis on the prevention of voltage loss. X-SH has shown that in addition to facilitating the escape of electron-hole pairs from energetically low lying iCT states, a reduction in Eb also signifies increased energetic overlap between the lower exciton band edge states and niCT states. 2. Large electronic couplings to generate delocalised CT states, HHT km > λHT i /2, HET ln > λET i /2. It has also been shown experimentally that increased carrier mobility can increase OSC efficiency69,70. X-SH finds that an additional effect of the increased CT-bandwidth is an improvement of its overlap with the exciton band edge resulting in a higher density of delocalised hybrid XT-CT states and a more effective transition from the latter states to niCT and CS states. 3. Large excitonic couplings to generate delocalised XT states, HXT kl > λXT i /2. In addition to the increased density of hybrid XT-CT states, similarly as in 2. above, and a more effective transition to these states from pure excitonic states. 4. Structurally ordered donor/acceptor aggregates. The introduction of an interfacial defect that disrupts electron transfer in the PDI acceptor phase near the interface resulted in the cold exciton dissociation pathway gaining significance with respect to the hot pathway, leaving electron-hole pairs more likely to geminately recombine after exciton dissociation. 1.-4. will generally suppress the formation of kinetically trapped iCT states and increase the contribution of hot exciton dissociation. Lastly, a small energy gap between the LUMO and higher lying unoccupied orbitals of the acceptor can further assist free charge generation as suggested in ref. 71. In the present system the energy gap between LUMO and LUMO+1 was relatively large (1.8eV), which is why it was not included in our simulations. However, when this gap is small, the entropic advantage is expected to further enhance the beneficial effects of strong electronic coupling and a small exciton binding energy, in comparison to the results reported in this study.

Outlook
X-SH strikes an optimal balance between theoretical rigour and computational efficiency. An important advantage over other model Hamiltonian approaches is the atomistic resolution of the heterojunction, whose electronic Hamiltonian is parameterised with respect to DFT, and then propagated on-the-fly using ultrafast site energies and electronic couplings. The thermal fluctuations and the transport physics of organic semiconductors are naturally included in the literature. This facilitates a realistic and physically intuitive comparison of the dynamics between different OSC materials, as any difference in performance can be directly attributed to electronic parameters intrinsic to the molecular species or its morphology. The current computational expense of X-SH still limits its practical use to about 40 molecules per donor or acceptor phase due to the quadratic scaling of the size of the CT block with number of molecules. While this is commensurate for the modelling of large 1D or small 2D systems it currently prevents the modelling of 3D systems. This limitation will be addressed in future work. The current X-SH simulations are focused on the 10ps timescale. Non-radiative recombination to the ground state has not been included in the present simulations because, according to simple rate-based estimates, the recombination kinetics is expected to occur on significantly longer timescales2. However, as the simulation time that can be accessed by X-SH increases, the recombination to the ground state needs to be included such that the competition between charge separation and recombination can be leveraged to yield IQEs for atomistic heterojunction structures.

Methods
Excitonic State-based Surface-hopping (X-SH)
As described in the section Results and Discussion, the total electronic Hamiltonian in X-SH is expressed as a sum of Frenkel exciton (XT), charge transfer (CT) and XT-CT coupling Hamiltonians,

Hb (R) = HbXT (R) + HbCT (R) + HbXT-CT (R), (2)

where the dependence on the (classical) nuclear coordinates R is explicitly stated. The XT block is given by

HbXT = ∑ M XT k=1 HXT kk (R) |ϕX k T ⟩⟨ϕX k T | + ∑ M XT k=1 ∑ M XT l≠k HXT kl (R) |ϕX k T ⟩⟨ϕX l T | (3)

where HXT kk (R) and |ϕXT k = |ϕXT k (R) respectively refer to the (quasi-)diabatic site energy and wavefunction of XT state k, in which molecule k is in the lowest singlet excited state (Frenkel exciton) and all other molecules are in the ground state, and HXT kl (R), k ≠ l, refers to the electronic coupling between two diabatic XT states, where the excitons are localised on different molecules k and l. Here, the number of exciton states, M XT, is equal to the number of acceptor molecules, M A, M XT = M A.

The CT block is given by

HbCT = ∑ M D k=1 ∑ M A l=1 HC kl T (R) |ϕC kl T ⟩⟨ϕC kl T | + ∑ M D k,m=1 ∑ M A l,n=1 HC kl T,mn (R) |ϕC kl T ⟩⟨ϕC m T n | (4)

where HC kl T (R) and |ϕC kl T (R) refer to the diabatic site energy and wavefunction of CT state kl with the excess hole on molecule k of the M D donor molecules and the excess electron on molecule l of the M A acceptor molecules, and HC kl T,mn (R) refers to the electronic coupling between CT states kl and mn. The total number of CT states is then M CT = M A M D. The primes on the right hand side of Eq.(4) indicates that the diagonal term (k = m, l = n) is excluded.

Lastly, the XT-CT block contains the electronic couplings between XT and CT states,

HbXT-CT = ∑ M XT k=1 ∑ M D m=1 ∑ M A n=1 HX k, T m n CT (R) |ϕX k T ⟩⟨ϕC m T n | + ∑ M D m=1 ∑ M A n=1 ∑ M XT k=1 HX m T n,k CT (R) |ϕC m T n ⟩⟨ϕX k T | (5)

Insertion of the electronic wavefunction
Ψ(t) = ∑ M XT k=1 uX k T (t) |ϕ k XT ⟩ + ∑ M D m=1 ∑ M A n=1 uC m T n (t) |ϕ m CT n ⟩ (6)

into the time-dependent Schrödinger equation, gives the time evolution of the electronic states in terms of the diabatic expansion coefficients u i,

iℏ d u i (t) / dt = ∑ M j=1 u j (t)(H ij (R) − iℏ d ij ), (7)

where for ease of notation we have introduced u i and |ϕ i, u i ∈ {uX k T}, {uC m T n} and |ϕ i ∈ {|ϕX k T}, {|ϕC m T n}, k= 1,…,M XT , m= 1,…,M D , n=1,…,M A , i=1,…,M, M=M XT + M D M A. H ij denotes the matrix elements of the total electronic Hamiltonian (Eq.(2)); the non-adiabatic coupling elements between the quasi-diabatic states, d ij = ⟨ϕ i | d/dt |ϕ j ⟩, are generally small (see ref. 40 for a discussion) and are neglected.

In surface hopping, the nuclear force on only a single active adiabatic electronic state (index a) is required, FI,a, which is given by ref. 72

FI,a = −[U†(∇ I H)U] aa (8)

where ∇ I = ∂/∂ R I, and U is the unitary transformation diagonalising H in Eq. (2). At each MD time step, the probabilities for a surface hop from the active electronic state a to the other adiabatic states j are calculated according to Tully’s expression73, which contains the non-adiabatic coupling element (NACE) between the adiabatic electronic states, d ja = ⟨ψ j | ∇ R ψ a ⟩ · Ṙ, where the NACE between the quasi-diabatic states (d kl ≈ [D] kl) are neglected. Energy conservation is ensured by only allowing trajectories to hop if their kinetic energy along the direction of the non-adiabatic coupling vector (NACV) is equal to or larger than the difference between the two potential energy surfaces (E j − E a), with the nuclear velocities being reversed in the direction of the NACV if the hop is rejected74,75. We adopt the same modifications to the standard FSSH algorithm as in our previous works76 to ensure internal consistency and (approximate) detailed balance in addition to total energy conservation. This includes an energy-based decoherence correction77, a state-tracking algorithm to prevent trivial crossings, and the elimination of spurious long-range exciton and charge transfer. We refer to ref. 76 for a detailed discussion of these modifications and their physical underpinnings.

Site energies
In practice, we define all site energies with respect to the lowest energy state of the electronic Hamiltonian, which is the iCT state |ϕC 11 T . Notice that the electronic ground state (S0) is not included in our model because it will be populated on timescales that are much longer than the present X-SH simulations. The site energy of XT states is then given by

HXT kk (R) = ϵXT + ΔEX k T (R), (9)

where ϵXT is the electronic energy of the Frenkel exciton at the interface |ϕXT k relative to the one for the iCT state |ϕCT 1 11 in the minimum energy configuration of the electronic ground state, R GS, ϵXT = E X k T (R GS) − E C 11 T (R GS). ϵXT is the electronic energy of the Frenkel exciton at the interface relative to the one for the iCT state in the minimum energy configuration of the electronic ground state, R GS. ΔEX k T (R) = EX k T (R) − EX k T (R GS). The latter term takes into account the thermal fluctuations of the site energy related to the local electron-phonon coupling. ϵXT is calculated for a gas-phase interfacial a6T:PDI dimer, whose geometry was extracted from the crystal structure of the heterojunction. The adiabatic excited state energies of the dimer were first calculated with TDDFT, and the energy difference between the localised PDI exciton and interfacial CT-state was subsequently obtained with the MS-FED-FCD78,79 diabatisation procedure. This gave ϵXT = 390meV. For calculation of EX k T (R), the electronic ground state potential energy for a6T molecules and for PDI molecules l≠k were described by the standard Generalised Amber Force Field (GAFF)80, while for the excited PDI molecule k the equilibrium bond lengths were adjusted to reproduce the internal reorganisation energy for S1 exciton transfer between two PDI molecules from (TD)DFT calculations, λS , S = XT. All other parameters for the excited PDI molecule k remain the same as in the electronic ground state. The internal reorganisation energy was calculated in the gas phase using the 4-point scheme:

λS i = [E A (RB) + E B (RA)] − [E A (RA) + E B (RB)] (10)

where E A and E B denote the potential energy of the ground state (A=C) and excited state (B=S1) of a PDI molecule in vacuum, respectively, and RA and RB are the nuclear coordinates at the minimum of the S1 and S0 potential energy surfaces, respectively. The ground and excited state calculations for ϵXT and λXT i were carried out using the CAM-B3LYP81 functional and the 6-31G(d,p) basis set as implemented in the Gaussian1682 electronic structure software package.

The site energy of CT states is given by

HC kl T (R) = − ϵ 1 r kl − 1 r 11 − |E|(r kl − r 11) + ΔEC kl T (R) (11)

where the first term on the right hand side of Eq.(11) describes the Coulomb interaction between the electron-hole pair located on PDI molecule l and a6T molecule k at configuration R, relative to the Coulomb interaction of the iCT state (l=1, k=1) at configuration R GS. r kl is the electron-hole distance between the centres of mass of the respective molecules and ϵ r is the dielectric constant. The second term on the right hand side of Eq. (11) describes the interaction of the electron and hole with a static electric field in the direction orthogonal to the interface and of magnitude equal to |E|, relative to the electric field interaction in the iCT state, again taken at configuration R GS. The field term is introduced to account for the bias due to the presence of external electrodes in OSC devices. The magnitude of the electric field was set to 105Vcm−1, a value that is typical for OSC devices20,52. The third term on the right hand side of Eq. (11) accounts for thermal fluctuations of the site energies and is defined as the electronic energy of state |ϕC kl T at a nuclear configuration R relative to the energy of that state at R GS. ΔEC kl T (R) = EC kl T (R) − EC kl T (R GS). Note that since the iCT state is used as a baseline, HCT 11 (R) = ΔECT 11 (R) throughout X-SH dynamics. Analogous to XT states, for the calculation of EC kl T (R), the electronic ground state potential energy for a6T molecules and PDI molecules (m≠k, n≠l) was described by the standard Generalised Amber Force Field (GAFF)80, while for the positively charged a6T molecule k and the negatively charged PDI molecule l the respective equilibrium bond lengths were adjusted to reproduce the internal reorganisation energy for hole transfer (HT) between two a6T and electron transfer (ET) between two PDI molecules from DFT calculations, respectively. All other parameters for the charged a6T and PDI molecules remain the same as in the electronic ground state. The internal reorganisation energy was calculated in the gas phase using the 4-point scheme (Eq.(10)) where S = HT, ET. Here, E A and E B denote the potential energy of the ground state for the charged molecule (A=C) and neutral molecule (B=N) in vacuum, respectively, and RA and RB are the nuclear coordinates at the minimum of the potential energy surfaces in the C and N state, respectively. The DFT calculations for λHT and λET were carried out in CP2K83 with the B3LYP functional84 and the 6-31G(d,p) basis set. The thermal fluctuations of the first two terms in Eq.(11) (which are part of ΔEC kl T (R)) are very small due to negligible centre of mass motion (r kl) of the molecules in our assembly and they are thus neglected.

The nuclear gradients of the site energies, required for the calculation of nuclear forces on the adiabatic electronic states (Eq.(8)) and NACVs, are obtained from the force field. Note that in Eqs.(9) and (11) only the electron-phonon terms, ΔEX k T (R) and ΔEC kl T (R), respectively, depend on nuclear coordinates whereas the electronic terms are coordinate independent.

Coupling matrix elements
Excitonic couplings between XT states are calculated using the transition charge approach54,

HXT kl (R) = ∑ K ∑ L qT K qT L / |RK − RL| (12)

where K and L label atoms in PDI molecules k and l, respectively, qT K and qT L are the TrESP transition charges and RK, RL the atomic positions. The TrESP charges for PDI were taken from our previous work40 where they were obtained from CAM-B3LYP calculations in Gaussian using the 6–31G(d,p) basis set. The transition charges were frozen during X-SH dynamics. We have shown in previous work that the TrESP approach gives very accurate results when compared to full excitonic couplings (that also include exchange and overlap contributions). Moreover, it was shown that the small fluctuations of transition charges during dynamics can be safely neglected. The nuclear gradients of the excitonic couplings, which are required for the calculation of the nuclear force on the adiabatic electronic states (Eq.(8)) and for NACVs, are obtained analytically from Eq.(12).

As to the calculation of electronic coupling between CT states, we note that in a given CT state |ϕC kl T, the hole and electron are located on distinct molecules, a6T molecule k and PDI molecule l, respectively. Hence, we assert that the CT-state can be written as a single product of localised electron and hole wavefunctions, |ϕC kl T ≈ |ϕh k ⟩ |ϕe l ⟩. The coupling matrix elements between two CT states |ϕC kl T and |ϕC m T n are then approximated as follows,

HCT kl,mn = ⟨ϕCT kl | HbCT |ϕCT mn ⟩ (13)

≈ ⟨ϕh k ϕe l | HbCT |ϕh m ϕe n ⟩ (14)

≈ ⟨ϕh k | Hbh |ϕh m ⟩ δ ln (1 − δ km) + ⟨ϕe l | Hbe |ϕe n ⟩ δ km (1 − δ ln) (15)

= HHT km δ ln (1 − δ km) + HET ln δ km (1 − δ ln) (16)

According to Eq.(15) only those CT states couple that differ in the occupation of either the hole (k≠m, l = n) or the electron (k = m, l≠n) but not in the occupation of both particles. The matrix element between the latter CT states (k≠m, l≠n) is very small compared to the former cases and is neglected85. The integration yields the electronic Hamiltonian of the donor phase hosting the hole, Hbh, and the electronic Hamiltonian of the acceptor phase hosting the electron, Hbe. The electronic couplings for hole transfer in the donor phase, HHT km, and for electron transfer in the acceptor phase, HET ln, are calculated using the ultrafast Analytic Overlap Method (AOM)55,56. In AOM, electronic couplings are estimated via the linear relationship H ab = C′ S′ ab where S′ ab is the orbital overlap between the relevant frontier orbitals of electron-donating and electron-accepting molecules in vacuo (HOMO-HOMO for hole transfer, LUMO-LUMO for electron transfer) and C′ is a constant of proportion determined by fitting the above linear relationship to electronic coupling values obtained from DFT calculations and projection operator-based diabatisation (POD(PBE), 6-31G(d,p))86.

AOM is highly efficient because S′ ab is calculated analytically along the non-adiabatic dynamics trajectories for each molecular pair using frozen orbital expansion coefficients, see ref. 26. This method has been extensively benchmarked in our previous works and was found to give good accuracy on a test set from the HAB79 dataset87 with a mean relative unsigned error of 24.7% (for the range of couplings relevant to this work (10–100meV)), while reducing the computational cost by 4–5 orders of magnitude compared with DFT (POD) calculations56. The scaling constants C′ HT and C′ ET were taken to be the scaling constant (9.463 and −9.463eV, respectively) obtained from fitting against the HAB79 dataset87. The goodness of fit was tested by sampling 100 configurations along a 1ps X-SH run at 300K and comparing the resultant AOM couplings against DFT(POD) calculations. Both couplings displayed good agreement between POD and AOM for the chosen constant (R2=0.904 for C′ HT and R2=0.963 for C′ ET, Supplementary Fig. 14).

The electronic coupling between XT and CT-states, H k X , T m n CT, is limited to the donor-acceptor molecular pair at the interface owing to its fast exponential decay. All other elements in the HbXT-CT block are also calculated at each timestep, but are either close or equal to zero. The XT-CT coupling corresponds to an excited state hole transfer from the interfacial S1PDI molecule to the interfacial a6T molecule and can thus be estimated using AOM, similarly to the procedure described above. The constant of proportion, C′ XT-CT, is obtained by fitting a linear relationship between the overlap of the HOMO orbitals of interfacial a6T and PDI to excited state electronic coupling values obtained from TDDFT calculations and MS-FED-FCD diabatisation78,79 (CAM-B3LYP, 6-31G(d,p)). A fit on 25 configurations sampled along a 1ps X-SH trajectory at 300K gives a value C′ XT-CT = 17.06 eV (R2 = 0.618), see Supplementary Fig. 14. The nuclear gradients of all electronic coupling matrix elements are obtained numerically via finite differences using AOM. We finally note that the signs of the off-diagonal elements in the electronic Hamiltonian are, by construction, consistent with the phases of the corresponding molecular orbitals, which are in-turn dependent on the atomic ordering of such molecules. We refer to ref. 38 for a detailed discussion.

Simulation protocol
An atomistic model for the a6T:PDI interface is constructed from the experimental crystal structures of a6T and PDI. Our aim here is to build a plausible and thermally stable interface structure where molecules adopt the single crystal packing structure some distance away from the interface. To this end, initial atomic coordinates of the a6T and PDI single crystals were taken from the Cambridge Crystallographic Data Centre with structural identifiers ZAQZUM01 and LENPEZ01, respectively. The PDI supercell (containing 240 molecules) was constructed from 20×6×1 unit cells, whilst the a6T supercell (300 molecules) was constructed from 10×15×1 unit cells. The interface was then created by the manual alignment of each supercell along the b and c-directions, followed by their combination along their respective a-directions, such that the π−π stacking between interfacial molecules is maximised, resulting in a distance between the centers of mass of interfacial a6T and PDI molecules of 9.8 Å. Periodic boundary conditions were then applied along the b and c-directions using lattice constants that are very close to the ones for the respective pristine crystals thereby compensating for the small lattice mismatch. To preserve the heterojunction structure along the non-periodic a-direction, we applied weak restraints to the geometric centres of only the upper/lower rings of each molecule, in both the equilibration and production runs. The dimensions of the simulation cell were 300×300×300 Å3 and kept constant during equilibration and production runs. The initial interface structure was equilibrated for 500ps to 300K in the electronic ground state (GS) using a Nosé-Hoover thermostat. Whilst the orientations of the interfacial molecules slightly deviate from the one in the pristine crystal, the molecules further apart adopt the orientation of the respective pristine crystal as indicated by a low RMSE of 0.228 Å with respect to the experimental crystal structures. To proceed, 1000 nuclear configurations were equidistantly sampled along a subsequent 1ns run in the NVE ensemble and used as initial geometries for X-SH. We note in passing that the generated structure may not be the energetically most stable interface structure - more sophisticated structure generation tools would be required to look at different interface structures but this is beyond the intended scope of the current work.

X-SH simulations were carried out for 7 different parameter regimes summarised in Table 2. The electronically active region in which XT and CT states are propagated was chosen to be a one-dimensional chain of 20 a6T and 20 PDI molecules along the a-direction of the above equilibrated interface. All other molecules in the supercell were treated as electronically inactive and interacted with the active region only via non-bonded interactions. Photo-excitation was modelled by a 10fs Gaussian laser pulse. Each X-SH trajectory was initiated in an excited eigenstate randomly taken from a distribution of eigenstates with a total probability that is proportional to the probability for the incoming photon energy being equal to the energy gap between the ground and excited eigenstate, the corresponding transition dipole moment and the probability of the excited eigenstate being sampled along an electronic ground state MD trajectory of the a6T:PDI interface at 300K. These MD trajectories were used to sample the initial excited state populations. In each of the 500 X-SH trajectories, the chosen initial electronic eigenstate was assigned the initial active adiabatic state, ψ a(t = 0), and taken to be the initial electronic wavefunction, Ψ(t = 0) = ψ a(t = 0). The trajectories were propagated in the NVE ensemble for 10ps using a nuclear timestep Δt = 0.05fs. The electronic equation (7) was integrated with the 4th order Runge-Kutta algorithm (RK4), using an electronic timestep δt = 0.01fs. All X-SH simulations were carried out with our in-house implementation of X-SH in the CP2K simulation package83.

Property calculation
The XT, iCT, niCT and CS populations shown in Fig. 4 are calculated as

P XT = ⟨ ∑ M XT k=1 |uX k T (t)|2 ⟩ N (17)

P iCT = ⟨ |uC 11 T (t)|2 ⟩ N (18)

P niCT = ⟨ ∑ M D k=1 ∑ M A l=1 H(r kl − r 11)H(r max − r kl) |uC kl T (t)|2 ⟩ N (19)

P CS = ⟨ ∑ M D k=1 ∑ M A l=1 H(r kl − r max) |uC kl T (t)|2 ⟩ N (20)

where H(x) is the Heaviside step function modified such that H(x) = 0 when x≤0, r 11 is the centre of mass distance between the interfacial a6T-PDI molecular pair, r max is the distance at which the Coulomb potential (first two right-hand-side terms of Eq.(11)) is at its maximum, and ⟨…⟩ N denotes the average over N = 500 X-SH trajectories.

The inverse participation ratio (IPR) is a measure of the wavefunction delocalisation and is usually defined as IPR=1/(∑ k |u k |4), where u k are diabatic expansion coefficients. The presence of different particle-types in X-SH (XT and CT) precludes the direct use of this definition because the squares of the diabatic expansion coefficients of a given particle type no longer sum to unity, which underestimates the purity of the total electronic wavefunction Ψ(t) with respect to each particle type. Thus we first renormalise the XT contribution to Ψ(t) and insert this in the usual definition for IPR, giving

IPR XT (t) = [ ∑ M XT k=1 |uX k T (t)|2 ]2 / ∑ M XT k=1 |uX k T (t)|4 (21)

To calculate the IPR of the electron in the acceptor phase, we first obtain the probability of an electron being found on an acceptor molecule l and the hole on any donor molecule k, P l = ∑ D k=1 |u kl (t)|2, and then renormalise the CT contribution to Ψ(t) giving

IPR e (t) = [ ∑ M A l=1 P l ]2 / ∑ M A l=1 P l 2 (22)

The IPR of the hole in the donor phase is calculated similarly,

IPR h (t) = [ ∑ M D k=1 P k ]2 / ∑ M D k=1 P k 2 (23)

The exciton binding energy is given by:

Eb = Ef − Eopt (24)

where Ef refers to the energy of the uncorrelated electron-hole pair and Eopt refers to the vertical excitation energy between the ground state and the lower exciton (S1) band edge. Both quantities in Eq.(24) are calculated relative to the iCT state energy, where Ef naturally varies with respect to the Coulomb barrier height. Further details about how each quantity is calculated are included in Supplementary Discussion 3.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an interesting theoretical study that is appropriate for publication in Nature Communications after the authors consider my minor comments about the context of the work. The paper, by studying parameters in a theoretical model (electronic coupling, dielectric constant), shows convincingly how exciton dissociation is much more efficient when pathways that avoid the bound electron-hole pair state are exploited. This is an important design principle, as noted in the paper.\n\nThe paper is appropriate for this journal because the topic is at the current forefront of the field. For example, transient delocalization of excitons has been posited as a new design principle to aid exciton transport: https://pubs.acs.org/doi/10.1021/acs.jpclett.2c01133\n\nThe work, while new and novel in terms of the approach, does echo prior studies, which the authors may consider noting more explicitly to support their motivation and conclusions. Specifically, it is well recognized that the interfacial charge-separated state is a kinetic trap. In this paper https://pubs.acs.org/doi/10.1021/nn700179k the basis for a band of CT states that provide pathways with diminished exciton binding is presented. That picture is close to the basic starting point for the present study. The present study significantly expands the scope by predicting dynamics.\n\nThis well-known paper https://www.science.org/doi/10.1126/science.1246249 gave experimental evidence for the advantages of initial charge separation producing a delocalized anion state in the fullerene domain.\n\nThis paper https://pubs.acs.org/doi/abs/10.1021/jp5098102 is closely related to the present study and points out the entropic advantage of a large density of CT states for preventing recombination.", "reviewer_2": "The authors present a detailed theoretical study of nanoscale donor-acceptor interfaces in a prototypical a6T:PDI system, using the X-SH framework to simulate non-adiabatic molecular dynamics. The results indicate that efficient dissociation of excitonic states into free carriers occurs via transiently delocalized hybrid exciton-charge transfer states. Strong electronic couplings between the donor and acceptor molecules promote ‘hot exciton’ dissociation, a similar effect to increasing the dielectric constant of the heterojunction. The authors also discuss the experimental implications of their findings. Additionally, the article introduces an improved X-SH methodology (accounting for off-diagonal disorder) that could benefit future applications to various materials. Given the significance of these findings, I recommend this work for publication in Nature Communications following a revision addressing the comments below. No further review is necessary.\n\n1. The authors consider for their simulations idealized crystalline interfaces between donor and acceptor phases. Indeed, while experimentally it was found that crystalline interfaces work better than amorphous ones for charge separation efficiency, I wonder what the role of defects is. To what extent some typical defects, e.g. faults in the stack of both materials, affect the efficiency of charge transfer mechanisms. These modification can be easily incorporated and tested on the level of X-SH model. This discussion would be an important addition to the paper and point to experimentalists how defect tolerant the charge separation is, for example, for a hot-exciton pathway.\n\n2. There are lots of experimental and theoretical studies in the field pointing to the important of vibronic coherence that can fundamentally affect many photoinduced processes and excited state dynamics. The authors use surface hopping dynamics which fundamentally lacks vibrational wave functions, and decoherence corrections are just fixing the non-adiabatic algorithms. The authors need to include a brief discussion on this phenomenon (which to me could be very important for hot exciton dissociation pathway) and how it can alter the emerging picture.\n\n3. Cold exciton dissociation. Is it correct to assume that the simulations are limited by 10ps of simulations time. Within this timescale, there was no observed dissociation of the iCT states. The subsequent dynamics (i.e. “charge separation or recombination of iCT state”) involves longer timescales and this dynamics is not included in the kinetic models presented. Perhaps some clarifying remarks would be helpful here.\n\n4. Hot exciton dissociation. I would question validity of the statement “a small increase in dielectric constant from 3.5 to 5 (Table 2, simulation 4) results in a decrease of exciton binding energy from 260 to 10 meV.” This is a huge claim and I do not thinks that just referring to Eq. 24 provides justice. Exciton binding energy is a delicate quantity to calculate. I would request the authors to substantiate their claims. Experimental evidence is indirect (ref. 61).\n\n5. Regarding the abstract statement: “Our results highlight the importance of electronic couplings, an often overlooked property, for opto-electronic charge generation.” I believe that this is not true and tuning of electronic coupling has been an active subject of the discussion in the community for decades across both theoretical and experimental studies. Perhaps this statement and similar tone in the main body of the article need to be softened.\n\n6. The present approach considers singlet states only. However, I would argue that weakly bound triplet charge transfer states are iso-energetic with their singlet counterparts and can be easily form even in the absence of strong spin-orbit coupling. Their presence and participation in the process can add additional entropic driving force further improving charge separation. I would suggest to add discussion of this scenario.\n\n7. Code availability statement says: “The custom codes used for this study are available from the corresponding authors upon request.” It would be really helpful to clarify releases of the code for public, for example, if any of the standard version of X-SH are available and what are the plans for realizing the improved version. This is an important advance of this study.\n\nI do not think that X-SH code used for this study was released, see my point #7"}, "author_response": "RESPONSES TO REVIEWER COMMENTS\nThis is a document detailing our point-by-point response to each comments by the reviewers.\nText in black corresponds to the reviewers’ comments repeated verbatim. Blue text corresponds\nto our responses to these comments, while red text corresponds to adjustments that were\naccordingly made in the manuscript.\nReviewer #1 (Remarks to the Author):\nThis is an interesting theoretical study that is appropriate for publication in Nature\nCommunications after the authors consider my minor comments about the context of the work.\nThe paper, by studying parameters in a theoretical model (electronic coupling, dielectric\nconstant), shows convincingly how exciton dissociation is much more efficient when pathways\nthat avoid the bound electron-hole pair state are exploited. This is an important design\nprinciple, as noted in the paper.\nWe thank the reviewer for the very positive remarks on our work and for providing additional\ncontext to our simulations.\nThe paper is appropriate for this journal because the topic is at the current forefront of the field.\nFor example, transient delocalization of excitons has been posited as a new design principle to\naid exciton transport:\nhttps://pubs.acs.org/doi/10.1021/acs.jpclett.2c01133\nThis work has been cited in the original manuscript (Ref. 37).\n1. The work, while new and novel in terms of the approach, does echo prior studies, which the\nauthors may consider noting more explicitly to support their motivation and conclusions.\nSpecifically, it is well recognized that the interfacial charge-separated state is a kinetic trap. In\nthis paper\nhttps://pubs.acs.org/doi/10.1021/nn700179k\nthe basis for a band of CT states that provide pathways with diminished exciton binding is\npresented. That picture is close to the basic starting point for the present study. The present\nstudy significantly expands the scope by predicting dynamics.\nThe suggested paper considers a static Hamiltonian in the basis of localised Frenkel exciton\nand charge-transfer states and indeed highlights the possible advantage of the energetic\nalignment of the exciton and CT bands, due to a complementary combination of the strength of\nCoulombic electron-hole interactions and electronic couplings. However, as the reviewer\npointed out, it is a purely static perspective and lacks the important quantum dynamics that we\nprovide in our paper. We have therefore expanded our Introduction section on page 3:\n“In this respect, a study by Scholes considered the energetics of delocalised eigenstates of a\nstatic site-based Hamiltonian, and concluded that the energetic alignment of the exciton and\nniCT states promotes free charge generation, and avoids the occupation of the kinetically\ntrapped iCT state.[43] Further support for this finding requires the prediction of dynamics,\nspecifically the propagation of time-dependent populations of the electronic states involved in\ncharge generation, which has since been attempted with a variety of methods with varying\ndegrees of accuracy. Gelinas et al. modelled their experimental evidence of ultrafast charge\nseparation in a fullerene acceptor heterojunction by propagating the electronic state\npopulations of a similar site-based Hamiltonian with a truncated state space, and concluded\nthat charge generation in such a blend is assisted by the transient occupation of high-energy\nniCT states with a delocalised anion.[14] A notable development was made by Vukimirovic and\nJankovic [44], who designed an approach to describe the time-dependent population dynamics\nof delocalised electronic states, starting from a Hamiltonian constructed from the full manifold\nof localised singlet exciton/CT-states in a model heterojunction.”\n2. This well-known paper https://www.science.org/doi/10.1126/science.1246249\ngave experimental evidence for the advantages of initial charge separation producing a\ndelocalized anion state in the fullerene domain.\nWe thank the reviewer for this reference. This article was referred to in the original manuscript\nwhen we referenced experimental evidence of the hot exciton dissociation pathway. The\nreviewer is however correct that the article should be mentioned with more emphasis due to\ntheoretical models therein directly building on those constructed by Scholes and explicitly tying\nits experimental observations to the delocalisation of the anion in the fullerene acceptor phase.\nWe have included an extended description of the findings in Ref. 14 in the Introduction section\nas shown above.\n3. This paper https://pubs.acs.org/doi/abs/10.1021/jp5098102isclosely related to the present\nstudy and points out the entropic advantage of a large density of CT states for preventing\nrecombination.\nWe thank the reviewer for this suggestion. We now refer to this paper as a further point to\nconsider following our suggested design rules in the Discussion section on page 15:\n“Lastly, a small energy gap between the LUMO and higher lying unoccupied orbitals of the\nacceptor can further support free charge generation by increasing the entropic advantage of CT\nvs XT states, as suggested in Ref. [71]. In the present system the energy gap between LUMO and\nLUMO+1 was relatively large (1.8 eV), which is why it was not included in our simulations.\nHowever, when this gap is small, the entropic advantage is expected to further enhance the\nbeneficial effects of strong electronic coupling and a small exciton binding energy, in\ncomparison to the results reported in this study.”\nReviewer #2 (Remarks to the Author):\nThe authors present a detailed theoretical study of nanoscale donor-acceptor interfaces in a\nprototypical a6T:PDI system, using the X-SH framework to simulate non-adiabatic molecular\ndynamics. The results indicate that efficient dissociation of excitonic states into free carriers\noccurs via transiently delocalized hybrid exciton-charge transfer states. Strong electronic\ncouplings between the donor and acceptor molecules promote ‘hot exciton’ dissociation, a\nsimilar effect to increasing the dielectric constant of the heterojunction. The authors also\ndiscuss the experimental implications of their findings. Additionally, the article introduces an\nimproved X-SH methodology (accounting for off-diagonal disorder) that could benefit future\napplications to various materials. Given the significance of these findings, I recommend this\nwork for publication in Nature Communications following a revision addressing the comments\nbelow. No further review is necessary.\nWe thank the reviewer for their appreciation of our work and for recommending publication after\nrevision. Below, we provide a point-by-point response to the remaining concerns.\n1. The authors consider for their simulations idealized crystalline interfaces between donor and\nacceptor phases. Indeed, while experimentally it was found that crystalline interfaces work\nbetter than amorphous ones for charge separation efficiency, I wonder what the role of defects\nis. To what extent some typical defects, e.g. faults in the stack of both materials, affect the\nefficiency of charge transfer mechanisms. These modification can be easily incorporated and\ntested on the level of X-SH model. This discussion would be an important addition to the paper\nand point to experimentalists how defect tolerant the charge separation is, for example, for a\nhot-exciton pathway.\nWe thank the reviewer for introducing this aspect of stacking defects and its potential effect on\ncharge generation. We have examined this on the level of the X-SH model by introducing an\ninterfacial defect into the PDI acceptor phase, which disrupts its stacking with the a6T donor\nphase. This allows us to probe the effect of a local defect on exciton dissociation and the\ndynamics of the CT-states. With an interfacial defect now included in the heterojunction\nstructure, we run an extra X-SH simulation for 10ps where we select the same exciton binding\nenergy and scaling of the electronic coupling as simulation 5 in Table 2, since the CS yield here\nafter 10ps is significant (33.1%). We first refer to these new findings in the Introduction in page\n4:\n“We lastly examine the role played by structural disorder by selecting parameters that led to\nefficient charge generation in the crystalline a6T:PDI interface, and re-simulating the exciton\ndissociation process with an interfacial defect introduced to disrupt the donor/acceptor\nstacking pattern.”\n“Additionally, for a parameter combination where the hot path way dominates, the insertion of\nan interfacial defect into the a6T:PDI heterojunction significantly reduces the rate of exciton\ndissociation.”\nWe include the resulting population dynamics and analysis in a new panel in Figure 4, Figure 4d,\nand the Results section (page 12), respectively:\n1: New version of Fig. 4 in the manuscript to include the population dynamics of the a6T:PDI heterojunction with the\ndefect, panel (d) shows the comparison of electronic populations of simulation 5 with(dotted) and without(solid) a\ndefect.\n“Effect of Interface Defect. So far, we have focussed on the dynamics of free charge\nphotogeneration in a perfectly ordered heterojunction, as a crystalline morphology of\nthe donor/acceptor aggregates has proven beneficial to charge generation efficiency [4] and\nincreased the role played by delocalisation. [26–37] Here, we introduce a prototypical interfacial\ndefect into the PDI acceptor phase (see Fig. S15), to examine its potential effect on exciton\ndissociation and the dynamics of the CT-states. We focus on the dynamics of a system where\nfree charges are being mostly generated by hot exciton dissociation (simulation 5 in Table 2), as\nwe expect the hot pathway to be most affected by the introduction of a defect. The change in\npopulation dynamics upon the introduction of the defect is displayed in Fig. 4d. Exciton\ndissociation has slowed down markedly, as the disruption in a6T-PDI stacking at the interface\nhas reduced the coupling between the XT and CT states and reduced the density of hybrid XT-CT\nstates, which are necessary for exciton dissociation (see Fig. 3). This results in a significantly\nlower occupation of CT-states throughout the dynamics. Moreover, a smaller fraction of\nexcitons now dissociate via the hot pathway, whilst a more significant fraction dissociate via the\ncold pathway that is prone to non-radiative recombination. The interfacial defect stifles the hot\npathway by disrupting electron transfer between PDI molecules near the interface. This reduces\nthe density of hybrid XT-CT states with charges delocalised further away from the interface,\nwhich are instrumental to hot exciton dissociation (see Fig. 5).”\nWe then translate these findings into a fourth design rule in the Discussion section (page 15):\n“4. Structurally ordered donor/acceptor aggregates. The introduction of an interfacial defect\nthat disrupts electron transfer inthe acceptor phase near the interface resulted in the cold\nexciton dissociation pathway gaining significance with respect to the hot pathway, leaving\nelectron-hole pairs more likely to geminately recombine after exciton dissociation. 1.-4. will\ngenerally suppress the formation of kinetically trapped iCT states and increase the contribution\nof hot exciton dissociation.”\nWe also discuss the creation of the defect in a new section in Supplementary Information:\n2: Snapshot of a subset of the a6T:PDI heterojunction (Fig. S15) with an interfacial defect introduced into the PDI\nphase. Coloured molecules are those within the active region in X-SH.\n“The a6T:PDI heterojunction with the interfacial defect is shown in Fig. S15. The defect was\nmanually created by taking the final geometry of the NVT equilibration run of the crystalline\na6T:PDI heterojunction, deleting the second-most interfacial PDI molecule, and then rotating\nthe interfacial PDI molecule until the ET and XT-CT couplings involving this molecule were\nreduced by a factor of 4 or more (compared to those reported in Table 1 in the main text). The XT\ncoupling involving the interfacial PDI molecule was only reduced by ≈ 25%, since the Coulombic\ninteraction between the S1 transition densities is less sensitive to distance than the orbital\noverlap. The restraints were then adjusted to reflect the new positions of the upper/lower rings\nof the interfacial PDI molecule (see Methods for information about the application of restraints).\nThe adjusted heterojunction was then propagated in the neutral ground state for 1.1ns in the\nNVE regime, with the first 100ps allowing the system to equilibrate (in case the intermolecular\nforces were significantly perturbed by the transformation of the molecule) whilst 1000\ngeometries were sampled each ps from the remaining 1ns of simulation time. The fluctuation of\nthe temperature around 300 K was very similar to the NVE run of the fully ordered\nheterojunction.”\n2. There are lots of experimental and theoretical studies in the field pointing to the important of\nvibronic coherence that can fundamentally affect many photoinduced processes and excited\nstate dynamics. The authors use surface hopping dynamics which fundamentally lacks\nvibrational wave functions, and decoherence corrections are just fixing the non-adiabatic\nalgorithms. The authors need to include a brief discussion on this phenomenon (which to me\ncould be very important for hot exciton dissociation pathway) and how it can alter the emerging\npicture.\nWe thank the reviewer for bringing up this interesting point. As he/she correctly notes, surface\nhopping (SH) is a quantum-classical method and therefore lacks explicitly quantized vibronic\nwave functions. Nevertheless, SH represents the coupled electron–nuclear wavepacket through\nan ensemble of classical trajectories, and the nuclear (vibrational) motion is fully considered\nwithin this representation. Many nonadiabatic phenomena—such as wavepacket branching,\nelectronic population transfer, and partial delocalization of the electronic wave function—are\ndescribed with good accuracy by SH across a wide range of systems [T. R. Nelson, et al.Nat\nCommun 9, 2316 (2018).]. Decoherence is indeed a necessary algorithmic correction whose\ninclusion, along with others such as the trivial crossing correction and velocity rescaling after a\nsurface hop, is essential to recover the correct dynamics. These aspects have been thoroughly\ntested in our previous works [Carof, A., et. al.Phys. Chem. Chem. Phys.21, 26368–26386\n(2019), Giannini, S., et al. J. Phys. Chem. Lett.9, 3116–3123 (2018). Carof, A. et al.J. Chem.\nPhys.147, 214113 (2017).].\nThat said, compared to fully quantum-dynamical methods, quantum-classical dynamics may\nbe expected to capture less accurately ultrafast relaxation processes, where the quantum\ncharacter of high-frequency vibrational modes in some cases affects the wavepacket dynamics\n(as evidenced by coherent beatings in the populations [Popp, W., Polkehn, M., Binder, R. &\nBurghardt, I. J. Phys. Chem. Lett.10, 3326–3332 (2019)] with recurrences of periods shorter than\n20 fs). Nevertheless, as we have shown previously for a small model donor-acceptor system\nthrough direct comparison of our preliminary X-SH version and fully quantum dynamical\nsimulation at the level of MCTDH (using the same electronic Hamiltonian in both methods),\nsuch initial relaxation dynamics are captured with reasonable accuracy by our algorithm, even\nthough the high-frequency modes are treated classically in X-SH [Peng, W. et al. J. Phys. Chem.\nLett. 13, 7105–7112 (2022)].\nMoreover, we argue that since charge separation occurs on a longer picosecond timescale (for\nthe NFA cells probed by experiments referenced in this study [Bakulin, A.A. et al., Science, 335,\n1340–1344 (2012); Qian, D. et al., Nat. Mater., 17, 703-709, 2018]), the vibrational energy in the\nexcited state has been distributed over most vibrational modes. In systems with many degrees\nof freedom on this longer timescale, the initial vibrational coherences are strongly damped,\nand what remains is a delocalized (but incoherent) electronic wavefunction. In this regime, we argue\nthat the crucial feature of an algorithm is its ability to distribute excess energy among the\nvarious vibrational modes and to capture electronic detailed balance at equilibrium. Our\nalgorithm does so, as demonstrated in our previous work [Giannini, S. et al. Nat. Commun.13,\n2755 (2022); Giannini, S. et al.Nat. Commun.10, 3843 (2019).], yielding transport properties\nsuch as exciton diffusion and charge-carrier mobility in very good agreement with experimental\nresults.\nWe have included a paragraph on the implications of using quantum-classical dynamics in\nsimulating charge separation in the Discussion section (page 13) of the main text:\n“Besides dKMC, fully quantum-dynamical methods have been employed to study charge\nseparation in simplified donor–acceptor interface models. Previous studies have shown that the\nquantum nature of high-frequency vibrational modes can, in some cases, influence the\nwavepacket dynamics, as evidenced by coherent beatings in the populations [60]. These\ncoherent oscillations can affect relaxation processes and charge separation on ultrafast\ntimescales. From our previous comparisons with fully quantum-dynamical simulations (using\nMCTDH for the same electronic Hamiltonian), we found that such initial dynamics are captured\nwith reasonable accuracy by X-SH, even though the high-frequency modes are treated\nclassically. An important benefit of our approach over a fully quantum-mechanical treatment,\nbeyond the significant computational efficiency, is that our scheme correctly recovers the long-\ntime detailed balance, a fundamental requirement for processes occurring on extended\ntimescales where essentially all vibrational modes are in thermal equilibrium and vibronic\ncoherences are significantly damped [48].”\n3. Cold exciton dissociation. Is it correct to assume that the simulations are limited by 10ps of\nsimulations time. Within this timescale, there was no observed dissociation of the iCT states.\nThe subsequent dynamics (i.e. “charge separation or recombination of iCT state”) involves\nlonger timescales and this dynamics is not included in the kinetic models presented. Perhaps\nsome clarifying remarks would be helpful here.\nWe thank the reviewer for suggesting clarification of this point. Each X-SH simulation in this\nstudy was indeed limited to 10ps due to its computational expense, so we focussed on\ncomparing the X-SH population dynamics to the first 10ps of the kinetic model. For X-SH\nsimulations where cold exciton dissociation dominates (e.g. Table 2, simulation 1), we did not\nobserve any dissociation of the iCT states. This indicates that, for the parameter regimes herein,\nX-SH predicts that the iCT states either recombine or dissociate on timescales longer than\n10ps.\nWe have introduced an additional statement when explaining the observed cold pathway in the\nResults section (page 8), to clarify that we investigate charge generation within the first 10ps due\nto the length of the X-SH simulations:\n“The subsequent dynamics, that is charge separation or recombination of the iCT state, involves\nlonger time scales. This dynamics is not included in the present simulations as these longer\ntimescales cannot be presently accessed. The mechanism described above is referred to as\n“cold exciton dissociation” because all excitons relax to the lowest excited (“cold”) state of the\nsystem, i.e. iCT state, before they acquire enough thermal energy to separate to free charges or\nrecombine to the ground state.”\n4. Hot exciton dissociation. I would question validity of the statement “a small increase in\ndielectric constant from 3.5 to 5 (Table 2, simulation 4) results in a decrease of exciton binding\nenergy from 260 to 10 meV.” This is a huge claim and I do not thinks that just referring to Eq. 24\nprovides justice. Exciton binding energy is a delicate quantity to calculate. I would request the\nauthors to substantiate their claims. Experimental evidence is indirect (ref. 61).\nWe thank the reviewer for recommending clarification of this point. We agree that the\ncalculation of very accurate values of the exciton binding energy is a non-trivial task. First, we\nwould like to note that Eq. 24 is a widely accepted definition for exciton binding energy [Savoie,\nB.M. et al.Acc. Chem. Res.47, 3385–3394 (2014); Yufan, Z. et al.Adv. Energy Sustainability Res.\n3, 2100184 (2022)]. Second, for the purpose of this work it is sufficient to calculate the terms in\nEq. 24 with some common approximations, as outlined in the Methods section and section S3\nof the SI. The exciton binding energies quoted in the main text are valid within these\napproximations.\nWe have introduced additional comments in the Results section of the main text (page 9)\nhighlighted by the reviewer, to clarify the approximations we have made upon calculating the\nexciton binding energy of the a6T:PDI heterojunction:\n“Alternatively, a small increase in dielectric constant from 3.5 to 5 (Table 2, simulation 4) shifts\nthe niCT and CS bands down in energy to align with the XT band (Fig. 2c). Adopting the common\nexpression for exciton binding energy [53] (Eq. 24) and approximating the formation energy of\nfree charges from the iCT state by a screened Coulomb potential, we find that the increase in\ndielectric constant corresponds to a lowering of the exciton binding energy from 260 to 10 meV\n(however, we note that this is an approximation as the electron-hole interaction energy at short\nrange may not exactly follow a screened Coulomb interaction). In both cases, with an increase\nof the electronic coupling or dielectric constant, the resulting close energetic alignment also\nincreases the density of hybrid XT-CT states.”\n5. Regarding the abstract statement: “Our results highlight the importance of electronic\ncouplings, an often overlooked property, for opto-electronic charge generation.” I believe that\nthis is not true and tuning of electronic coupling has been an active subject of the discussion in\nthe community for decades across both theoretical and experimental studies. Perhaps this\nstatement and similar tone in the main body of the article need to be softened.\nWe agree the effect of electronic coupling on optoelectronic charge generation has been\ndiscussed before (e.g. by Balzer et al. in Ref. 34). Here, we intended to highlight that this\nimportant materials property has received much less attention for optimisation of materials for\norganic solar cells than it has for organic electronics. In response to the reviewers suggestion\nwe have softened the wording and changed the sentence in the abstract to:\n“Thus, our results highlight the importance of the electronic coupling and dielectric constant for\nopto-electronic charge generation.”\n6. The present approach considers singlet states only. However, I would argue that weakly\nbound triplet charge transfer states are iso-energetic with their singlet counterparts and can be\neasily form even in the absence of strong spin-orbit coupling. Their presence and participation\nin the process can add additional entropic driving force further improving charge separation.\nI would suggest to add discussion of this scenario.\nWe thank the reviewer for this interesting point. We have added a discussion on triplet formation\nin the Discussion section (page 14):\n“As mentioned earlier, we only considered singlet states in this work, disregarding any potential\nformation of low-lying triplet states. Gillett et al. provided a comprehensive experimental\nassessment of the role of triplet CT states in low-offset and NFA-based organic solar cells. [64]\nThey can be close in energy to their singlet counterparts because the exchange energy\ndecreases rapidly with increasing electron–hole separation. [65] Triplet CT states may form via a\ngeminate spin-mixing mechanism from singlet CT states on nanosecond timescales; however,\nGillett et al. found no experimental evidence for this geminate pathway in the NFA blends they\nexamined. Instead, they identified a dominant non-geminate mechanism, in which the\nrecombination of free charges generates triplet CT excitons that undergo fast back charge\ntransfer to Frenkel exciton triplets. However, the non-geminate recombination pathway\nrequiring free charge carriers happens on much longer time scales than the current X-SH\nsimulations can access. Yet, depending on the material under consideration, triplet formation\nmay occur faster, in which case it provides additional entropic driving force for free charge\ngeneration [66,67].\"\n7. Code availability statement says: “The custom codes used for this study are available from\nthe corresponding authors upon request.” It would be really helpful to clarify releases of the\ncode for public, for example, if any of the standard version of X-SH are available and what are\nthe plans for realizing the improved version. This is an important advance of this study.\nWe agree with the reviewer that a public release of the X-SH code is beneficial for the\nadvancement of this study. The code utilised in the X-SH simulations presented herein is\nincorporated into the latest version of our in-house modified CP2K-V4.1 package. This is a\ndevelopment version, and the latest fully tested version of the code is publicly available on\nhttps://github.com/blumberger/flavoured-cptk-public-copy.git. We have also amended our\nCode Availability statement accordingly:\n“The latest tested version of the custom code used for this study is publicly available on\nhttps://github.com/blumberger/flavoured-cptk-public-copy.git.”\nReviewer #2 (Remarks on code availability):\nI do not think that X-SH code used for this study was released, see my point #7\nPlease see the answer above.\nAdditional changes:\nWe have slightly amended two sentences to make better use of supplementary figures in\nsupporting the conclusions drawn in the main text (page 11):\n“This novel state is a superposition of delocalised XT, iCT and niCT states. In such states, the\nexcitons and charges are delocalised to a greater extent (see Fig. S6) and can be centred further\nfrom the donor/acceptor interface (Fig. S8).”\n“A remarkable feature of the hot exciton dissociation mechanism is that a significant proportion\nof excitons no longer dissociate at the interface, as traditional models of excitonic solar cells\nwould assume[2], but several lattice spacings / nanometres away from it, owing to the\ndelocalised nature of hybrid XT-CT states in this regime (see Figs. S4 and S10, using the center\nof the wavefunction as a distance measure).”"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have thoughtfully revised their paper, which was already excellent. I enthusiastically recommend acceptance.", "reviewer_2": "The authors have thoroughly addressed all reviewer's comments, with clear attention to detail and a deep understanding of the subject matter. Their responses are thoughtful, comprehensive, and convincing, and the corresponding revisions in the manuscript significantly strengthen the work. Importantly, the authors have also released the relevant code, further supporting the transparency and reproducibility of their results. Given these efforts, I strongly recommend publication of the manuscript in Nature Communications in its present form."}, "author_response": ""}]

Decision letter