Even a little delocalization produces large kinetic enhancements of charge-separation efficiency in organic photovoltaics.

In organic photovoltaics, charges can separate efficiently even if their Coulomb attraction is an order of magnitude greater than the available thermal energy. Delocalization has been suggested to explain this fact, because it could increase the initial separation of charges in the charge-transfer (CT) state, reducing their attraction. However, understanding the mechanism requires a kinetic model of delocalized charge separation, which has proven difficult because it involves tracking the correlated quantum-mechanical motion of the electron and the hole in large simulation boxes required for disordered materials. Here, we report the first three-dimensional simulations of charge-separation dynamics in the presence of disorder, delocalization, and polaron formation, finding that even slight delocalization, across less than two molecules, can substantially enhance the charge-separation efficiency, even starting with thermalized CT states. Delocalization does not enhance efficiency by reducing the Coulomb attraction; instead, the enhancement is a kinetic effect produced by the increased overlap of electronic states.

INTRODUCTION

The precise mechanism of how charges in many organic photovoltaics (OPVs) separate with near-perfect efficiency remains unclear, especially considering that their Coulomb attraction can be more than an order of magnitude larger than the available thermal energy (, ). OPVs typically contain a heterojunction of electron donor and electron acceptor materials. Upon excitation, an exciton is created in either material and diffuses to the donor-acceptor interface, where charge transfer can form an interfacial charge-transfer (CT) state (left of Fig. 1A). For the photogeneration process to continue, the charges must escape their Coulomb attraction and separate (right of Fig. 1A). In organic semiconductors with the typical dielectric constants of 3 to 4, charges separated by 1 nm across the interface experience a Coulomb attraction of 360 to 480 meV, more than an order of magnitude larger than the thermal energy kBT = 25 meV.

Fig. 1.

Mechanism of delocalization-enhanced charge separation.

(A) Charge separation is typically modeled (purple insets) by assuming that the charges are localized onto individual molecules or sites (spheres) by disorder (different shades). When both charges are localized at the donor-acceptor interface, they have a small electron-hole separation reh and a large Coulomb attraction U(r). It has been proposed that delocalization of charges across many molecules (green insets) facilitates their separation by increasing their initial separation and decreasing their Coulomb attraction. (B) While delocalization does decrease the Coulomb binding energy (∆U), it increases the overall binding energy (∆E), meaning that a reduction in the CT state binding energy is not the cause of delocalization enhancements. Instead, delocalization increases the overlaps between electronic states, allowing delocalized charges (green clouds) to hop further and faster than localized ones (purple dots).

The efficient separation is commonly attributed to charge delocalization, an argument supported by the experimental observation of efficient and fast separation via delocalized CT states (–). On this view, the delocalization increases the initial separation of charges in the CT state, reducing the Coulomb attraction and the energetic barrier that needs to be overcome (Fig. 1A). However, the concept of a kinetic barrier to charge separation is problematic, because separation is not a one-step process from a CT state to a separated state but a multiple step process through many states, where each transition has its own energetic barrier (Fig. 1B). Furthermore, the barrier can completely disappear in free energy when entropy and disorder are considered () and actually increase when delocalization is included (). Therefore, any charge-separation efficiency improvement caused by delocalization must come from nonequilibrium kinetic effects, not purely energetic considerations (, ).

Many kinetic models of charge separation have included charge delocalization, with most of them finding an increased efficiency when delocalization is considered (–). These approaches have ranged from quantum-mechanical descriptions of delocalization in disordered materials (–) to phenomenological treatments that include delocalization in an effective way (–). However, different approaches have used approximations that limit their range of applicability in three important ways. First, approaches that do not include static disorder are applicable to organic crystals but not disordered organic semiconductors. Second, approaches that treat delocalization in an effective way tend to not include the complete quantum-mechanical description required for accurately predicting the effects and extent of delocalization. Last, approaches that do not adequately describe the coupling of the charges to their environment can fail to describe the formation of polarons, which can localize the states significantly (), meaning that those approaches can overestimate delocalization and efficiency enhancements. Therefore, an accurate kinetic model of charge separation in OPVs must include these three important ingredients: disorder, a quantum-mechanical treatment of delocalization, and polaron formation.

The most complete kinetic models of charge separation have included all three important ingredients (–). However, they remain computationally expensive, limiting their application to small systems, low dimensions, or short times. Tamura and Burghardt used an atomistic approach, solving the time-dependent Hamiltonian using multiconfigurational time-dependent Hartree (, , ) to predict the ultrafast separation of an interfacial exciton over about 100 fs, attributing the efficient separation to vibronically hot CT states and to charge delocalization reducing the Coulomb attraction. While atomistic approaches do not require adjustable parameters, they are expensive and therefore limited to relatively few states and short times. Bittner and Silva () and Janković and Vukmirović () both parameterized model Hamiltonians, which allows for longer simulations on more states, as fewer degrees of freedom are tracked. While their approaches are computationally cheaper, they are still limited to two- and one-dimensional simulations, respectively. Both approaches include the formation of polarons using the polaron transformation. Bittner and Silva () studied the separation of charges in a two-dimensional heterojunction using a time-dependent Schrödinger equation. They observe fast separation into free polarons, due to charges tunneling from the initial exciton in under 35 fs (rather than proceeding through CT states), due to resonance between excitons and free polarons. However, the environment was treated as a bath of extended, shared phonons (which is common for crystalline systems), rather than as local molecular vibrations, which is more appropriate for disordered molecular systems (). Janković and Vukmirović () modeled charge separation in a one-dimensional system over long times to study the competition of delocalization, disorder, and polaron effects. They found that separation proceeds slowly through thermalized CT states; after excitons form CT states in 1 to 10 ps, the latter separate in about 1 ns. However, rates were calculated using modified Redfield theory, which requires small off-diagonal system-bath couplings, an assumption that is not met in many materials, including some organic semiconductors, where strong system-bath coupling alone can localize the electronic states ().

The best kinetic models have remained computationally expensive for two main reasons. First, correctly describing delocalization requires a quantum-mechanical treatment, which is difficult in disordered materials, where periodic boundary conditions must be replaced with large simulation boxes. Second, charge separation is a two-body problem involving the correlated motion of an electron and a hole, meaning that the Hilbert space size is roughly the square of a single-body calculation. As a result, a quantum-mechanical treatment has so far proved intractable in three dimensions (). To study charge separation over long times, in large systems, and in three dimensions, we require a more efficient approach that still meets all requirements.

A logical starting point for understanding delocalization in charge separation is extending kinetic models that describe the transport of partially delocalized single particles (–). However, even this task is difficult when disorder, delocalization, and polaron formation are required (). The best recent approaches to tackle the single-particle problem are all effective Hamiltonian methods and include fragment orbital-based surface hopping (–), where the Hamiltonians are parameterized using an atomistic approach requiring no adjustable parameters, as well as adaptive hierarchy of pure states equations (), which emphasizes the nonperturbative treatment of the system-environment couplings using hierarchy equations of motion. The computational cost of these approaches has limited their application to short times in two-dimensional systems or large one-dimensional systems.

Recently, we introduced delocalized kinetic Monte Carlo (dKMC), the first three-dimensional model of charge and exciton transport that includes disorder, delocalization, and polaron formation (). dKMC is a computationally improved version of secular polaron-transformed Redfield theory (sPTRE) (), which is the stationary and secular limit of multichromophoric coherent resonance energy transfer (–), and one example of a second-order polaronic master equation (–). Because sPTRE is in the polaron frame, which changes as a function of system-bath coupling, it can describe polaron transport in the intermediate regime while also reproducing the well-known hopping and band conduction extremes (). dKMC shows that transport in moderately disordered materials is that of charges hopping between partially delocalized electronic states and that even a small amount of delocalization can increase carrier mobilities markedly ().

Here, we extend dKMC to make the charge-separation problem computationally accessible, allowing the first three-dimensional simulation of the dynamics and efficiency of charge separation in the presence of disorder, delocalization, and polaron formation. We use it to show that small amounts of delocalization can produce large efficiency enhancements, even for thermalized CT states. Contrary to the common hypothesis, these delocalization enhancements are not a consequence of a reduction in the initial Coulomb binding energy. Rather, delocalization actually increases the total binding energy, and the efficiency enhancements are a kinetic effect caused by greater overlaps between electronic states, which allow charges to move further and faster (Fig. 1B). All of the approximations used in this work are conservative, chosen to underestimate the extent and role of delocalization; in the end, we discuss ways to relax the approximations, which could lead to an even greater role for delocalization.

RESULTS

Model

Our approach is based on dKMC (). Here, we summarize the changes used to extend dKMC to the two-body charge-separation problem, which are detailed in Methods.

The charge-separation problem is described by the Hamiltonianwhose components describe the system HS, the bath HB, and the interaction between them HSB.

We model the system using a tight-binding model of a d-dimensional cubic lattice containing N/2 donor sites next to N/2 acceptor sites, each representing a molecule or part of a molecule (Fig. 2A). Each site is assigned a highest occupied molecular orbital (HOMO) and a lowest unoccupied molecular orbital (LUMO) energy, so that an electron occupying that site will have the LUMO energy and a hole will have the HOMO energy. These energies are assumed to be disordered to model the different local environments around different molecules, which arise from static variations in the spacing and orientation of molecules. The HOMO and LUMO energies are drawn randomly from Gaussian distributions and , where the energetic disorder σ is assumed to be equal for the HOMOs and the LUMOs ().

Fig. 2.

Components of the dKMC model.

(A) A heterojunction is modeled as a lattice of sites with disordered energies (different shades), representing acceptor (orange) and donor (blue) molecules. Each site is coupled to an environment (motion lines) and to its neighbors. (B) Delocalization of the electronic states is found by diagonalizing the system’s Hamiltonian, with the eigenstates representing the simultaneous position of the electron (in the acceptor, blue) and the hole (in the donor, orange). (C) Polaron formation further localizes the states. Inverse participation ratios (IPRs) are shown for J = 30 meV and σ = 150 meV.

In the two-body problem, the Hilbert space consists of ordered pairs of sites, where ∣m, n⟩ represents an electron on donor site m and the hole on acceptor site n. The energy of this pair is , where is the LUMO energy of site m, is the HOMO energy of site n, and U(r) is the Coulomb potential of charges separated by a distance rwhere e is the elementary charge, ε0 is the vacuum permittivity, and ε is the dielectric constant (here, taken to be ε = 3.5).

A pair of sites is electronically coupled to other pairs of sites that can be obtained by moving either the electron or the hole, but not both, giving the system HamiltonianIn general, the couplings J and J can be disordered or long range; however, we assume only constant nearest-neighbor couplings with strength J.

To model the bath, every site is assumed to have an identical, independent bath of harmonic oscillators representing molecular vibrations, which is the usual assumption for disordered molecular materials (, ). Then, the bath Hamiltonian is given bywhere ω is the frequency of the kth bath mode at the nth site, and and b are the corresponding creation and annihilation operators.

The system-bath interaction is described by a coupling of every site to its bath modes. We assume a linear coupling of strength g between site n and bath mode k, so that the interaction Hamiltonian is

The formation of polarons [quasi-particles containing a charge and the distortion it causes to the bath (, )] is described by applying the polaron transformation to H, which displaces the bath modes using the state-dependent displacement operator ()Using the polaron transformation has two benefits. First, it reduces the system-bath coupling by absorbing most of the interaction into the polaron itself, allowing the residual interactions to be treated perturbatively (). Second, it reduces the electronic couplings within the system, giving smaller polarons () and reducing the complexity of dKMC calculations ().

The results of applying the polaron transformation to H are given in Methods. The polaron-transformed system Hamiltonian is then diagonalized to find the joint polaron states of the electron and the hole. The extent of the delocalization of polaron state ν can be quantified using the inverse participation ratio (IPR)which describes roughly how many pairs of sites (m, n) the state ν is delocalized over. The IPRs calculated using Eq. 7 depend on all the parameters of the model; in particular, IPRs increase with electronic coupling J and decrease with disorder σ and system-bath coupling g. As shown in Fig. 3, the Coulomb interaction stabilizes (Fig. 3A) and localizes (Fig. 3, B and C) states in which the electron and the hole are close together, such as CT states.

Fig. 3.

Properties of the polaronic states.

In (A) to (C), the dots represent the polaronic states of a two-dimensional (2D) heterojunction with J = 75 meV and σ = 150 meV. (A) The Coulomb interaction stabilizes states with small electron-hole separations (the orange line is the Coulomb potential and the dashed lines are ±2σ on either side). At any reh, including for CT states, there are states with a wide range of energies. (B) States with small reh (and, therefore, large Coulomb attractions) are more localized. (C) Low-energy states are more localized, while the more delocalized states lie closer to the middle of the density of states. (D to F) The average separation, IPR, and energy of the three kinds of initial CT states, as a function of the electronic coupling J. Overlap CT states are chosen based on each state’s overlap with interfacial CT site pairs, producing states that are the most separated, delocalized, and highest in energy. Random CT states are chosen uniformly from the states with the smallest reh. Thermalized CT states are chosen from the states with the smallest reh in proportion to their Boltzmann factor, producing initial states that are the least separated, least delocalized, and lowest in energy. The error bars are the SEMs.

We begin dKMC simulations with the charges in a CT state. Because of disorder, there are many CT states of various energies and multiple plausible ways they could be chosen, as reflected in the debate about whether separation proceeds through low-energy, thermalized CT states or high-energy, hot CT states. To examine various possibilities, we choose the initial CT state in one of three ways. Random CT states are chosen uniformly at random from the NCT polaron states with the smallest electron-hole separations, where NCT is set equal to the number of interfacial pairs on the lattice. Thermalized CT states are chosen from the same NCT most closely bound polaron states but in proportion to their Boltzmann factors. Last, overlap CT states are chosen from all the polaron states, in proportion to their overlap with interfacial pairs, ∣∑(⟨m, n∣ν〉∣2. This choice mimics states prepared by charge transfer from an interfacial exciton, which would be most strongly coupled to interfacial pairs. Among the three initial conditions, thermalized CT states are, on average, lowest in energy, most bound, and least delocalized, while overlap CT states are at the opposite extreme (Fig. 3, D to F).

Starting with any of these initial CT states, dKMC is a form of kinetic Monte Carlo (KMC) that uses polaron-frame Redfield rates, computed in Methods, to evolve trajectories through the polaron states. As in regular KMC [which instead uses Marcus () or Miller-Abrahams () rates], the next state is chosen with a probability proportional to the rate of transfer to the state (, , , , ). The hopping continues until a terminating condition is reached. If the polarons are separated by more than a chosen separation distance rsep, which we set to 5 nm, the charge separation is considered successful. By contrast, the charge separation is considered failed if the polarons recombine or if the number of hops exceeds the limit nhops, which we set to 2000. The latter assumption avoids infinite loops, such as hopping between two low-lying states, where the polarons would probably eventually recombine. Last, the internal quantum efficiency (IQE) is calculated as the fraction of trajectories where the polarons separate, averaged over many realizations of disorder.

Delocalization increases charge-separation efficiency

Figure 4 shows that delocalization increases the efficiency of charge separation by comparing the IQE calculated with localized polarons, using a standard KMC approach, to that calculated with partially delocalized polarons, using dKMC. When the electronic coupling J is low, the electronic states are localized, meaning that KMC and dKMC agree. As J increases, the states become more delocalized and the delocalization enhancement predicted by dKMC increases regardless of what initial CT states are chosen.

Fig. 4.

Delocalization and dimensionality increase charge-separation efficiency.

(A) IQEs of charge separation, for a system with σ = 150 meV and varying electronic coupling J, modeled using regular KMC (blue) and dKMC (orange). When the states are localized (small J), dKMC and KMC agree, but as J increases, delocalization significantly enhances the dKMC efficiency over the classical KMC hopping. For all J, overlap CT states separate slightly more efficiently than random CT states and significantly more efficiently than thermalized ones. (B) IQEs of charge separation starting in random CT states, with σ = 150 meV and varying J, modeled by KMC and dKMC in each dimension. Delocalization enhances IQE in all dimensions but more strongly in higher dimensions. The 3D dKMC line stops at a relatively small value of J because of computational cost. The error bars in both panels are the standard errors of the mean (SEMs).

Figure 4A shows significantly different efficiencies for thermalized, random, or overlap initial CT states. In both dKMC and KMC, the thermalized CT states separate with a lower efficiency than random CT states, suggesting that, for the parameters studied here, a randomly excited CT state is able to separate before it thermalizes. This result supports the idea that hot CT states can improve charge separation in OPVs, as has been observed experimentally (–). However, our results are also consistent with observations of highly efficient separation from thermalized CT states (–), because dKMC can reproduce high IQEs from thermalized states, especially if they are assisted by delocalization.

These large enhancements can be caused by small amounts of delocalization. Figure 5A shows the IQE of random and thermalized CT states as a function of their IPR, illustrating how quickly the efficiency grows with delocalization. At the highest coupling, J = 75 meV, thermalized and random CT states have IPRs of only 1.5 and 2.6, but these are enough to produce fivefold and twofold IQE enhancements, respectively. Therefore, just as how small amounts of delocalization can markedly improve polaron mobilities (), they can also facilitate the separation of charges experiencing a significant Coulomb attraction, a process that is underestimated by standard KMC.

Fig. 5.

Kinetic origin of delocalization enhancements.

(A) The charge-separation efficiency (IQE) increases with the delocalization of the initial CT states (IPR, shown for J from 0 to 75 meV). Results are for a 2D heterojunction with σ = 150 meV. (B) The efficiency increase with delocalization is typically attributed to larger initial electron-hole separation reh and the resulting reduction in the Coulomb binding energy ∆U. Although ∆U does decrease with reh, the total binding energy of the initial state (∆E) actually increases with reh due to delocalization stabilization. (C) Therefore, the efficiency increases with the total binding energy of the initial state, contrary to the common hypothesis that delocalization enhancements are simply a consequence of a reduction in binding energy.

Kinetic effects are more important than the reduction in Coulomb attraction

Delocalization has usually been argued to enhance the IQE by increasing the initial electron-hole separation in the CT state, thus reducing the Coulomb attraction and binding energy. However, this standard argument is incorrect; instead, delocalization enhancements are a kinetic effect chiefly caused by the increased overlap of electronic states (Fig. 5). It is true that the Coulomb binding energy ΔU = U∞ − UCT = −UCT decreases with initial reh, as shown in Fig. 5B. However, the Coulomb binding energy is not the sole energetic penalty of charge separation, and Fig. 5B also shows the total binding energy ΔE as a function of the reh of the initial CT states. It is calculated as ΔE = E∞ − ECT, where E∞ is the equilibrium energy of two separated polarons, each calculated independently as the thermal expectation value of the energy of a single polaron in a box of size . As shown in Fig. 5B, ΔE increases with increasing initial delocalization and reh.

The increase in total binding energy due to delocalization was predicted by Gluchowski et al. (), who attributed it to two mechanisms. First, level repulsion between coupled pairs of sites will generally stabilize the lower-energy state, further increasing the binding energy of already most tightly bound states. Second, the lowest-lying CT states generally remain localized as J increases because they are unlikely to have near-resonant neighbors. By contrast, higher-lying CT states become more delocalized at high J, increasing their reh until they no longer enter into the averaging for ΔU and ΔE. Overall, at higher J, ΔE increases because of the greater contribution from low-lying localized traps, which are themselves stabilized by level repulsion.

Overall, Fig. 5C shows that the IQE increases with the binding energy, showing that the delocalization does not enhance efficiency by reducing the energetic barrier to charge separation. Instead, delocalization enhancements must be a purely kinetic (as opposed to energetic) effect, caused by the increased overlaps between the partially delocalized electronic states, which help the polarons move apart faster and further in fewer hops.

Delocalization enhancements are greater in higher dimensions

Previously, we found that delocalization enhances the transport of polarons more in higher dimensions (). In particular, polarons are more delocalized in higher dimensions because of the increased number of neighboring sites, which increase the likelihood of having a neighboring site with similar energy to allow delocalization.

Higher dimensions are also important for accurately modeling charge separation. Figure 4B shows that IQEs computed using both regular KMC and dKMC are higher in higher dimensions. In both cases, higher dimensions increase the number of nearest neighbors that polarons can move to, increasing the likelihood of a fast separation pathway. This effect becomes more pronounced when delocalization is described using dKMC. When the polarons are delocalized, they can move further in one hop, and the number of possible destinations increases with dimension faster than in localized hopping, especially considering the greater delocalization lengths in higher dimensions ().

Delocalization enhancements increase with moderate disorder

dKMC’s computational savings allow parameter scans that independently probe the roles of microscopic parameters. Figure 6 shows the IQE as a function of both the electronic coupling J, which increases delocalization, and energetic disorder σ, which localizes states. In the limit of small J, when the electronic states are localized onto individual sites, dKMC and KMC agree. As J increases, causing delocalization, the efficiency of separation is always greater when delocalization is included using dKMC.

Fig. 6.

Delocalization enhancements increase with moderate disorder.

Parameter scan of IQEs of charge separation in 2D for (A) random and (B) thermalized CT states, as a function of electronic coupling J and disorder σ, using KMC (blue) and dKMC (orange). dKMC agrees with KMC in the low-J limit (when the states are localized) and always predicts IQEs greater than KMC as J increases. For both random and thermalized CT states, higher J always leads to higher IQE; the behavior as a function of disorder is more complicated, although modest amounts of disorder generally increase the IQE.

Figure 6 reveals the subtle effect of disorder on IQE. It is known that a modest amount of disorder can be beneficial for the separation of localized charges (, , ) because it provides energetically downhill pathways for the separation of many initial CT states. Figure 6 shows that delocalization extends the regime over which disorder is beneficial for charge separation, because it increases the overlap between electronic states allowing for polarons to move further and faster. For extreme disorders, the delocalization enhancement eventually decreases as the states become completely localized, giving agreement between dKMC and KMC.

DISCUSSION

Computational limitations

Despite the savings enabled by dKMC, it remains limited by computational cost. Especially in three dimensions, dKMC is limited to smaller values of the electronic coupling (Fig. 4B) because the polaron states become more delocalized at higher J, requiring the diagonalization of large Hamiltonians. This problem is particularly acute in the two-particle Hilbert space, where Hamiltonian diagonalization scales as in three dimensions. However, the results in Fig. 4B already reveal significant enhancements at low J, especially compared to lower dimensions, and we expect that the enhancement would continue to increase at higher couplings. Future computational advances may extend the results to higher J.

Possible extensions

In the future, dKMC could be extended to address other important questions in the physics of disordered materials. We envisage both extensions that relax some of our assumptions and applications to previously unexplored physical problems.

Most assumptions in the current version of dKMC could be relaxed. Doing so would likely lead to an increase in the predicted role of delocalization, because all the assumptions in this work are conservative, chosen to underestimate the extent and effects of delocalization.

The simplest extension would be the inclusion of non–nearest-neighbor couplings J and J in Eq. 3. Doing so would not require any alterations to the dKMC algorithm but would result in more delocalized polaron states, assuming the nearest-neighbor couplings are kept constant.

More generally, while we studied parameter ranges where the underlying sPTRE master equation is valid, more general master equations could be used in other cases. In particular, sPTRE can be inaccurate for systems weakly coupled to slow baths, when the secular and Markov approximations fail (, , , ). In those cases, electronic dynamics can be faster than bath relaxation, meaning that the bath modes are not fully relaxed, as assumed in Eq. 6. dKMC could be applied to those systems by replacing the fully displaced polaron transformation with its variational counterpart (, –), which would also allow ohmic or sub-ohmic baths to be treated.

Similarly, the secular approximation inherent in sPTRE could be relaxed. The secular approximation leads to a neglect of coherences between polaronic states, which reduces the number of elements of the density matrix being tracked from O(N2) to O(N). In the regime we considered, this is a good approximation (), because coherences between charge states are rarely induced and almost never survive long enough to affect long-time charge transport. Nevertheless, dKMC could be expanded to include coherences if desired and if the computational resources allow it.

A final pair of assumptions in our treatment is that both the static disorder (Eq. 3) and the system-bath coupling (Eq. 5) are diagonal in the site basis. More generally, the Hamiltonian could be modified to include both off-diagonal disorder and off-diagonal system-bath couplings. The latter are particularly important in organic crystals, where the fluctuations of electronic couplings are an additional source of localization of the electronic states (, ). Of the two changes, off-diagonal disorder could be implemented without any changes to the dKMC algorithm by replacing the constant J with a site-dependent one. By contrast, off-diagonal system-bath coupling would require modifications to the equations of motion. Polaron theories generally rely on the ability of the polaron transformation to eliminate diagonal system-bath couplings; therefore, they are unable to fully eliminate off-diagonal ones, and the remainder would have to be included in the perturbative treatment.

Two-body dKMC could also be extended beyond the separation of CT states. In particular, adding excitonic states to the simulation would give a quantum treatment of exciton dissociation both in the bulk and at interfaces. It could also describe exciton–to–CT state transfer and predict the nature of CT states that mediate charge separation, obviating the need for the comparison between the random, thermalized, and overlap CT states introduced above. Doing so may explain and unite the seemingly disparate experimental observations of efficient separation requiring hot CT states and those occurring from thermalized CT states. The inclusion of excitons would, however, pose computational challenges, because it would require replacing nearest-neighbor electronic couplings with long-range excitonic ones. A process similar to exciton dissociation is singlet fission, and we expect that an excitonic dKMC code could be adapted to treat that problem as well, replacing the separated polarons with separated triplets.

Last, we anticipate that it will be possible to use dKMC to parameterize drift-diffusion simulations by computing the necessary dissociation and recombination rates. Doing so would allow a multiscale simulation of a broad range of processes in organic semiconductors in a way that takes quantum effects into account.

In conclusion, our results explain how even small amounts of delocalization substantially assist polarons in escaping their large Coulomb attraction in OPVs. This is largely a kinetic effect, caused by increased overlap of electronic states, as opposed to the common hypothesis of a reduction in Coulomb binding. For instance, the initial delocalization across less than two molecules can suffice to make the separation of thermalized CT states five times more efficient. Furthermore, higher-dimensional effects, including greater delocalization, are particularly important to capturing the physics of charge separation. These results have been made possible by dKMC, a robust and computationally efficient technique for modeling charge-separation dynamics that includes all of the important features: three dimensions, delocalization, disorder, and polaron formation. Overall, dKMC opens new avenues in the exploration of organic semiconductors by allowing a more comprehensive theoretical description than has been possible.

METHODS

Here, we extend dKMC from the one-particle mobility calculation () to the two-particle charge-separation problem by applying the polaron transformation to the Hamiltonian, establishing the sPTRE master equation for the two-particle picture and then establishing the dKMC procedure for charge separation.

Polaron transformation

Applying the polaron transformation to the total Hamiltonian (Eq. 1) yields

The polaron-transformed system Hamiltonian iswhere , and the electronic couplings are renormalized bywhere β = 1/kBT and we take T = 300 K. The bath Hamiltonian is unchanged, , while the polaron-transformed interaction Hamiltonian iswhere

We apply two standard approximations to reduce the computational cost of summing over many bath modes: All sites couple to their baths with the same strength, g = g, and the spectral density is replaced with a continuous function. In principle, one could choose structured spectral densities for specific organic molecules; here, we choose the widely used super-ohmic spectral density (, , , , ), where λ = 100 meV is the reorganization energy and ωc = 62 meV () is the cutoff frequency.

The two approximations simplify the renormalization factors in Eq. 10 to

Because κ < 1, it always reduces the electronic coupling in Eq. 9, meaning that the polaron transformation reduces the delocalization of the electronic states ().

Secular Redfield theory

We use the sPTRE master equation (). As the polaron transformation reduces the system-bath coupling, we apply the second-order perturbative Redfield theory to (). Following the secular approximation, which is accurate for most disordered materials of interest (, ), we obtain the sPTRE master equationwhich describes the evolution of the populations of polaron states, found by diagonalizing . This evolution is determined by the secular Redfield tensor Rννwhere the damping rates arewith ων = Eν − Eμ andThe bath correlation function, where the hats denote the interaction picture, is given by ()where λ = δ − δ + δ − δ and

dKMC

In principle, the full evolution of the density matrix of both the electron and hole polarons could be tracked using the sPTRE master equation. However, doing so is expensive for three reasons. First, finding the polaron states requires diagonalizing , which scales as O(N6). Second, calculating all components of the Redfield tensor (Eq. 15) requires rates between all pairs of states, of which there are O(N4). Last, calculating each of these rates requires computing the damping rates (Eq. 16), each of which includes a sum over O(N5) terms. Therefore, propagating the full sPTRE master equation would scale roughly as O(N6) + O(N9).

To reduce the complexity, we use the four approximations developed for the original dKMC method (). Summarized in Fig. 7 and Algorithm 1, they are as follows:

Fig. 7.

Computational components of dKMC for charge separation.

(A) The full sPTRE master equation tracks the time-dependent evolution of the full charge density of the electron and hole through all polaron states. To avoid the computational cost of doing so, we apply four approximations (B to E). (B) KMC: We average many trajectories, formed probabilistically from sequential hops. (C) Hopping radius rhop: We only calculate rates to polaron states where the center (black dots) of the electron and hole is close enough. (D) Overlap radius rove: When calculating transfer rates between polaron states, we only consider pairs of sites (lattice points) that are close to both states. (E) Diagonalizing on the fly: Rather than diagonalizing the full Hamiltonian, we diagonalize a subsystem surrounding each polaron. Once either polaron moves too close to the edge of its box, we diagonalize a new subsystem.

dKMC maps the master equation onto KMC (Fig. 7B). Instead of propagating the full density matrix, it tracks, and averages over, ntraj stochastic trajectories. The trajectories are found by hopping from one state to another, with the target chosen probabilistically in proportion to the corresponding Redfield rate. This procedure is continued until the polarons recombine, become separated by a distance rsep, or hop more than nhops times. The trajectory approach means that only outgoing rates need to be calculated before each hop, rather than between all pairs of states, reducing the number of rate to be calculated from O(N4) to O(N2nhopsntraj). The IQE is the percentage of the trajectories where polarons separate, averaged over niters realizations of disorder.

Instead of calculating all outgoing rates from the current state, dKMC only calculates rates to states where the combined distance that the electron and holes hop is less than a hopping radius rhop (Fig. 7C). This approximation is justified by the small spatial overlap, and therefore Redfield rate, between states where either the electron or the hole (or both) hops far away. To calculate the combined hopping distance, we consider the positions of the electron and hole in state ν to be the expectation values and and only calculate rates from state ν to states ν′ if . Doing so reduces the number of rates calculated at each hop from O(N2) to . We use rhop values benchmarked for single-particle dynamics (). The approximation is controlled by increasing or decreasing rhop to achieve the target accuracy.

When calculating each rate, instead of summing over N5 quintuples (m ∈ D and p, q, p′, q′ ∈ A, or n ∈ A and p, q, p′, q′ ∈ D) in Eq. 16, we only sum over a truncated set (Fig. 7D), because Anderson localization of the wave functions makes overlaps between distant states small. To do so, we neglect all terms in Eq. 16 that contain overlaps of the form ⟨α ∣ x, y⟩ such that , where α is one of μ, ν, μ′, or ν′, and x and y are any of m, n, p, q, p′, or q′. This reduces the number of terms in each rate calculation from O(N5) to . The accuracy is controllable by adjusting rove; we use rove values calculated in ().

Instead of diagonalizing the entire Hamiltonian to generate polaron states, we diagonalize on the fly smaller Hamiltonians encompassing pairs of sites within boxes of size surrounding both polarons (Fig. 7E). The polarons can move within their boxes until either one gets too close to the edge of its box, when the boxes are moved and rediagonalized. As the computational bottleneck is often diagonalizing large Hamiltonians, we make the boxes as small as possible, choosing Nbox = (rhop + rove), and rediagonalize when either polaron is within Nbox/2 of the edge of its box. This reduces the diagonalization cost from O(N6) to for nrediag rediagonalizations.

Overall, dKMC reduces the computational complexity from O(N6) + O(N9) to .

Given parameters N, d, σ, J, λ, ωc, T, Rrecomb, rsep, nhops, niter, and ntraj:

1. Calculate rhop and rove as described in ().

2. For niter realizations of disorder:

a. Generate N lattice of sites, with N/2 random donor HOMO energies and N/2 random acceptor LUMO energies.

b. Create a polaron-transformed containing pairs of sites within a box of size at the center of the lattice. Diagonalize to find the polaron states, their energies, and the positions of electrons and holes in every state.

c. For ntraj trajectories:

i. Set nsep ← 0 and choose an initial state ν in any of the ways described in the text.

ii. For nhops hops:

A. Create a list L of all states ν′ such that .

B. Calculate Rνν for all ν′ ∈ L using Eq. 15, neglecting all terms in Eq. 16 that contain overlaps of the form ⟨α ∣ x, y⟩ such that .

C. Calculate using Eq. 21 and append g to L.

D. Set for all ν′ ∈ L and set T ← ∑νSν.

E. Find ν′ such that Sν < uT < Sν, for uniform random number u ∈ (0,1], and update ν ← ν′.

F. If ν = g, exit the for loop.

G. If , set nsep ← nsep + 1 and exit the for loop.

H. If or is within Nbox/2 of the edge of the current boxes, diagonalize a new containing pairs of sites within two boxes of size centered at and .

d. Calculate IQE = nsep/ntraj.

3. Calculate mean IQE by averaging all IQEs.

Recombination

In standard KMC, CT state recombination is often described with a constant rate Rrecomb, occurring when the electron and hole are on adjacent sites across the donor-acceptor interface. For dKMC, we calculate the recombination rate of CT state ν using Fermi’s golden rule,where the sum goes over interfacial site pairs and ρrecomb is the density of states. Assuming the interfacial pairs are coupled to the ground state with constant coupling ⟨m, n∣H∣g⟩ = Jrecomb, the rate becomeswhere Rrecomb = 2π∣Jrecomb∣2ρrecomb. Recombination therefore occurs at the ordinary Monte Carlo rate modified by the square of the sum of the CT state’s amplitudes on all interfacial site pairs. This result agrees with generalized Marcus theory for transfer between weakly coupled delocalized states () and with previous recombination studies (). In calculations, we used Rrecomb = 1010s−1.