Identifying pitfalls when using the Miller–Abrahams rate in kinetic Monte Carlo simulations

Abstract

The conventional Miller–Abrahams (MA) rate is a frequently applied rate for modeling the hopping transport of charges or excitons in organic semiconductors. However, the expression known as Miller–Abrahams rate is an approximation that has a more limited range of validity than the original, full expression. We study the effect of both rates in Kinetic Monte Carlo simulations on the resulting charge carrier mobility in OFET and OLED structures. We find significant differences for small disorders as well as for high electric fields. While the original, full expression predicts an increase and finally saturation of the mobility with temperature and field, the conventional (approximated) MA rate erroneously yields a decreasing mobility with temperature and field. We demonstrate that this results from the constant rate for energetically downward carrier jumps in the conventional, approximated MA rate in contrast to the full rate where downwards jumps accelerate with increasing energy difference and discuss the physical origin of this. This aspect becomes relevant when downward jumps are rate-limiting, e.g. for small disorders or high fields.

---

Introduction

A central parameter controlling the performance of organic semiconductor devices such as transistors, light-emitting diodes and solar cells is the charge carrier mobility.¹ Meanwhile, mobilities in the range of several cm² V⁻¹ s⁻¹ can be obtained for OFET or OLED applications when the molecular or polymeric films concerned are, first, free of extrinsic traps and, second, have a high degree of structural order.^2–5 This can be achieved in multiple ways. In solution-processed films, suitable processing methods, such as meniscus-guided coating, and the choice of solvent have been shown to increase the films’ order. Concerning the materials themselves, appropriate choices of material parameters, such as regioregularity and molecular weight are known to enhance molecular ordering. Finally, for vacuum deposited films modifying the evaporation rate and substrate temperature can similarly lead to ordered films.^6–9 Today's high charge carrier mobilities are a direct consequence of higher morphological order, e.g. in films with crystalline domains or semi-crystalline microstructures, and the accompanying low energetic disorder.

Such films show, however, not only higher values in the charge carrier mobilities, but also a temperature dependence that is no longer consistent with the commonly used description in terms of a thermally activated hopping transport that can be modelled using a conventional Miller–Abrahams rate.⁵ Some well-ordered organic semiconductors with high charge carrier mobility show a falling mobility with increasing temperature, consistent with a band-like transport of delocalized charge carriers. At the same time other experiments, e.g. optical-pump terahertz probe spectroscopy which yields the optical conductivity or charge modulation spectroscopy which gives the charge induced absorption, suggest a localized character of the charges.⁵ In these compounds, coupling to intermolecular phonon modes is strong and modifies the charge transfer process. The resulting mobility can then be described in terms of transient localization or delocalization mechanism.^5,10–13

Some crystalline materials, however, show a constant mobility over a large temperature range, at odds with both the conventional MA-type hopping transport and the description in terms of transient localization. Examples include transport in polyacenes along the c-axis, where the electronic coupling is weak, e.g. electron transport in naphthalene and anthracene.^14,15 For electron transport in naphthalene along the c-axis, the mobility decreases from about 50 to 100 K and then remains constant up to 325 K.¹⁴ Similarly, a completely temperature independent electron transport prevails in anthracene along the c-axis from 77 K to 375 K.¹⁵ Clearly, this is not compatible with a description in terms of band-like delocalized states. To understand why the description of a localized hopping type transport also seems to fail, it is worthwhile to take a closer look at the conventional Miller–Abrahams rate. While many different rate expressions have been developed and discussed over the last decades,^16–20 the Miller–Abrahams (MA) rate has been used particularly frequently.²¹ It is conventionally expressed in the following form:^22,23

Here, ν₀ is the attempt-to-hop-frequency, γ is the inverse localization radius, r is the distance and ΔE = ε_final − ε_initial is the energy difference between the final and initial site energy ε. Eqn (1) is frequently used in simulations and theoretical work to model charge transfer in organic semiconductors.^21,24–26 It forms the foundation for the Gaussian disorder model and its various modifications when it is combined with a master equation approach, a Gaussian density of states and notion that charge transport occurs through thermally activated hops from an equilibrium energy to a transport energy.^24,27–30 In the field of organic semiconductors, its popularity derives from the fact that it is simple, and that it was very successful in describing the field and temperature dependence of charge carrier mobility in the disordered semiconductors that had been used in the past.^24,27–30Eqn (1) has also been used as a starting point to develop expressions that, in addition to energetic disorder, include the reorganization energy λ through a weighting factor (C-parameter).^31–33 Furthermore, it is applied to describe charge transport in inorganic semiconductors,³⁴ hybrid perovskites,³⁵ or biomolecules such as DNA³⁶ or to model random walks of Brownian particles.³⁷ When combined with the Förster transfer rate, the resulting expression is used to describe singlet exciton transfer, e.g. in organic solar cells or light harvesting complexes.^38,39 An independently derived but similar expression is used in the general Metropolis Monte Carlo algorithm to calculate the acceptance probability of a random change to the modeled system.⁴⁰

The Miller Abrahams (MA) rate has already been discussed and compared to other possible hopping rates.^{17,20,21,41–46} For example, Vukmirović and Wang emphasized the need to consider the phonon density of states.²⁰ Also, at strong electric fields⁴⁷ or at very small disorder strengths of 3 meV,⁴⁸eqn (1) was observed to result in an inverse temperature dependence of the mobility that contradicts the original assumption of thermally assisted hopping. Obviously, there are cases in which eqn (1) has shortcomings.

Its long string of successes leads to a neglect of the fact that the MA rate as given in eqn (1) is, in fact, an approximation of a more general expression derived by Miller and Abrahams for weakly doped semiconductors. To derive the non-approximated, full expression, Miller and Abrahams considered that charge transfer is associated with the absorption (ΔE > 0) or emission (ΔE < 0) of a single acoustic phonon.²² Hence, ΔE should not exceed the maximum Debye energy of single acoustic phonons. This satisfies energy conservation, yet, as detailed beautifully by Fornari et al., the key aspect is that the single (non-local) phonon mode actually enables the transition between the two sites to occur, similar to a Herzberg–Teller coupling mechanism.⁴⁹ For uncorrelated disorder, the final rate expression of Miller and Abrahams can be written as⁵⁰

where ħ is the Planck constant and A is a dimensionless constant. In systems, in which the mean energy differences do not differ greatly, A is connected to the constant attempt-to-hop-frequency ν₀ from eqn (1) by ν₀ = A〈|ΔE|〉ħ. In the case of |E| ≫ 2k_BT, the conventionally applied version of the MA rate, i.e., eqn (1) is recovered^23,50 and shall henceforth be designated as “approximate MA rate”. This rate is still used in many simulations, even commercial ones, because it is simple and depends on only a few parameters.^21,25,26

Miller and Abraham's final rate expression in the form of eqn (2) was derived by Fishchuk et al. and compared to eqn (1).⁵⁰ Fishchuk's work is an analytic treatment based on the effective-medium approximation approach and the concept of transport energy. It was motivated by the intention to explain the unusually weak temperature dependence of the hole mobility in a film of the only weakly disordered polymer (MeLPPP).⁵⁰ Here we use kinetic Monte Carlo simulations to explore whether some of the above-mentioned shortcomings of the conventional, i.e. approximate MA rate (eqn (1)) may result from using it outside the parameter range for which its derivation is valid, so that the full rate should be used. A key insight is that for small disorder, and hence small ΔE, the energy-dependent expressions in eqn (1) and (2) become constant, so that the rate k is determined by the prefactor. This prefactor is constant for eqn (1) yet directly proportional to T for eqn (2). The hopping rate k directly yields the diffusivity D, from which the mobility μ is obtained as by means of the Einstein–Smoluchowski equation. A temperature-independent mobility, as experimentally observed in systems with a small energetic disorder, can therefore only be obtained when using the full rate, eqn (2), so that the temperature dependencies cancel. We focus on discussing the implications and the physical origin of differences between the full and the approximate rate.

Methods

In the following, we describe our simulation model, the geometry of the modelled devices and the procedure to obtain the charge carrier mobility from the simulation results. We simulate two different device geometries: a sandwich-type geometry, used typically for modelling organic light emitting diodes (OLEDs), with a simulation box size of n_x × n_y × n_z = 100 nm × 50 nm × 20 nm and a device, which differs from the other geometry by its height n_z = 3 nm and an additional gate voltage field along the z-axis, thus representing an organic field transistor (OFET), cf.Fig. 1. The small height in the OFET case is still reasonable since most charge carriers occupy only the lower few layers in an OFET at high gate voltages.^26,51 In both cases, we apply periodic boundary conditions along the x- and y-axes and model the site energies E as uncorrelated disorder drawn from a Gaussian distribution.²⁴ Therefore, the density of electronic states (DOS) is of the formwhere the standard deviation σ characterizes the energetic disorder, taking values in our simulations of 5–100 meV. The uncorrelated disorder was chosen to allow for direct comparison with the analytical calculation of Fishchuk et al. that is valid only for the uncorrelated disorder.⁵⁰ Results obtained for the correlated disorder are added for reference in the ESI‡ (Fig. S.1–S.3). They are qualitatively similar, yet the difference between using the full and the approximate rate is more pronounced.

Fig. 1 Model geometries used in our simulations for the (a) OLED and (b) OFET geometry. In both cases, only transport in the semiconductor (shown in blue) is modelled with periodic boundary conditions along x- and y-axis.

For both the OLED and the OFET cases, an external electric field F is introduced in the x-direction. The resulting potential drop over a single lattice spacing (a = 1 nm) is given by eFa. For each hop with distance a in the field direction, this potential difference is added to the energy difference ΔE between the final and initial site. For hops against field direction, ΔE is subtracted. Because the potential difference eFa is defined relative to the direction of motion, it remains valid even when a carrier crosses the periodic boundary along the x-axis. The coulomb interaction between the charges (carrier–carrier interactions) is considered explicitly by using the minimum image convention and by updating the changes in the coulomb potential, as it was suggested by Li and Brédas.⁵² The rates for all charges to hop to neighboring sites within a maximum hopping radius of are calculated by eqn (1) with ν₀ = 10¹³ s⁻¹, or eqn (2) with A = 1 ≈ 24.18 eV⁻¹ħν₀, respectively. The inverse localisation radius was γ = 5 nm⁻¹ for all cases. Then, one transfer rate is randomly chosen, and the corresponding event is executed. The corresponding time interval, in which this event occurs, is calculated by

where X is a random number.²⁶ The simulation terminates when the current of carriers along the x-axis has converged. We run simulations with every parameter set until at least 100 charge carriers were used in total. The mobility then was calculated from data after convergence bywhere Δx is the net distance all carriers travelled in field direction in the time interval Δt, F is the field created by the applied source–drain voltage and N is the number of carriers, which is set constant to N = 10 in the OLED-type device and proportional to the gate voltage in the OFET device.¹

Results

We simulated the charge carrier mobility in systems with OFET as well as OLED geometry using the full and approximate MA rates, parametric in disorder, temperature and charge concentration. We shall use μ_full and μ_approx to differentiate the mobilities obtained with both rates.

Temperature dependence

OFET-structure model. First, we focus on the transfer behaviour simulated with the OFET-structure model. Fig. 2 shows the dependence of the mobility on the charge carrier concentration and temperature for different disorder strengths σ = 5, 10, 25, and 50 meV, comparing μ_full on the left-hand-side to μ_approx on the right-hand-side. At high disorders, such as 50 meV, the mobilities obtained from both rates exhibit a similar behavior, i.e., the mobility increases with temperature and with charge carrier concentration. The mobility increase with temperature reflects the increased available thermal activation energy. The increase in carrier density results from the increase in the Fermi level with gate voltage, which reduces the activation energy needed for charge transport. As the disorder reduces, the temperature dependence also reduces since less activation energy is required for transport. We call to mind that a value of σ = 25 meV corresponds to the thermal energy at room temperature. For very small disorders, e.g. σ < 10 meV, however, μ_approx starts to deviate significantly from μ_full. While for μ_full, the temperature dependence of the mobility reduces further to yield a nearly temperature independent mobility, for μ_approx the temperature dependence of the mobility reverses, which would imply an increasing mobility upon cooling. We also find that the curves at 400 K happen to agree for both rates; however, this is a coincidence due to the choice of the constants A and ν₀.

Fig. 2 Dependence of the carrier mobility on the gate voltage in an OFET geometry, parametrized with different temperatures, shown for several disorders (from top to bottom: 5 meV, 10 meV, 25 meV, and 50 meV), derived with the full rate (left-hand-side, μ_full) or the approximate rate (right-hand-side, μ_approx).

OLED-structure model. We found the same temperature evolution of the mobility in simulations with the OLED-geometry. The mobility obtained from both rates at different disorder strengths and temperatures at a field strength of E = 0.1 MV cm⁻¹ is displayed in an Arrhenius form in Fig. 3(a). For low temperatures at 75 and 100 meV disorder, the simulations did not produce reliable data in a reasonable computational time. Analogous to the OFET structure, at high disorder values, mobilities obtained from both rates exhibit a similar temperature dependence, i.e. increasing mobility with temperature, albeit the values differ slightly (cf. ESI,‡ Fig. S.4). Like for the OFET case, μ_full is almost independent of temperature for disorder values at 10 and 5 meV, while μ_approx decreases with temperature, thus showing an inverse temperature dependence.

Fig. 3 Temperature dependence of the mobility for different σ = 5, 10, 25, 50, 75, and 100 meV in an OLED-type device using either the full Miller–Abrahams hopping rate or the approximate expression. At high mobilities, the confidence intervals are smaller than the symbol size. (a) Results in an Arrhenius representation. The auxiliary line at 10⁻⁴ cm² V⁻¹ s⁻¹ serves to differentiate the different temperature dependencies for 5 and 10 meV disorder for both rates. (b) Mobility plotted against the relative disorder (σ/k_BT)². Beware that some data points pertaining to different values of σ coincide. The grey line shows an exponential decay with a decay constant of (2/3)².

To investigate the differences between both rates in their dependence on disorder and temperature, the data are additionally plotted against the quadratic relative disorder (σ/k_BT)² in Fig. 3(b). For μ_full, the data points obtained for different disorder strengths coincide when plotted against (σ/k_BT)² and the resulting curve closely follows an exponential decay , indicated as a straight grey line (see also Fig. S.5, ESI‡). The deviation at large values of (σ/k_BT)², i.e. for σ/k_BT larger than 3, is due to the finite size of the simulation box. States at the very bottom of the density of states are, by definition, rare. For a finite box size, these very low probability sites are no longer present in the distribution. This prevents relaxation of the charge carriers into tail states, so that equilibrium is no longer established.⁵³

In contrast, the results obtained with the approximate rate deviate strongly from this behaviour. Importantly, the σ/k_BT scaling is no longer valid, and it makes a difference whether the temperature or the disorder is changed. The effect is stronger at smaller disorders.

Field dependence

We also compared the dependence of the mobility on the applied field and disorder strength. The mobility μ_full, shown in the upper plot of Fig. 4, increases with field until it saturates and converges towards a maximum mobility for high fields. This implies that the drift velocity is increasing linearly with the field. In contrast, μ_approx, shown in the lower plot of Fig. 4, first increases at high disorders yet then decreases with 1/F when the drain field is increased further.

Fig. 4 Dependence of the mobility on the electric field in an OLED-type device at T = 300 K for the full (top) and approximate (bottom) Miller–Abrahams rate.

Discussion

The key results of the current simulation study are as follows: (i) in contrast to the full rate, the approximate MA rate predicts a decrease of the mobility with increasing temperature for small disorder values, (ii) the mobility scales with σ/k_BT for the full MA rate yet not the approximate rate, and (iii) the full rate predicts an asymptotic approach to a finite maximum mobility with the electric field, while the approximate rate suggests the mobility to decrease at high electric fields. Such a decrease has indeed been observed yet there are several ways to account for it.⁵⁴

In the following, we shall address the conceptional difference between the full and the approximate MA rates and the implication for analyzing charge transport in a disordered organic semiconductor. In Fig. 5, we illustrate the difference between the jump rates for a charge carrier in a hopping system as a function of (a) energetic difference between charge donating and charge accepting states and (b) temperature.

Fig. 5 The transfer rate k, divided by the prefactor ν₀exp(−2γr), for the full MA rate (eqn (2)) in green and the approximate MA rate (eqn (1)) in blue, shown for jumps downwards in energy (left) and upwards in energy (right), (a) dependent on the energy difference between the initial and final site at k_BT = 25 meV, and (b) as a function of temperature at ΔE = 5 meV.

Regarding the energy dependence, the upward jumps are thermally activated, tending to a Boltzmann-term for both rates at large energy differences, since in eqn (2) sinh(x) ≈ 0.5e^x for x ≫ 1. The key difference lies in the downward jumps. In the approximate rate, any energy dependence is ignored so that jumps downward in energy are all taken to be equally likely, while downward jumps accelerate with the energy gap for the full rate. This is quite in contrast to the energy gap law,⁵⁵ as already noted by Fornari et al.⁴⁹

The temperature dependence for both rates is also strikingly different. For the full rate, the upward and downward jump rates both increase roughly linear with temperature. They differ mainly at very low temperatures where the downward jump rate remains finite in contrast to the vanishing upward jump rate. For the approximate rate, however, the downward jump rate remains temperature-independent while the upward jump rate merely reflects the temperature evolution of the Boltzmann factor. In the high temperature limit, the rates for both jump directions of the approximate rate become constant.

To extract the difference in the underlying physics, we turn to the derivation of a very general rate expression that was developed by Fornari et al.⁴⁹ They emphasize that the coupling between the initial and final site is brought about by the coupling to inducing vibrational modes and derive a general rate expression in which the transfer rate k depends on the product of the spectral density J(ω) and the phonon occupation number N(ω), see the ESI‡ Section D. The basic idea is that phonon modes m with energy ω^I_m modulate the coupling and induce the transition. |M_12,m|² is connected to the transition matrix element V₁₂ by the displacement q_mvia V₁₂ = |M_12,m|²q_m. Fornari et al. define . The spectral density J(ω) is therefore a measure for the coupling strength per mode. The phonon occupation number N(ω) is given by the Bose–Einstein statistics.

Within the framework of Fornari et al., eqn (1) and (2) can be obtained if one assumes a certain form for the spectral density J(ω). The full MA rate, eqn (2), implies a spectral density that increases linearly with energy. In contrast, a spectral density that saturates underlies the approximate MA rate (cf. ESI‡ Section D). The linearly increasing rate with the downward energy gap for the full rate may hence be associated with the increasing spectral density implicit in the full rate, while the constant downward jump rate in the approximate rate may be seen as connected to the saturating spectral density.

Trends in the dependence on temperature and disorder

Having considered how the rates differ between the full and the approximate MA model, we can discuss the observed opposite trends of the mobility obtained from both rates. Let us first consider the temperature dependence shown in Fig. 2 and 3. We recall that for low or moderate applied fields (in contrast to the high field limit discussed in the context of Fig. 4), we found a similar temperature-dependence of μ_approx compared to μ_full at high energetic disorders, yet an inversed one for small disorders. Moreover, curves derived for different disorder strengths coincide when plotted against (σ/k_BT)² for μ_full over the entire parameter range yet not for μ_approx. The reduction of mobility upon increasing the temperature for μ_approx at small disorders had previously already been noted by Chatten et al.⁴⁸ in KMC simulations carried out for σ = 3 meV. The failure of μ_approx to scale with (σ/k_BT)² at low disorder values was already highlighted by Fishchuk et al.⁵⁰ in their comparison of μ_full with μ_approx that was performed using effective-medium theory.

To understand these related phenomena, it is useful to start by reflecting on the case where energetic disorder is large. For large disorder, full and approximate rates yield very similar results. Large uncorrelated disorders trivially imply large energy gaps between adjacent sites, i.e. |ΔE| ≫ 2k_BT, so that the full and the approximate rate are both limited by the upward jumps, i.e., by the Boltzmann term. More precisely, charge transport in a disordered density of states (DOS) is limited by upward jumps from a temperature-dependent equilibrium energy to a transport level slightly below the center of the DOS. This is the origin of the exponential dependence of mobility on (σ/k_BT)² that is observed experimentally in systems with moderate to large disorder.^24,56 This behavior is also reflected in our data for the regime (cf. grey line in Fig. 3b). For larger values of , the finite simulation box precludes full relaxation of the simulation data to the expected grey line as detailed above.

Mobilities obtained with both rates begin to differ in temperature evolution and value as the disorder reduces. The reason is that thermal activation as well as the concept of the transport level both lose their relevance in the case of low disorder since there is nearly no difference between the energies of neighboring sites. For low disorder, upward and downward jumps attain similar weight. The transport is then well described as normal diffusion on a three-dimensional grid, in which case the diffusion coefficient D is given as:⁵⁷

Here, a is the lattice constant and k is the transfer probability. Although this expression was derived for random-walk steps only to directly neighboring sites, one can similarly show that D ∼ 〈k〉 also for the case of jumps to neighbors further away (cf. ESI,‡ Section E). For diffusive motion, mobility and diffusivity are related according to Einstein–Smoluchowski by^58,59

eD = μk_BT (7)

Hence the mobility follows μ ∼ 〈k〉/T, with 〈k〉 being the average transfer probability. In the high temperature limit, the full rate depends linearly on T for both upward and downward jumps (hence, 〈k_full〉 ∼ T), while the approximate rate is constant or converges to a constant value for downward and upward jumps, respectively (implying 〈k_approx〉 ∼ const.). This is illustrated in Fig. 5(b) for ΔE = 5 meV and in the ESI‡ for larger values of ΔE (Fig. S.7, ESI‡). With this, one obtains for the temperature dependence of the mobility:

Hence, for small disorders, μ_approx decreases with increasing temperature, while μ_full is approximated as temperature independent.

Thus, overall, ln μ ∼ (σ/k_BT)² applies for the full and the approximate rate in the limit of high disorders, while eqn (8) and (9) represent the limit of low disorders. In between we expect a continuous transition in the evolution of μ(σ,T). Evidently, μ_full scales with σ/k_BT, so that it does not matter if a KMC simulation was conducted with, say, half the value of σ or twice the value of T, and curves obtained by varying the one or the other all coincide (cf.Fig. 3). Similarly, other curves obtained for different values of σ, yet same σ/k_BT also coincide. Finally, the slope of the ln μ_full ∼ (σ/k_BT)² curve is −(2/3)² for all values of (σ/k_BT)². We recall that the slope is indicated by the grey line in Fig. 3b, with the deviation of the simulation from this for large values of (σ/k_BT)² being due to the energy relaxation of charge carriers not reaching equilibrium. A scaling of mobility with (σ/k_BT)² is not the case for μ_approx due to the 1/T dependence (eqn (9)) which becomes critical when performing KMC simulations for low disorders. The relationship has been established on the basis of analytical calculations and kinetic Monte Carlo simulations with the approximate MA rate. It was found to agree very well with experiments for several decades and has thus become an established reference for the description of charge transport in systems with moderate to high disorder and negligible polaronic contributions.⁵⁴ The fact that the use of the full MA rate preserves this relationship also for the regime of low disorders is therefore reassuring.

The observations made by Chatten et al. and Fishchuk et al. can thus be understood in the context of eqn (8) and (9).^48,50 There are two points in the work by Fishchuk et al. that deserve further discussion. First, Fishchuk does not find a universal slope of (2/3)² ≈ (0.66)² for the ln μ_full ∼ (σ/k_BT)² curve for all values of (σ/k_BT)² like we do. Instead, he obtains two limiting slopes of (4/9)² ≈ (0.44)² and (3/5)² = (0.60)² for small and large values of σ at constant temperature, which might be due to limitations of his effective-medium-approximation approach. For his calculations, Fishchuk was using an effective-medium theory in combination with the concept of an effective transport level. Evidently, this is a good approximation to the transport situation for well disordered systems, yet the effective transport level concept becomes progressively less suitable for more ordered semiconductors.

The second point worth mentioning is that, over the entire parameter range, Fishchuk et al. find μ_approx to have a stronger dependence on (σ/k_BT)² than μ_full, while we find the opposite. More specifically, they noted a stronger temperature dependence for μ_approx, even though his mobility values for the temperature dependence were not obtained by variating the temperature but by variating the disorder. The approximate rate does not scale with σ/k_BT but depends differently on both parameters. Consequently, it is not possible to directly compare the temperature dependencies. μ_approx can exhibit a stronger or weaker temperature dependence compared to μ_full, depending on the disorder considered.

Trends in the dependence on electric field

We now turn to the field dependence of the mobility, shown in Fig. 4. In summary the mobility derived from using the full rate, μ_full, increases monotonously with field, following for not too low fields, and eventually converges towards a maximum mobility. When using the approximate rate, the resulting μ_approx shows a similar increase and values for large disorders and low fields, yet at high fields or low disorders, μ_approx decreases with increasing field following a μ ∼ 1/F dependence. The dependence is experimentally well verified and extends also to lower fields when using correlated disorder (cf. ESI,‡ Fig. S.3) instead of the non-correlated disorder we employed here.^28,60 As discussed by Poole and Frenkel, it results from the lowering of the energy barrier by the electric field.¹ In this regime of moderate fields, and hence moderate ΔE, the transport is predominantly limited by the uphill jumps. Since the full and approximate uphill rate converge with increasing ΔE (cf.Fig. 5), it is not surprising that the electric field has the same effect on transport simulated with both rates.

More remarkable is the difference in the mobility for high fields. The reduction of mobility with increasing field has occasionally been observed experimentally^61–65 albeit mostly for field strengths of 0.01–1.44 MV cm⁻¹. It was phenomenologically explained by the notion that higher fields would prevent charge carriers from circumnavigating a particularly high energy barrier by a sequence of low barrier jumps.²⁴ This notion was supported by the KMC simulations based on the approximate MA rate that explicitly considers positional disorder through a parameter Σ that modified the prefactor of the jump rate. However, even when the positional disorder was low (Σ < 1.5), μ(F) was found to reduce with the field. To obtain agreement between the results of the KMC simulations and the analytical expression, an empirical, additional constant of 2.25 was introduced.^24,66 The field dependence of μ_approx was then described as

with C̃ being a constant that depends on the lattice constant a.²⁴ The simulations carried out here with the full rate yield a constant mobility for all disorders at high fields. In light of the full-rate simulations, the ad hoc introduction of the constant factor of 2.25 for low positional disorder (eqn (10.ii)) seems no longer justified. Obviously, the role of positional disorder in modifying the coupling between sites is an aspect requiring further attention. To which extent eqn (10.i) correctly describes charge transport in the presence of high fields and positional disorder is an aspect requiring further attention that is beyond the scope of the current manuscript. Moreover, experimental data for field strengths above 1 MV cm⁻¹ for films with low disorder are required, which are at present not available.

What causes the different field dependences when using the full or the approximate rate? Let us consider the limit of high fields (eaF ≥ σ, where a is the lattice constant). For σ = 50 meV, a high field would be F ≥ 0.5 MV cm⁻¹, i.e. MV^1/2 cm^−1/2 (cf. ESI‡ Section G). In this regime, charge carrier motion by diffusion is negligible and, instead, drift, i.e. directed transport driven by the electric field, predominates. The field lowers the energy barrier between sites in field direction, so that in the extreme case virtually all hops are downhill and directed towards the collecting electrode. Hence, the drift velocity v_Drift is directly proportional to the average jump rate 〈k〉, i.e. v_Drift ∼ 〈k〉. Moreover, upward jumps cease to exist in this extreme limit, so that only the downward jumps are rate-limiting.

For example, when a field of 2 MV cm⁻¹ acts on the sample, the energy between two adjacent sites in the field direction is reduced by 200 meV, so that any former disorder-related energy barrier with less than 200 meV vanishes. For large downward jumps, the full rate accelerates nearly linearly with the energy gap between adjacent sites, ΔE, while the approximate rate is constant (see Fig. 5(a)). Furthermore, on average, ΔE can be assumed to be proportional to the field F for large fields. Therefore, one finds for the average rates at high fields 〈k_full〉 ∼ F and 〈k_approx〉 = const.; and with the definition of the drift mobility, μ = v_Drift/F, accordingly,

Therefore, in the regime where directed transport towards the drain electrode predominates, (eaF ≥ σ), both rates lead to different dependencies of the mobility on the field. Eqn (12), i.e. the decrease of mobility with the field in the approximated MA rate arises from the implicit limited phonon density of states. Eqn (11) and (12) agree with the analytical results of Fishchuk et al. in the large field limit.⁵⁰ We should point out a caveat of our treatment. We did not investigate the effect of positional disorder, which modifies the coupling between sites and might possibly reduce the difference between the results from full and approximate rates. This aspect needs to be addressed in future work.

Finally, we briefly comment on the downward rate in both the approximate and full expressions from Miller and Abrahams. As noted by Fornari et al., a constant as well as an increasing rate with increasing energy gap may seem at variance with the energy gap law. The energy gap law is well established and predicts decreasing rates for increasing downward energy gaps.^49,55 However, in the energy gap law, the transfer of energy is limited by the dissipation of excess energy through a multiphonon emission process. This becomes increasingly less likely the more phonons of the same energy need to be emitted. The downhill (or uphill) charge transfer process considered in the Miller–Abrahams model, in contrast, involves the emission (or absorption) of a single phonon. The rate is limited by the electronic coupling between the two states rather than by the dissipation of energy. As detailed by Fornari et al., the emission (or absorption) of a phonon not only satisfies energy conservation requirements but also modulates the electronic coupling between the two states involved. Hence, the rate rises with the availability of phonon states.

Conclusions

The choice of hopping rate is a critical question when modelling charge transport in organic semiconductors.²¹ Inspired by the analytical treatment from Fishchuk et al., we highlighted that the hopping rate which is conventionally referred to as Miller–Abrahams rate is an approximation of a more general (full) rate. This approximation holds for the case where the site energy difference is large compared to thermal energies.⁵⁰ Using KMC simulations, we illustrate how the common, approximated rate impacts notably μ(T) and μ(F). Mathematically, the approximation is valid for |ΔE| ≫ 2k_BT. Investigating the temperature and field dependence of the charge carrier mobility with the full and the approximated rate suggests that the condition σ ≥ k_BT is already sufficient to obtain μ_approx(T) ≈ μ_full(T). However, we point out that μ_approx(F) ≈ μ_full(F) additionally requires eaF < σ. In simple words, the approximate rate works well for sufficiently large disorders and low to moderate electric fields.

We analyzed the differences between rates and resulting charge carrier mobilities for the full and the approximated equations as a function of the energy gap between adjacent sites and as a function of temperature. In essence, the differences in the mobilities obtained from both rates may be attributed to different implicit shapes of the spectral density J(ω), with J(ω) saturating for large energy gaps for the approximate expression, while it increases for the full expression. J(ω) may be seen as a measure for the coupling strength per phonon mode.

Previous discussions on the deficiencies of the MA rate considered only the approximate rate.^{17,20,41–48,67} Some of the issues raised resolve themselves when the full rate is considered in a parameter range where the approximated rate is no longer valid. The use of the full rate becomes important when quantitatively analyzing the mobilities in systems with weak disorder, such as highly aligned or crystalline material. In fact, models that describe charge transport in very weakly disordered systems through the concept of transient (de)localization focus strongly on the explicit role of phonon modes in the charge transport.^10,13,68 Using the full MA rate may help to conceptually link the description of charge transport in disordered and in well-ordered organic semiconductors.

Author contributions

All authors together conceived and conceptualized the research project. MSD developed and wrote the software, performed all simulations, formally analyzed them and visualized the results. She also wrote the first draft of the manuscript. All authors critically discussed the results together and edited and reviewed subsequent versions of the manuscript. AKö acquired funding, provided resources and supervised the project.

Data availability

The data that support this article are available from the corresponding author upon reasonable request.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

AKö and HO acknowledge financial support from the Bavarian State Ministry for Science and Arts through the Initiative “Solar Technologies go Hybrid (SolTech)”, the German Science Foundation through the international research training group IRTG2818 “OPTEXC”, and the German Science Foundation through the CRC 1585 MultiTrans. HO further thanks the German Science Foundation for a Heisenberg professorship (OB42519-7). MSD was supported by the elite network of Bavaria (Study Program Biological Physics) and the University of Bayreuth Graduate School. AKa also thanks the IRTG2818 “OPTEXC”.

References

1. A. Köhler and H. Bässler, Electronic Processes in Organic Semiconductors: An Introduction, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2015.
2. H. T. Nicolai, M. Kuik, G. A. H. Wetzelaer, B. de Boer, C. Campbell, C. Risko, J. L. Brédas and P. W. M. Blom, Nat. Mater., 2012, 11, 882–887.
3. H. Sirringhaus, Adv. Mater., 2014, 26, 1319–1335.
4. G.-D. Ye, R. Ding, S.-H. Li, L. Ni, S.-T. Dai, N.-K. Chen, Y.-F. Liu, R. Guo, L. Wang, X.-B. Li, B. Xu and J. Feng, Light: Sci. Appl., 2024, 13, 136.
5. S. Fratini, M. Nikolka, A. Salleo, G. Schweicher and H. Sirringhaus, Nat. Mater., 2020, 19, 491–502.
6. O. Yildiz, Z. Wang, M. Brzezinski, S. Wang, Z. Li, J. J. Michels, P. W. M. Blom, W. Pisula and T. Marszalek, Adv. Funct. Mater., 2024, 34, 2314131.
7. H. Sirringhaus, P. J. Brown, R. H. Friend, M. M. Nielsen, K. Bechgaard, B. M. W. Langeveld-Voss, A. J. H. Spiering, R. A. J. Janssen, E. W. Meijer, P. Herwig and D. M. Leeuw, Nature, 1999, 401, 685–688.
8. A. L. T. Khan, P. Sreearunothai, L. M. Herz, M. J. Banach and A. Köhler, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 085201.
9. M. Reichenberger, D. Kroh, G. M. M. Matrone, K. Schötz, S. Pröller, O. Filonik, M. E. Thordardottir, E. M. Herzig, H. Bässler, N. Stingelin and A. Köhler, J. Polym. Sci., Part B: Polym. Phys., 2018, 56, 532–542.
10. S. Giannini, A. Carof, M. Ellis, H. Yang, O. G. Ziogos, S. Ghosh and J. Blumberger, Nat. Commun., 2019, 10, 3843.
11. S. Giannini and J. Blumberger, Acc. Chem. Res., 2022, 55, 819–830.
12. R. A. Marcus, Angew. Chem., Int. Ed. Engl., 1993, 32, 1111–1121.
13. A. Troisi and G. Orlandi, Phys. Rev. Lett., 2006, 96, 086601.
14. L. B. Schein, C. B. Duke and A. R. McGhie, Phys. Rev. Lett., 1978, 40, 197–200.
15. K. H. Probst and N. Karl, Phys. Status Solidi A, 1975, 27, 499–508.
16. N. Tessler, Y. Preezant, N. Rappaport and Y. Roichman, Adv. Mater., 2009, 21, 2741–2761.
17. Z. G. Yu, D. L. Smith, A. Saxena, R. L. Martin and A. R. Bishop, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 085202.
18. M. P. A. Fisher and A. T. Dorsey, Phys. Rev. Lett., 1985, 54, 1609–1612.
19. K. Asadi, A. J. Kronemeijer, T. Cramer, L. J. A. Koster, P. W. M. Blom and D. M. de Leeuw, Nat. Commun., 2013, 4, 1710.
20. N. Vukmirović and L.-W. Wang, Appl. Phys. Lett., 2010, 97, 043305.
21. K. Zojer, Adv. Opt. Mater., 2021, 9, 2100219.
22. A. Miller and E. Abrahams, Phys. Rev., 1960, 120, 745–755.
23. V. Ambegaokar, B. I. Halperin and J. S. Langer, Phys. Rev. B, 1971, 4, 2612–2620.
24. H. Bässler, Phys. Status Solidi B, 1993, 175, 15–56.
25. J. Lee, J. Phys. Appl. Phys., 2021, 54, 143002.
26. H. Li and J.-L. Brédas, Natl. Sci. Rev., 2021, 8, nwaa167.
27. H. Bässler, Phys. Status Solidi B, 1981, 107, 9–54.
28. Yu. N. Gartstein and E. M. Conwell, Chem. Phys. Lett., 1995, 245, 351–358.
29. M. Bouhassoune, S. L. M. van Mensfoort, P. A. Bobbert and R. Coehoorn, Org. Electron., 2009, 10, 437–445.
30. W. F. Pasveer, J. Cottaar, C. Tanase, R. Coehoorn, P. A. Bobbert, P. W. M. Blom, D. M. de Leeuw and M. A. J. Michels, Phys. Rev. Lett., 2005, 94, 206601.
31. I. I. Fishchuk, A. Kadashchuk, S. T. Hoffmann, S. Athanasopoulos, J. Genoe, H. Bässler and A. Köhler, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 125202.
32. J. Cottaar, L. J. A. Koster, R. Coehoorn and P. A. Bobbert, Phys. Rev. Lett., 2011, 107, 136601.
33. J. Cottaar, R. Coehoorn and P. A. Bobbert, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 245205.
34. P. Muralidharan, S. M. Goodnick and D. Vasileska, J. Comput. Electron., 2019, 18, 1152–1161.
35. M. Baranowski, J. M. Urban, N. Zhang, A. Surrente, D. K. Maude, Z. Andaji-Garmaroudi, S. D. Stranks and P. Plochocka, J. Phys. Chem. C, 2018, 122, 17473–17480.
36. S. Nokhrin, M. Baru and J. S. Lee, Nanotechnology, 2007, 18, 095205.
37. S. Banerjee, R. Biswas, K. Seki and B. Bagchi, J. Chem. Phys., 2014, 141, 124105.
38. P. K. Watkins, A. B. Walker and G. L. B. Verschoor, Nano Lett., 2005, 5, 1814–1818.
39. M. Onizhuk, S. Sohoni, G. Galli and G. S. Engel, J. Phys. Chem. Lett., 2021, 12, 6967–6973.
40. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. Phys., 1953, 21, 1087–1092.
41. H. Houili, E. Tutiš, I. Batistić and L. Zuppiroli, J. Appl. Phys., 2006, 100, 033702.
42. Ya. V. Burdakov and V. R. Nikitenko, J. Phys.: Conf. Ser., 2017, 938, 012067.
43. J. Stephan, S. Schrader and L. Brehmer, Synth. Met., 2000, 111–112, 353–357.
44. N. J. van der Kaap, I. Katsouras, K. Asadi, P. W. M. Blom, L. J. A. Koster and D. M. de Leeuw, Phys. Rev. B, 2016, 93, 140206.
45. A. V. Nenashev, F. Gebhard and S. D. Baranovskii, Phys. Rev. B, 2020, 102, 066201.
46. S. T. Hoffmann, S. Athanasopoulos, D. Beljonne, H. Bässler and A. Köhler, J. Phys. Chem. C, 2012, 116, 16371–16383.
47. I. Katsouras, A. Najafi, K. Asadi, A. J. Kronemeijer, A. J. Oostra, L. J. A. Koster, D. M. de Leeuw and P. W. M. Blom, Org. Electron., 2013, 14, 1591–1596.
48. A. J. Chatten, S. M. Tuladhar and J. Nelson, in Conference Record of the Thirty-first IEEE Photovoltaic Specialists Conference, 2005, 2005, pp. 31–36.
49. R. P. Fornari, J. Aragó and A. Troisi, J. Chem. Phys., 2015, 142, 184105.
50. I. I. Fishchuk, D. Hertel, H. Bässler and A. K. Kadashchuk, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 125201.
51. H. Li, Y. Li, H. Li and J.-L. Brédas, Adv. Funct. Mater., 2017, 27, 1605715.
52. H. Li and J.-L. Brédas, J. Phys. Chem. Lett., 2017, 8, 2507–2512.
53. R. Saxena, V. R. Nikitenko, I. I. Fishchuk, Ya. V. Burdakov, Yu. V. Metel, J. Genoe, H. Bässler, A. Köhler and A. Kadashchuk, Phys. Rev. B, 2021, 103, 165202.
54. P. M. Borsenberger and D. S. Weiss, Organic Photoreceptors in Xerography, Marcel Dekker, Inc., New York, 1998, vol. 49.
55. R. Englman and J. Jortner, Mol. Phys., 1970, 18, 145–164.
56. S. Heun and P. M. Borsenberger, Chem. Phys., 1995, 200, 245–255.
57. R. LeSar, Introduction to Computational Materials Science: Fundamentals to Applications, Cambridge University Press, Cambridge, 2013.
58. M. von Smoluchowski, Ann. Phys., 1906, 326, 756–780.
59. A. Einstein, Ann. Phys., 1905, 322, 549–560.
60. S. V. Novikov, D. H. Dunlap, V. M. Kenkre, P. E. Parris and A. V. Vannikov, Phys. Rev. Lett., 1998, 81, 4472–4475.
61. A. Peled and L. B. Schein, Chem. Phys. Lett., 1988, 153, 422–424.
62. A. J. Mozer and N. S. Sariciftci, Chem. Phys. Lett., 2004, 389, 438–442.
63. A. J. Mozer, N. S. Sariciftci, A. Pivrikas, R. Österbacka, G. Juška, L. Brassat and H. Bässler, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 035214.
64. V. Kažukauskas, M. Pranaitis, V. Čyras, L. Sicot and F. Kajzar, Thin Solid Films, 2008, 516, 8988–8992.
65. M. Novo, M. van der Auweraer, F. C. de Schryver, P. Borsenberger and H. Bässler, Phys. Status Solidi B, 1993, 177, 223–241.
66. H. Bässler and A. Köhler, in Unimolecular and Supramolecular Electronics I: Chemistry and Physics Meet at Metal-Molecule Interfaces, ed. R. M. Metzger, Springer, Berlin, Heidelberg, 2012, pp. 1–65.
67. X. de Vries, P. Friederich, W. Wenzel, R. Coehoorn and P. A. Bobbert, Phys. Rev. B, 2018, 97, 075203.
68. S. Fratini and S. Ciuchi, Phys. Rev. Lett., 2009, 103, 266601.

---

Footnotes

† In memoriam I. I. Fishchuk (1945–2024).

‡ Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5tc01487e

---