Delocalised kinetic Monte Carlo for simulating delocalisation-enhanced charge and exciton transport in disordered materials

Abstract

Charge transport is well understood in both highly ordered materials (band conduction) or highly disordered ones (hopping conduction). In moderately disordered materials—including many organic semiconductors—the approximations valid in either extreme break down, making it difficult to accurately model the conduction. In particular, describing wavefunction delocalisation requires a quantum treatment, which is difficult in disordered materials that lack periodicity. Here, we present the first three-dimensional model of partially delocalised charge and exciton transport in materials in the intermediate disorder regime. Our approach is based on polaron-transformed Redfield theory, but overcomes several computational roadblocks by mapping the quantum-mechanical techniques onto kinetic Monte Carlo. Our theory, delocalised kinetic Monte Carlo (dKMC), shows that the fundamental physics of transport in moderately disordered materials is that of charges hopping between partially delocalised electronic states. Our results reveal why standard kinetic Monte Carlo can dramatically underestimate mobilities even in disordered organic semiconductors, where even a little delocalisation can substantially enhance mobilities, as well as showing that three-dimensional calculations capture important delocalisation effects neglected in lower-dimensional approximations.

---

Charge and exciton transport is fundamental to materials science, particularly in applications for energy storage and conversion, including photovoltaics, batteries, light harvesting systems, lighting and electrocatalysts. However, many next-generation materials that promise significant functional improvements are disordered and noisy, making them difficult to treat mathematically and improve computationally. The clearest example of disordered electronic materials are organic semiconductors (OSCs),¹ but we expect that much of what we say here also applies to materials such as hybrid perovskites, conductive metal–organic frameworks and quantum dots.

The difficulty with disordered materials is that they sit in the intermediate regime between the well-understood extremes of band conduction and hopping conduction (Fig. 1a).^1,2 In perfectly ordered crystals, charges move through Bloch waves, wavefunctions that are delocalised over infinitely many sites (individual atoms or molecules). By contrast, in extremely disordered materials (including some OSCs), electronic wavefunctions are localised to one molecule and charges move by thermally assisted hops to neighbouring sites. A theory of the intermediate transport regime must bridge these two qualitatively different extremes.

Fig. 1 (a) While the extreme regimes of transport, coherent band conduction through extended states and incoherent hopping between localised states, are well understood, materials such as organic semiconductors lie in the poorly understood intermediate regime. In this regime, charges hop between partially delocalised electronic states, where the rate of hopping depends on the overlap between the states. (b) We study a system that is a regular lattice of sites, with disordered energies (different colours), that are coupled to neighbours with an electronic coupling J and to an environment (motion lines). (c) Delocalisation of the electronic states is found by diagonalising the system's Hamiltonian. Electronic couplings tend to delocalise the states, while disorder and environmental noise localise them. The resulting delocalisation of the electronic states can be quantified by the inverse participation ratio (IPR) (eqn (24)) or the delocalisation length (l_deloc) (eqn (23)), which are included here for a material with J = 75 meV and σ = 150 meV. When the formation of polarons is accounted for, the states become less delocalised.

In disordered materials, two mechanisms localise electronic states away from infinite Bloch waves (Fig. 1c). The first is Anderson localisation, which is caused by static disorder.³ The second is the formation of polarons due to a carrier's interaction with the environment.^4–6 Either mechanism, if strong enough, can localise states onto individual sites (giving so-called small polarons), but in the intermediate regime the localisation is not complete, and polaron states can be delocalised over multiple sites.

Many OSCs fall into the intermediate regime, which can be most clearly seen from failures of conventional simulations. Transport in most small-molecule OSCs (apart from organic crystals, where the importance of delocalisation has long been recognised⁶) is usually modelled as fully localised hopping, typically through a Gaussian density of states.⁷ The simplest hopping-rate expression is the Miller–Abrahams equation,⁸ which neglects polaron formation, while Marcus theory accommodates both disorder and polaron formation.^9–12 These microscopic theories can be connected to measurable mobilities using kinetic Monte Carlo (KMC) simulations, a probabilistic approach based on averaging stochastic trajectories.¹² KMC often significantly underestimates mobilities, requiring unphysically fast hopping rates (faster than molecular vibrations) to reproduce experimental results. This underestimation occurs because the assumption of completely localised states fails if inter-molecular couplings are comparable to the disorder or to the system-environment coupling, allowing polarons to delocalise across multiple sites.^13–18 An accurate description of intermediate-regime transport in OSCs is particularly important for understanding organic photovoltaics, where delocalisation has been proposed as the key explanation for how charges overcome their Coulomb attraction to achieve rapid and efficient charge separation.^19–29

Nevertheless, describing intermediate-regime transport has proven difficult, particularly in statically disordered systems.² The challenge is that accounting for delocalisation requires a quantum-mechanical treatment, whose computational cost can balloon when disorder prevents periodic boundary conditions being used and forces large simulation boxes instead. In addition, mobility in disordered systems is often governed by deep traps, and long simulation times can be required to reach converged mobilities.³⁰

Existing methods can be broadly divided into atomistic ones and those based on effective Hamiltonians. Atomistic calculations track the dynamics of both the nuclear degrees of freedom (usually using molecular mechanics) and the electronic ones (using quantum equations of motion).^31–37 Atomistic simulations do not have adjustable parameters, but they suffer from the considerable cost of tracking the atomic motion. As a result, taking OSCs as an example, the best atomistic simulations of charge transport are limited to about a thousand molecules, tracked for around 1 ps.^34–37 These capabilities enable remarkable simulations of layered organic crystals, which admit a two-dimensional simulation and are ordered enough that mobilities converge rapidly. However, the same approach cannot be applied to a three-dimensional disordered material that may require a nanosecond-long simulation. By contrast, effective-Hamiltonian models track fewer degrees of freedom,^17,38–43 allowing for larger and longer simulations. These approaches parametrise model Hamiltonians, which then govern the time evolution. Even if the parametrisation is accomplished using atomistic simulations, the important distinction is that an effective-Hamiltonian approach no longer tracks information about the individual atoms after the parametrisation. The weakness of these approaches is that they can neglect important phenomena if they are not included in the model Hamiltonian. For example, effective Hamiltonians designed to treat ordered organic crystals cannot treat disordered materials.

Several effective-Hamiltonian theories incorporate the three ingredients critical to describing the intermediate regime: delocalisation, disorder, and polaron formation.^22,29,43,44 Most of these approaches describe polarons using the polaron transformation, which reduces the otherwise strong system-environment coupling, enabling a perturbative treatment of the remaining interactions. The differences between the approaches are in the details of the effective Hamiltonian and of the perturbative corrections, which give them strengths and weaknesses in different regimes. The Bittner–Silva theory²² is similar to what we do below, except that they study a bath of extended, shared phonons, an assumption that may not be appropriate in molecular systems where phonons are better described as local molecular vibrations.⁴⁵ The Janković–Vukmirović approach²⁹ uses modified Redfield theory as its perturbation theory, making it valid when off-diagonal system-environment couplings are small. However, this condition is not guaranteed in some materials that interest us, including some OSCs; furthermore, disorder is the only localising feature in modified Redfield theory, meaning that the method cannot describe ordered systems where strong system-environment coupling results in small polarons.⁴⁶ The approach of Varvelo et al.⁴³ differs from those above in that it uses a non-perturbative treatment of the system-environment couplings based on the hierarchy equations of motion. In principle, this technique is more accurate in the intermediate regime, but it is computationally more expensive and has so far only been used for one-dimensional chains.

We follow and extend the secular polaron-transformed Redfield equation (sPTRE),⁴⁴ which has several key advantages. Most importantly, sPTRE is entirely in the polaron frame, which changes as a function of the system-environment coupling, allowing sPTRE to characterise intermediate-regime transport as well as exactly reproducing both the band-conduction and hopping-conduction extremes.^44,47–52 The up-front use of the polaron transformation also reduces the delocalisation of the electronic states¹⁸ (Fig. 1c), making mobility calculations easier. Its main limitation has been its computational cost, with sPTRE only ever applied to one-dimensional systems.⁴⁴

Here, we overcome computational roadblocks that have limited sPTRE to one-dimensional systems to present the first three-dimensional description of partially delocalised carriers—whether charges or excitons—in intermediately disordered materials, over times as long as nanoseconds. Our results reproduce sPTRE in one dimension and hopping transport in the low-coupling limit, before showing that even small amounts of delocalisation can dramatically increase mobilities in two and, especially, three dimensions. We also show that these quantum-mechanical enhancements increase at low temperatures due to increasing polaron delocalisation.

I Secular polaron-transformed Redfield equation

Our approach is based on sPTRE,⁴⁴ which we review in this section.

Hamiltonian

We wish to describe an open quantum system, whose total Hamiltonian

H_tot = H_S + H_B + H_SB (1)

consists of components describing the system (H_S), the bath (H_B) and the interaction between them (H_SB). All of the parameters introduced below that enter into H_tot could, in principle, be computed using atomistic simulations that combine molecular mechanics and quantum chemistry.^{1,12,23,29,45}

Our system is a tight-binding model of a d-dimensional cubic lattice of N^d sites, such as molecules or parts of molecules (Fig. 1b). To represent disorder, the energy E_n of each site n is independently drawn from the Gaussian distribution whose standard deviation σ is the disorder of the material. The energetic disorder could arise from static variations in the orientation and spacing of molecules, producing a unique local environment around each molecule. The sites are assumed to be electronically coupled to nearest neighbours with coupling J, which enables delocalisation. We assume a constant nearest-neighbour coupling, although this assumption could easily be relaxed to allow off-diagonal disorder. Overall, the system Hamiltonian is therefore

where |n〉 represents the charge or exciton localised on site n.

We treat the environment as an independent, identical bath on every site, consisting of a series of harmonic oscillators, which can be thought of as vibrations of bonds in the molecules. The bath Hamiltonian is, therefore,

where the kth bath mode attached to the nth site has frequency ω_nk, with creation and annihilation operators b_nk and . Assuming a local bath is common in describing disordered molecular materials;^1,45 in crystalline systems, extended phonons that can couple different sites would be more appropriate.

The interaction between the system and the environment is treated by coupling every site to its bath with couplings g_nk, so that

This linear coupling model is a standard approximation, based on keeping the leading term in the Taylor expansion of a general system-bath interaction.

Polaron transformation

Many materials have electronic (J) or system-bath couplings (g_k) that are too large to be treated as small perturbations. The polaron transformation reduces the system-environment coupling by absorbing it into the polaron itself, permitting the model to be treated using Redfield theory. Polaron formation is described using the state-dependent displacement operator⁶Applying it to the Hamiltonian incorporates lattice distortions into the system by displacing the environmental modes; the polaron-transformed Hamiltonian, indicated by tildes, is then

H̃_tot = e^SH_tote^−S = H̃_S + H̃_B + H̃_SB. (6)

Here, the system Hamiltonian becomes

where and the coupling between sites is renormalised by the factor κ_mn,The bath Hamiltonian remains unchanged, H̃_B = H_B, while the system-bath Hamiltonian becomeswith the new operator

Summing many discrete vibrational modes is computationally costly, so we make two standard simplifications. First, we assume that system-bath couplings are identical at all sites, g_nk = g_k. Second, we assume that the spectral density is a continuous function.⁴⁵ The renormalisation factor then becomes

Because κ < 1, the polaron transformation has two computational advantages: first, it reduces the electronic coupling J, making Redfield theory applicable, and second, it reduces the delocalisation of the electronic states obtained from H̃_S by diagonalisation¹⁸ (Fig. 1c). Here, for concreteness, we adopt the widely used super-ohmic spectral density where λ is the reorganisation energy and ω_c is the cutoff frequency.^53–56 Unless specified otherwise, we use λ = 100 meV, ω_c = 62 meV (ref. 44) and T = 300 K. More generally, the equations above could be used for the more structured spectral densities of organic molecules.

Secular Redfield theory

The polaron transformation reduces the system-bath coupling, allowing H̃_SB to be treated as a perturbation to the system.⁴⁴ Redfield theory is a second-order perturbative approach that, when applied in the polaron frame, results in the polaron-transformed Redfield equation (PTRE). It describes the evolution of the polaron-transformed reduced density matrix [small rho, Greek, tilde] in the basis |μ〉 of polaron states found by diagonalising H̃_S:where ω_μν = E_μ − E_ν. The Redfield tensordescribes the bath-induced relaxation in terms of damping rateswhereis the half-Fourier transform of the bath correlation function⁵⁴where λ_mn,m′n′ = δ_mm′ − δ_mn′ + δ_nn′ − δ_nm′ and

The PTRE of eqn (12) can be further simplified using the secular approximation to give the secular PTRE (sPTRE). The polaron-transformed density matrix consists of diagonal populations and off-diagonal coherences. The evolution of [small rho, Greek, tilde]_μν is controlled by the Redfield tensor, containing terms that transfer populations, dephase coherences, transfer coherences, and mix populations and coherences. In the interaction picture, these terms oscillate with a combined frequency of ω_νμ − ω_ν′μ′. If this frequency is much greater than the inverse of the time frame Δt over which the PTRE is solved, the oscillations are so rapid that the influence of these terms averages out to zero.⁴⁵ The only terms that survive are those for which ω_νμ − ω_ν′μ′ ≪ Δt⁻¹. This condition is met for population transfers and coherence dephasing, and the secular approximation is the assumption that only those terms survive. The result is sPTRE, in which populations and coherences are decoupled.

Furthermore, only populations are relevant for charge transport,⁴⁴ and they are invariant under the polaron transformation, ρ_s(t) = [small rho, Greek, tilde]_s(t), leaving

where the Redfield tensor is now only two dimensional, containing only population transfer termsThe secular approximation reduces computational cost by reducing the number of density-matrix elements from N^2d to N^d and the number of Redfield-tensor elements from N^4d to N^2d.

The secular approximation does not significantly reduce the accuracy for the disordered systems we are studying, for two reasons. First, for parameters typical of the ones we survey, sPTRE agrees well⁴⁴ with a time-convolutionless second-order polaron master equation⁵⁷ that does not use the secular or Markov approximations. Second, any inaccuracies will be minimal for long-time and long-range mobility calculations; at long times, the rate-limiting step in carrier diffusion is the thermal de-trapping from traps that have large energy differences from their neighbours.

Calculating mobilities

The full time evolution of polaron-state populations is given by sPTRE (eqn (18)). Fig. 2a illustrates sPTRE evolution, where the charge density can spread to all other eigenstates in continuous time and in proportion to the corresponding Redfield rate. Eqn (18) has the solution

ρ(t) = exp(Rt)ρ(0), (20)

which can be used to calculate the expectation value of the mean-squared displacement of the charge, 〈r²(t)〉 = tr(r²ρ(t)), at any time.

Fig. 2 The four approximations underlying dKMC. (a) Full sPTRE master equation: the charge density can spread continuously throughout all polaron states. This approach is too expensive in more than one dimension. (b) Kinetic Monte Carlo: individual trajectories are formed from discrete, sequential hops, and are eventually averaged. (c) Hopping radius r_hop: hops are only calculated for states whose centres (black dots) are close enough. (d) Overlap radius r_ove: only sites (grid points) that are close to both the initial and final polaron states are considered in calculating the hopping rate. (e) Diagonalising on the fly: instead of the whole Hamiltonian, only a subsystem of size N^d_box is diagonalised at a time. As the charge moves too close to the boundary, a new Hamiltonian is re-diagonalised centred at the new location of the charge.

To calculate the mobility, 〈r²(t)〉 is averaged over many realisations of disorder (n_iters), i.e., disordered energy landscapes generated using the same microscopic parameters. The resulting average determines the diffusion constant

Finally, for a carrier of charge q, the mobility is given by the Einstein relation

II Delocalised kinetic Monte Carlo

sPTRE has only been applied to one-dimensional systems,⁴⁴ because of three computational hurdles. First, generating the polaron states by diagonalising the N^dH̃_S scales as O(N^3d), where N is the number of sites along each of the d dimensions. Second, tracking the population transfer between all pairs of polaron states involves calculating the full Redfield tensor R_νν′ (eqn (19)), which has N^2d elements. Lastly, each population transfer rate depends on the damping rates Γ (eqn (14)), calculating which involves a sum over N^4d sites to account for the spatial overlap of the polaron states. Therefore, sPTRE scales as O(N^3d) + O(N^6d) overall, which, for reasonably sized lattices, is manageable only for d = 1.

Our approach, dKMC, overcomes these limitations using four approximations (Fig. 2):

(1) Kinetic Monte Carlo

We reduce the number of Redfield rates that need to be calculated by mapping sPTRE onto kinetic Monte Carlo (KMC). Rather than tracking the time-dependent populations of all polaron states, we track stochastic trajectories through the polaron states (Fig. 2b), followed by averaging. This approach mirrors standard KMC, which also probabilistically integrates a large master equation. Individual trajectories are found by hopping to another polaron state with probability proportional to the corresponding Redfield rate for population transfer. Hopping continues until a pre-determined end time t_end. Therefore, instead of calculating Redfield rates for every pair of polaron states, as in sPTRE, KMC requires only calculating outgoing rates at each step. This reduces the number of Redfield rates that need to be calculated from N^2d to N^dn_hopn_traj, where n_hop is the number of hops, which depends on t_end, and n_traj is the number of trajectories, which controls the final averaging error.

To calculate 〈r²(t)〉 for eqn (22), we assume the charge occupying a particular polaron state (at a particular time) is located at its centre, defined as the expectation value of the position, C_ν = 〈ν|r|ν〉; 〈r²(t)〉 is then the average of the square of this displacement over all the trajectories.

(2) Hopping cutoff radius

We reduce the number of Redfield rates to be calculated by introducing a hopping cutoff radius, r_hop. For polaron states that are far away from the current state, the spatial overlaps, and therefore Redfield rates, are very small compared to polaron states that are close by. Therefore, we only calculate rates to polaron states whose centre C_ν lies within r_hop of the centre of the current state (Fig. 2c). The error in this approximation is tunable, because r_hop can be arbitrarily increased depending on the desired accuracy. We choose our r_hop by gradually increasing it by one lattice spacing until the total sum of outgoing rates to states with centres within r_hop converges, not changing by more than a target factor a_hop between increments. The hopping cutoff radius reduces the number of Redfield rates to be calculated at each hop from N^d to O(r^d_hop), thus reducing the total number of rates that need to be calculated in each random energetic landscape to r^d_hopn_hopn_traj.

(3) Overlap cutoff radius

We reduce the cost of calculating individual Redfield rates by introducing an overlap cutoff radius, r_ove. While eigenstates do, in principle, spread across the entire lattice, Anderson localisation predicts that their amplitude decreases exponentially with distance from their centre. Therefore, in calculating damping rates (eqn (14)), we only sum over sites that are simultaneously within a distance of r_ove from the centres of both polaron states (Fig. 2d). Again, the error in this approximation is tunable, and can be decreased arbitrarily by increasing r_ove. We choose our r_ove by gradually increasing it by one lattice spacing until the sum of outgoing Redfield rates calculated by only including sites within r_ove of both the initial and final states converges, not changing by more than a target factor a_ove between increments. The overlap cutoff radius reduces the number of sites summed over in each damping rate from N^4d to O(r^4d_ove).

In practice, we calculate r_hop and r_ove simultaneously, as outlined in Fig. 3, by progressively increasing both until the total sum of the outgoing Redfield rates converges onto a desired accuracy. For simplicity, we choose the two target accuracies to be equal, a_dKMC = a_hop = a_ove.

Fig. 3 The delocalised kinetic Monte Carlo algorithm.

(4) Diagonalising on the fly

We reduce the time required to calculate the polaron states by diagonalising the Hamiltonian on the fly. Rather than diagonalise the entire lattice, we only diagonalise a subset of the Hamiltonian of size N^d_box centred at the location of the charge (Fig. 2e). The charge moves within the box until it gets too close to the edge, which we set to be within r_hop + r_ove of the edge. This buffer ensures we can accurately describe the next hop, which requires a distance of r_hop to contain the centres of all relevant polaron states, as well as a further distance of r_ove to ensure the entirety of polaron states at the edge of r_hop are well defined. In three dimensions, diagonalisation is usually the computational bottleneck, so we make the box as small as possible, N_box = 2(r_hop + r_ove). Once the charge leaves the buffer, the Hamiltonian corresponding to a new box of size N^d_box centred at the new location of the charge is re-diagonalised. The landscape continues to be updated as the charge hops through the material, ultimately reducing the cost of calculating the polaron states from O(N^3d) to O(N^3d_boxn_hop).

Overall, the four approximations above transform sPTRE to dKMC and make it possible to model three-dimensional charge transport in disordered, noisy materials. The detailed steps involved in the algorithm are shown in Fig. 3. Overall, the scaling of the technique has been reduced from sPTRE's O(N^3d) + O(N^6d) to O(N^3d_boxn_hop) + O(r^4d_over^d_hopn_hopn_traj) for dKMC. For example, for a disordered 3D system with J/σ = 0.1 and N = 100, the scaling is reduced by at least 25 orders of magnitude.

III Results and discussion

Accuracy

All of the approximations in dKMC are controllable, meaning that the error can be arbitrarily reduced given additional computational resources. This accuracy can be demonstrated by comparing dKMC mobilities to those predicted by sPTRE in one dimension, and Fig. 4a shows the agreement increasing with the accuracy parameter a_dKMC. The agreement between sPTRE and dKMC at a_dKMC = 0.99 leads us to adopt that value throughout this paper. Cutoff radii always lead to an underestimation of delocalisation effects, meaning that the substantial delocalisation enhancements reported below are strictly lower limits.

Fig. 4 (a) Mobilities (μ) in one dimension at 100 ps predicted by both sPTRE and dKMC, with disorder σ = 150 meV and shown as a function of the electronic coupling J. dKMC can be made arbitrarily accurate by increasing the accuracy parameter a_dKMC (a_dKMC = 0.99 is used elsewhere in this paper). (b) Throughout this paper, mobilities—which are time-dependent³⁰—are calculated at 100 ps. The mobility enhancement due to delocalisation (μ_dKMC/μ_KMC) increases with the time (here shown for one dimension), meaning that longer-time delocalisation enhancements would always be at least as large as reported here.

It is also necessary to choose the time cutoff t_end, because mobilities in disordered materials are time dependent (or dispersive), and it can take a long time to converge on a steady-state mobility.³⁰ The ultimate choice will depend on particular applications; we report mobilities at t_end = 100 ps, which corresponds to charge transit times on typical length scales in OSCs (tens of nanometres). For our purposes, the important fact is that mobility enhancements due to delocalisation at longer times are always at least as large as at t_end (Fig. 4b). At longer times, charges are increasingly likely to get stuck in deeper traps, causing the mobility to decrease with time. Delocalisation allows the wavefunction to leak onto neighbouring sites, helping the detrapping and giving larger enhancements at longer times.

Importance of 3D effects

dKMC captures effects missing in lower-dimensional approximations, in particular the extent of delocalisation of the polaron states. We define the delocalisation lengthwhere we average the inverse participation ratios of the polaron states,which indicates the number of sites n that a state ν extends over.⁵⁸ Therefore, l_deloc describes the size of a state along each dimension, and is a way of comparing the extent of delocalisation across different dimensions.

Fig. 5a shows that l_deloc increases as a function of J in all three dimensions, as expected. More interestingly, l_deloc is significantly larger in three dimensions than in one or two, with the difference becoming larger with increasing J. Therefore, including all three dimensions is essential for modelling intermediately delocalised charge transport, and lower-dimensional models may significantly underestimate delocalisation effects.

Fig. 5 (a) Delocalisation length (l_deloc) in each dimension as a function of the electronic coupling (J), for disorder σ = 150 meV. All else being equal, the delocalisation length is considerably greater in three dimensions than in lower dimensions. (b) Mobilities (μ) at 100 ps calculated by dKMC and standard KMC as a function of J, with σ = 150 meV. We only include KMC in three dimensions for legibility, because KMC mobilities in the lower dimensions are similar in magnitude. KMC and dKMC agree in the low-coupling limit when states are localised, but as the electronic coupling and delocalisation increase, so too do the mobilities predicted by dKMC compared to those predicted by KMC.

Delocalisation enhances mobility

Our most important finding is that even modest delocalisation dramatically enhances mobilities, meaning that delocalisation is critical for explaining transport in the intermediate regime. As J increases, the increase in delocalisation increases overlaps between states and, therefore, the Redfield transfer rates and the ultimate mobilities (Fig. 5b). For values of J and σ that are reasonable for OSCs, including delocalisation using dKMC can increase mobilities by close to an order of magnitude above the localised-hopping of standard KMC, which helps explain why mobilities predicted by KMC are usually too low compared to experiment. Furthermore, these enhancements require only a small amount of delocalisation, less than 2 sites in each direction in 2D.

Fig. 6 The mobility enhancement (μ_dKMC/μ_KMC), shown here in two dimensions for J = 45 meV and σ = 150 meV, increases as the temperature is reduced, due to the increase in delocalisation of the polaron states (l_deloc).

Larger enhancements at lower temperatures

Finally, the delocalisation mobility enhancement is temperature dependent, with larger enhancements at low temperatures (Fig. 6). The higher mobilities at lower temperatures are caused by the increased delocalisation of the polaron states. Applying the polaron transformation reduces J by multiplying by a factor of κ (eqn (8)), where κ < 1. κ increases as T is lowered and, therefore, J is reduced less at low T than at higher T, allowing the polaron states to delocalise further and assisting their mobility.

Outlook

The most immediate area for future work is the application of the newly developed theory and the important trends identified in this paper to the prediction of experimental mobilities in concrete disordered materials. As we have shown, standard KMC approaches underestimate mobilities; therefore, we expect that including delocalisation will improve experimental agreement. However, connecting dKMC to experimental mobilities still faces some challenges, especially the difficulty of obtaining good estimates of disorder, whether experimentally⁵⁹ or using ab initio calculations.

The main limitation of dKMC remains computational cost; in particular, in Fig. 5b, 3D dKMC is limited to modest values of J. As the states become larger, the Hamiltonian box needs to be increased to accurately capture the states and allow them to move around. The large box is expensive to diagonalise, and it requires the calculation of more rates, with more sites contributing to the overlap calculations in each rate. Nevertheless, there is a clear trend, and we expect the importance of delocalisation to be even more pronounced at high J in 3D. In the future, it may be possible to reduce the computational cost further, either through additional approximations, or by identifying robust trends in the numerical results that can then be extrapolated.

dKMC is also limited by the approximations made in sPTRE, some of which could be relaxed using other master equations. sPTRE makes the Markov and secular approximations, and it can be inaccurate for systems that are weakly coupled to slow environments.^44,47,53 The parameter regimes we studied fall within the range of validity of sPTRE established in previous work.⁴⁴ However, extending dKMC to weakly coupled slow baths would require modifications to the underlying sPTRE, for example by incorporating a variational polaron transformation^53,57,60 instead of the fully displaced version used in sPTRE. The variational treatment would also enable the treatment of ohmic and sub-ohmic baths, unlike the super-ohmic ones assumed in sPTRE.

Furthermore, we expect that it will be possible to extend dKMC to other commonly encountered situations. The most straightforward extensions would include the prediction of mobilities at high charge densities, in the presence of external electric fields, on irregular or anisotropic lattices, or in spatially constrained domains. It may also be possible to extend dKMC to describe the more difficult problem of charge separation of excitons in organic photovoltaics. Because charge separation is a two-body problem involving the correlated motion of an electron and a hole, the computational difficulty is roughly the square of the single-body mobility calculation, meaning that a fully quantum-mechanical treatment has so far proved intractable in three dimensions.²³ We expect that dKMC will make this problem computationally accessible, allowing the first simulation of the full dynamics (and, therefore, efficiency) of charge separation in the presence of disorder, delocalisation, and noise. A complete kinetic model would help settle the debate about the main drivers of charge separation, and unite the proposed mechanisms including delocalisation,²⁰ entropy^61,62 and energy gradients.⁶³

IV Conclusions

dKMC is the first approach able to describe charge transport in intermediately disordered materials in three dimensions. It keeps the benefits of sPTRE—fully quantum dynamics, accurate treatment of polarons, and the ability to reproduce both extremes of transport—while overcoming computational obstacles that have prevented sPTRE from being used in more than one dimension. We have used dKMC to capture the effects of delocalisation and show that carrier mobilities are significantly higher than those predicted by standard KMC. Indeed, even small amounts of delocalisation—less than two sites—can increase mobilities by an order of magnitude. All of these mobility enhancements are greater at lower temperatures, due to the increased delocalisation of polaron states. In the future, we expect that dKMC can be extended to a wider range of systems, shedding even more insight into fundamental charge- and exciton–transport processes.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We were supported by a Westpac Scholars Trust Research Fellowship, a Westpac Scholars Trust Future Leaders Scholarship, an Australian Government Research Training Program Scholarship, and a University of Sydney Nano Institute Grand Challenge. We were supported by computational resources and assistance from the National Computational Infrastructure (NCI), which is supported by the Australian Government, and by the Sydney Informatics Hub and the University of Sydney's high-performance computing cluster Artemis.

References

1. A. Köhler and H. Bässler, Electronic Processes in Organic Semiconductors: An Introduction, Wiley, 2015.
2. H. Oberhofer, K. Reuter and J. Blumberger, Chem. Rev., 2017, 117, 10319.
3. P. W. Anderson, Phys. Rev., 1958, 109, 1492.
4. H. Fröhlich, Adv. Phys., 1954, 3, 325.
5. T. Holstein, Ann. Phys., 1959, 8, 343.
6. M. Grover and R. Silbey, J. Chem. Phys., 1971, 54, 4843.
7. H. Bässler, Phys. Status Solidi B, 1993, 175, 15.
8. A. Miller and E. Abrahams, Phys. Rev., 1960, 120, 745.
9. R. A. Marcus, J. Chem. Phys., 1956, 24, 966.
10. S. Athanasopoulos, J. Kirkpatrick, D. Martínez, J. M. Frost, C. M. Foden, A. B. Walker and J. Nelson, Nano Lett., 2007, 7, 1785.
11. I. I. Fishchuk, A. Kadashchuk, S. T. Hoffmann, S. Athanasopoulos, J. Genoe, H. Bässler and A. Köhler, Phys. Rev. B, 2013, 88, 125202.
12. V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey and J. L. Brédas, Chem. Rev., 2007, 107, 926.
13. M. Bednarz, V. A. Malyshev and J. Knoester, J. Chem. Phys., 2004, 120, 3827.
14. H. Oberhofer and J. Blumberger, Phys. Chem. Chem. Phys., 2012, 14, 13846.
15. F. C. Grozema and L. D. Siebbeles, Int. Rev. Phys. Chem., 2008, 27, 87.
16. H. Yang, F. Gajdos and J. Blumberger, J. Phys. Chem. C, 2017, 121, 7689.
17. C. Liu, K. Huang, W. T. Park, M. Li, T. Yang, X. Liu, L. Liang, T. Minari and Y. Y. Noh, Mater. Horiz., 2017, 4, 608.
18. B. Rice, A. A. Y. Guilbert, J. M. Frost and J. Nelson, J. Phys. Chem. Lett., 2018, 9, 6616.
19. H. Tamura and I. Burghardt, J. Am. Chem. Soc., 2013, 135, 16364.
20. S. Gélinas, A. Rao, A. Kumar, S. L. Smith, A. W. Chin, J. Clark, T. S. van der Poll, G. C. Bazan and R. H. Friend, Science, 2014, 343, 512.
21. S. Valleau, S. K. Saikin, M. H. Yung and A. A. Guzik, J. Chem. Phys., 2012, 137, 034109.
22. E. R. Bittner and C. Silva, Nat. Commun., 2014, 5, 3119.
23. S. Few, J. M. Frost and J. Nelson, Phys. Chem. Chem. Phys., 2015, 17, 2311.
24. M. Huix-Rotllant, H. Tamura and I. Burghardt, J. Phys. Chem. Lett., 2015, 6, 1702.
25. H. Bässler and A. Kohler, Phys. Chem. Chem. Phys., 2015, 17, 28451.
26. A. Gluchowski, K. L. Gray, S. N. Hood and I. Kassal, J. Phys. Chem. Lett., 2018, 9, 1359.
27. V. Janković and N. Vukmirović, J. Phys. Chem. C, 2018, 122, 10343–10359.
28. S. Athanasopoulos, H. Bässler and A. Köhler, J. Phys. Chem. Lett., 2019, 10, 7107.
29. V. Janković and N. Vukmirović, J. Phys. Chem. C, 2020, 124, 4378.
30. S. T. Hoffmann, S. Athanasopoulos, D. Beljonne, H. Bässler and A. Köhler, J. Phys. Chem. C, 2012, 116, 16371.
31. A. Troisi and G. Orlandi, Phys. Rev. Lett., 2006, 96, 086601.
32. S. Fratini, D. Mayou and S. Ciuchi, Adv. Funct. Mater., 2016, 26, 2292.
33. A. Heck, J. J. Kranz and M. Elstner, J. Chem. Theory Comput., 2016, 12, 3087.
34. S. Giannini, A. Carof and J. Blumberger, J. Phys. Chem. Lett., 2018, 9, 3116.
35. S. Giannini, A. Carof, M. Ellis, H. Yang, O. G. Ziogos, S. Ghosh and J. Blumberger, Nat. Commun., 2019, 10, 3843.
36. S. Giannini, O. G. Ziogos, A. Carof, M. Ellis and J. Blumberger, Adv. Theory Simul., 2020, 3, 2000093.
37. O. G. Ziogos, S. Giannini, M. Ellis and J. Blumberger, J. Mater. Chem. C, 2020, 8, 1054.
38. B. M. Savoie, K. L. Kohlstedt, N. E. Jackson, L. X. Chen, M. O. De La Cruz, G. C. Schatz, T. J. Marks and M. A. Ratner, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 10055.
39. N. E. Jackson, L. X. Chen and M. A. Ratner, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 8595.
40. Y. Jiang, X. Zhong, W. Shi, Q. Peng, H. Geng, Y. Zhao and Z. Shuai, Nanoscale Horiz., 2016, 1, 53.
41. T. Nematiaram, S. Ciuchi, X. Xie, S. Fratini and A. Troisi, J. Phys. Chem. C, 2019, 123, 6989.
42. A. Y. Sosorev, Mater. Des., 2020, 192, 108730.
43. L. Varvelo, J. K. Lynd and D. I. G. Bennett, 2020, arXiv:2008.06496.
44. C. K. Lee, J. Moix and J. Cao, J. Chem. Phys., 2015, 142, 164103.
45. V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, Wiley-VCH, Weinheim, 3rd edn, 2011.
46. S. Jesenko and M. Žnidarič, 2014, arXiv:1405.4156.
47. C. K. Lee, J. Moix and J. Cao, J. Chem. Phys., 2012, 136, 204120.
48. D. Xu and J. Cao, Frontiers of Physics, 2016, 11, 110308.
49. H. T. Chang, P. P. Zhang and Y. C. Cheng, J. Chem. Phys., 2013, 139, 224112.
50. S. Jang, Y. C. Cheng, D. R. Reichman and J. D. Eaves, J. Chem. Phys., 2008, 129, 101104.
51. A. Nazir, Phys. Rev. Lett., 2009, 103, 146404.
52. D. P. S. McCutcheon and A. Nazir, New J. Phys., 2010, 12, 113042.
53. F. A. Pollock, D. P. McCutcheon, B. W. Lovett, E. M. Gauger and A. Nazir, New J. Phys., 2013, 15, 075018.
54. S. Jang, J. Chem. Phys., 2011, 135, 034105.
55. S. Jang, J. Cao and R. J. Silbey, J. Phys. Chem. B, 2002, 106, 8313.
56. E. Y. Wilner, H. Wang, M. Thoss and E. Rabani, Phys. Rev. B, 2015, 92, 44.
57. E. N. Zimanyi and R. J. Silbey, Philos. Trans. R. Soc., A, 2012, 370, 3620.
58. J. M. Moix, M. Khasin and J. Cao, New J. Phys., 2013, 15, 085010.
59. S. Hood, N. Zarrabi, P. Meredith, I. Kassal and A. Armin, J. Phys. Chem. Lett., 2019, 10, 3863.
60. R. Silbey and R. A. Harris, J. Chem. Phys., 1984, 80, 2615.
61. T. M. Clarke and J. R. Durrant, Chem. Rev., 2010, 110, 6736.
62. S. N. Hood and I. Kassal, J. Phys. Chem. Lett., 2016, 7, 4495.
63. F. C. Jamieson, E. B. Domingo, T. McCarthy-Ward, M. Heeney, N. Stingelin and J. R. Durrant, Chem. Sci., 2012, 3, 485.

---