Extremely Non-Equilibrium Hopping Transport and Photogeneration Efficiency in Organic Semiconductors: An Analytic Approach

An analytical model of highly nonequilibrium hopping transport of charge carriers in disordered organic semiconductors has been developed. In particular, the initial time interval is considered when transport is controlled by hops down in energy. The model is applied to the calculation of the separation probability of geminate pairs in a semiconductor with a Gaussian energy distribution of localized states. This probability determines the photogeneration efficiency. The temperature dependence of the separation probability is obtained and shown to be much weaker than predicted by the classical Onsager model, in agreement with experiment and Monte Carlo simulations. The field dependence is taken into account using a modified effective temperature method.

−4 Despite the enormous diversity of structure and optoelectronic characteristics, organic semiconductors, including both low-molecular materials and polymers, share common transport properties, such as the hopping nature of transport of charge carriers and excitations, as well as structural and associated energy disorder.It is the latter that determines the strong temperature and field dependence of mobility in disordered organic semiconductors, according to the Gaussian disorder model. 5Within the framework of this model, many groups have carried out theoretical modeling of the mobility and diffusion coefficient of charge carriers, both by the Monte Carlo method 5−7 and by analytic methods; see, for example, refs 8−11.It is known that energy disorder leads to long thermalization times for the initially nonequilibrium energy distribution of carriers, for example, in the case of photogeneration.Thermalization occurs during transport, and this leads to anomalous transport characteristics, such as dispersive or extremely nonequilibrium transport. 5,8,12However, most studies considered the case of either quasi-equilibrium transport 9−11 or the case when highly nonequilibrium transport is controlled by thermally activated jumps.In these cases, the drift shift and diffusion spreading of the carriers in a weak field are connected by the Einstein relation. 13−16 The mode of jumps down in energy should be taken into account when analyzing the efficiency of photogeneration of charge carriers. 17,18There is a strong argument that the primary photoexcitations in organic semiconductors, including conjugated polymers, are molecular excitations (excitons). 19,20citon decay leads to the formation of electron−hole pairs coupled by Coulomb interaction (geminate pairs).Since the charges ("twins") must separate by a distance greater than r c (Onsager radius, r c > 13 nm, at temperature T < 300 K) which significantly exceeds the typical hopping length (on the order of a nanometer), the probability of pair separation, which governs the quantum charge photogeneration yield, is determined by the diffusion-drift motion of the "twins" in their Coulomb field and an external uniform electric field.Thus, the separation of geminate pairs is controlled by the transport of charge carriers.Transport can occur in an extremely nonequilibrium regime, especially at early times after a photoexcitation and at low temperatures (thermal energy is much smaller than the scale of energy disorder).However, the effect of energy disorder on the photogeneration of charge carriers has been studied much less than on transport in a uniform electric field.Often theoretical analysis is carried out within the framework of the classical Onsager model, which does not take into account energy disorder and the hopping nature of transport. 21,22Monte Carlo simulation 17 showed that the temperature dependence of mobility in a weak external field at low temperatures does not agree with this model, similarly to the experimental data. 23Even if the Onsager model is qualitatively correct (the temperature dependence is described by the Arrhenius law at sufficiently high temperatures), the question arises about the value of the key parameter of the model�the initial distance between the "twins" (initial separation), r 0 .Monte Carlo simulations showed that the effective (apparent) magnitude of the initial separation increases with decreasing temperature and is determined by the initial energy relaxation; i.e., downward energy jumps dominate. 17,18n a recent work, 18 we showed, using the model of the exponential energy distribution of hopping centers, that the analytic model of extremely nonequilibrium transport gives a weak temperature dependence in agreement with the Monte Carlo data 17 and experiment. 23In comparison with Monte Carlo, the analytic approach is not associated with specific parameter values and facilitates the analysis of the contribution of various processes to the phenomenon under study.It was shown that anomalously strong diffusion during the initial time interval (after pair generation) leads to an increase in the effective initial separation, which is greater the lower the temperature.However, simplifying assumptions were made in the work. 18In particular, the exponential energy distribution of hopping centers was considered, while organic semiconductors are characterized by a Gaussian distribution.
In this work, an analytic model is constructed for an arbitrary (rapidly decreasing in depth, for example, Gaussian) energy distribution of states (DOS), g(E), and it is shown that the relationship between drift and diffusion in the initial time interval and at low temperatures is determined by the type and width of the DOS, and at large times there is a transition to the Einstein relation.A characteristic time has been found after which the Einstein relation is valid and the Onsager model is applicable.An analytic model of drift and diffusion in an extremely nonequilibrium regime is applied to the problem of the separation probability of a geminate pair.Both temperature and field dependences are analyzed using the effective temperature method modified in this work.From the Master Equation to Equations of the Multiple Trapping and Release (MTR) Model.The left panel of Figure 1 shows the workflow of the derivation.We start from the well-known master equation for hopping transport 24 where f i ≪ 1 and f j ≪ 1 are occupation probabilities of states i and j and ν ij and ν ji are rates of hopping from i to j and vice versa, as defined by the Miller−Abraham (MA) model.We describe the localized states by their energy distribution (DOS function), g(E), and by means of this DOS, we rewrite the equation for the energy distribution of occupied states (ODOS), ρ(E, x, t) = g(E)f(E, x, t), passing from summation to integration where x is a macroscopic coordinate, r is a (microscopic) hopping distance, e is the elementary charge, s = cos θ, and θ is the angle between the direction of electric force, −eF ⃗ , and the direction of a jump, r⃗ .Equation 2 takes into account that the electron energy includes the electrostatic energy of the external field; see the right panel of Figure 1.
The concentration of charge carriers (electrons) is The rate of escape from a state with energy E is Unlike the hopping rates in eq 1, the hopping rate ω(E, E′, r) takes into account the possibility of going to states other than the final state, as well as the probability that there will be no return to the initial state; see the Supporting Information.For the small-field limit, one can rewrite eq 4 as At first, we consider eq 2 in the zero approximation, assuming F = 0 and neglecting the x-dependence of the function ρ(E, t), in order to analyze the energy relaxation itself.Equation 2 reduces to

The Journal of Physical Chemistry Letters
Considering sufficiently a nonequilibrium (not quasi-equilibrium) process, one has to note that states with a given energy E are predominantly filled by carriers from states with energies E′, such that the condition ω(E′)t ≫ 1 is satisfied (i.e., from "currently shallow" states).The time-dependent DOS for these states is g sh (E′, t) = g(E′)[1 − exp[−ω(E′)t]], according to Poisson's distribution. 25Oppositely, one can neglect transitions from the "currently deep" states; their distribution is Qualitatively, one can define the "currently shallow" and "currently deep" states by the conditions E′ ≥ E d (t) and E′ < E d (t), respectively.The demarcation energy, E d (t), is defined by the condition The detailed balance principle and expressing the ODOS in a separable form 10,16 as allow us to modify eq 6 to One can write a nonequilibrium energy distribution, ρ(E, x, t), in the form of eq 9 or 10, considering that (1) the length scale of variations with x considerably exceeds the hopping length; (2) , because the rate of downward jumps does not depend on the final energy (Miller−Abrahams model).
The release of a carrier from a state of the energy E by jumps downward and upward in energy prevails under the conditions E > E 1/2 and E < E 1/2 , respectively; 13,16,18 see part S1 of the Supporting Information (SI) about the energy E 1/2 .After the "segregation time" t s , which is defined by the condition E d (t s ) = E 1/2 , transport is controlled by thermally activated jumps, while downward (in energy) jumps prevail earlier. 8,13,26Taking the limit E d (t) → −∞ in the integral in eq 11 at t ≫ t s , using eqs 5, 10, and 11 and expressing one obtains or where c(E) = ω(E) exp(−E/kT)/N c .Equation 14is the balance equation of the MTR model. 27,28N c and n c (x, t) are the concentrations of "conductive" states (i.e., the states giving the main contribution to transport) and charge carriers in these states ("mobile carriers"), respectively.Equation 14 has been derived in the work 11 for the quasi-equilibrium conditions.
Returning to eq 2, we expand its right-hand side in a power series in small parameters eFrs and rs up to the second order, assuming that only carriers from "shallow" states can contribute to the transport process.Keeping only the first power of electric field (assuming weak field), we obtain after the integration over energy 2 sh diff dr (15)   where where a is the typical hopping distance; see below.
Using eq 9, one can rewrite the first (diffusion) term in the right-hand side of eq 15 as where is the concentration of "mobile" carriers and τ 0 is the typical escape time from "conductive states".Indeed, eq 18 follows from the conditions that the escape time from a state with the given energy E′ is smaller than τ 0 and the time τ 0 is rather small, ω(E′)τ 0 ≪ 1; hence the Poisson's probability that the state is "conductive" is approximately equal to ω(E′)τ 0 : The Journal of Physical Chemistry Letters i.e., it is a constant at t ≫ t s if the function g(E′) decreases rapidly with energy E′.Integration in parts in the second (drift) term in eq 15, using eqs 8−10, gives Thus, one can express eq 15 in the MTR model form where D c = a(t) 2 /(6τ 0 ) is the coefficient of diffusion of "mobile carriers".Equations 3, 14, and 22 form a complete set of MTR model equations.Since only a minority of charge carriers can be classified as "mobile" (ω(E′)τ 0 ≪ 1), their concentration is neglected in eq 3.One has to note that the value of τ 0 disappears from eqs 15 and 22; see eqs 18 and 19.
Strictly speaking, eq 14 is valid in a limited time interval, t ≫ t s ; however, one can consider also the initial time interval (t ≪ t s ) in the extremely nonequilibrium transport regime, assuming that the majority of carriers occupy "currently deep" states, i.e., n(x, t) ≈ f 0 (x, t)∫ dEg d (E, t)f F (E, t).For the "currently deep" states, one can neglect the second (release) term in eqs 11 and 14: Integration of eq 23 over energy E and then over time gives the relationship known in the theory of extremely nonequilibrium transport, 13,25,30 = n x t t t n x t ( , ) ( ) d ( , ) i.e., where The difference from the usual expression 13,25,30 is only in the appearance of the function Ω(t, E).In deriving the above equations, the adiabatic approximation was used, 13,25 according to which the function n c (x, t′) changes much faster than the function g d (E, t), due to fast hopping between "conductive" states.The approximation of "extremely non-equilibrium transport" is valid under the condition of rather strong disorder, kT ≪ E 0 , where E 0 is the energy scale of the DOS g(E), at a limited initial time interval, t ≪ t eq .For example, E 0 = 2 σ, t eq is defined by the condition 5 E d (t eq ) = −σ 2 /kT for the Gaussian function g(E) with the variance σ.
To obtain the equation for the total concentration of charge carriers, n(x, t), we analyze the time dependence of the characteristic energy E * (t) in the realistic case kT ≪ E 0 ; see eq in the first term of eq 21 and using eqs 5 and 23, one can express this term as (1/kT)[1 − τ 0 /τ(t)].In the long-time limit, t ≫ t s , g sh (E, t) ≈ g(E), except for a deep tail.Hence, the first term of eq 21 approaches 1/kT, while the second term approaches zero with the concentration of "currently deep" states.Vice versa, the second term of eq 21 prevails at the initial time period, t ≪ t s , and downward jumps from the narrow layer near the demarcation energy E d (t) control the transport. 14,16Hence, the typical release rate is ω(E d ) = t −1 .An estimation gives that the first term is small as kT/E 0 .Equation 9 approaches the step function, and one can estimate f 0 (x, t) ≈ n(x, t)/G[E d (t)] in eq 9, where G(E) ≡ ∫ −∞ E dEg(E).Since downward jumps prevail, i.e., E′ > E, . Estimating the second term in eq 21, multiplied by ∂n c (x, t)/∂x, and using eq 9, one obtains since one can estimate ∂f 0 (x, t)/∂t ≈ f 0 (x, t)ω[E d (t)] = f 0 (x, t)/t. 16,18ombining the above estimations and eqs 22 and 26, one obtains Integration over time gives The Journal of Physical Chemistry Letters where D c (t) = a 2 (t)/(6τ 0 ), and the energy E * (t) is estimated as Since the majority of jumps occur from "currently shallow" to "currently deep" states, ∫ dEg d (E, t)ω ̅ (E′, E) ≈ ω(E′) in eqs 24 and 27, in accordance with eq 5. Hence, In the opposite case, t ≫ t s , ω(E′) = ν 0 exp[−(E tr − E′)/kT] for the majority (E ≪ E 1/2 ) of occupied "currently shallow" states, 10,13  the same as in the usual MTR approach in the dispersive transport regime, where N t = G[E 1/2 ] is the total concentration of "traps", i. e., states below the transport energy E 1/2 , that is an analog of the mobility edge of the MTR model. 13eparation Probability of Geminate Pairs.On the basis of the above results, we can write a Smoluchowski equation for the dispersive transport regime of "geminies" 18,31 For the case of a negligible external field (the case of central symmetry, φ(r) = −e/4πεε 0 r) eq 35 takes the form ( ) provided that n(r, 0) = δ(r − r i )/4πr i 2 ; r i is the distance between the "twins" right after the exciton dissociation ("zero jump").
Solution of this equation, n(r, t), is obtained in part S2 of the SI for the case of a negligible external field using the Wentzel− Kramers−Brillouin (WKB) approximation.
In fact, an initial energy relaxation at t < t s is what is considered in the conventional Onsager model to be the thermalization of "hot" carriers; just this process results in an initial carrier separation at a distance r 0 .Then, at t > t s , one can apply the Onsager model.After the initial relaxation, E * (t) = kT = const, and transport kinetics (dispersive or quasiequilibrium) does not affect the separation probability, 31 Ω ∞ .We use the spatial distribution of carriers at time t s , n(r 0 , t s ) (see the SI), as the distribution over the true initial separations 18  In eq 37 Ω ∞ Ons (r 0 ) = exp(−r c /r 0 ), 21 where r c = e 2 /(4πεε 0 kT) is the Onsager radius (Coulomb radius) and ε is the relative dielectric permittivity.Examples.Equation 30 differs from the respective equation of the conventional MTR model 13,27 because of the following: (1) time dependences of functions τ(t), D c (t) are different at t ≪ t s (previously these equations were considered only for not too small times, t ≫ t s ); (2) the energy E * (t) in the drift term is time-dependent and differs from the (constant) thermal energy, kT.Since τ(t) is an increasing function of time, E * (t) → kT in eq 31 in the long-time limit, t > t s , while the first term dominates at short times, t ≪ t s : It is temperature-independent but time-dependent, hence violating the Einstein relation (one has to note, however, that the form of eq 30 differs from the conventional drift-diffusion equation).Equation 30 results from the approximation of extremely nonequilibrium transport; i.e., the majority of carriers occupy "currently deep" states. 25,30This approximation is valid at the limited time interval, t ≪ t eq , with t eq being the time of quasi-equilibration. 5 For the case of Gaussian DOS, Both times t s and t eq increase sharply with decreasing temperature, but the relation t s ≪ t eq is fulfilled; see Figure 2.

The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter
The Einstein relation is fulfilled if the time is rather long, t > t s ; i.e., transport is controlled by upward jumps from "currently shallow" (i.e., quasi-equilibrium) states to the states near the transport level E 1/2 , that contributes mostly to the transport. 13ppositely, i.e., if t ≪ t s , the thermal energy kT is replaced in eqs 30 and 36 by the temperature-independent but timedependent energy, This circumstance, i.e., the ratio of diffusive dispersion and drift shift scaled by energy E 0 instead of kT, was argued previously 14,15 for the case of exponential DOS, g(E) = (N/E 0 ) exp[E/E 0 ], E < 0, assuming that T = 0.The initial time period, t ≪ t s , is similar to the T = 0 case in the sense that jumps down in energy predominate.Our analysis shows that, for other (nonexponential) DOS, the energy E * depends on time even at the very initial time interval; see Figure 3. Figure 3 shows that the time interval when E * ≫ kT is rather short relative to the segregation time, t s , although it is long relative to the hopping time, t hop = ν 0 −1 exp[−γ(6/πN) 1/3 ], if T < 300 K; E * ≈ kT at t ≥ t s .The demarcation energy, E d (t), also approaches its asymptotic form, if t ≥ t s ; see Figure 4.The duration of the "anomalous" initial time period sharply increases with decreasing temperature.
In the work, 18 assuming an exponential DOS, the following approximation was used for the function E * (t) (combination of asymptotes): E * (t) = E 0 , t ≤ t s ; E * (t) = kT, t > t s .That is in apparent contradiction with the results of calculations for the Gaussian DOS; see Figure 3.However, calculations of the temperature dependence of separation probability (see eq 37) are in qualitative agreement with the results of the work. 18amely, the separation probability is much larger than the prediction of the Onsager model at low temperatures; see Figure 5.
The separation probability depends on hopping parameters.Namely, it increases with increasing localization parameter, loc = 2γN −1/3 , and with increasing disorder, σ, at a given temperature.One can introduce an effective value of the initial pair separation, 17,18 r 0 ef .Hence, one can calculate the geminate recombination quantum yield as Ω ∞ (T) = exp(−e 2 / 4πεε 0 kTr 0 ef (T)).This expression does not take into account an external field.Meanwhile, an external field is usually present in both Monte Carlo simulations 17 and experiments. 23he effective temperature concept (ETC), modified in this work, can be used for a joint analysis of the temperature and field dependence of the separation probability.The concept of field-dependent effective temperature T eff was proposed initially by Shklovskii, 32 and it was repeatedly applied later to various disordered materials; see, for example, refs 33 and 34.The essence of ETC is that the electric field F reduces the activation energy of jumps by an amount eFl, increasing their frequency and, accordingly, mobility.In a disordered medium, in the case of a uniform field, the length scale, l, is the localization radius, γ −1 .An increase in temperature has the same qualitative effect.In the limiting case of T = 0 one has T eff (0, F) = λeFl/k, where λ is a numerical constant.In general, interpolation expressions such as are used.Usually, n = 2, λ = 0.67. 33,34The application of the ETC suggests that the quantum yield of geminate recombination can be calculated as The common ETC assumes that the applied field is uniform and release from the initial state is a one-step process.However, in reality, separation of a geminate pair is a multistep process, and the electric field is nonuniform.The carrier

The Journal of Physical Chemistry Letters
undergoes a diffusion−drift motion in the Coulomb potential well.Therefore, one should not expect quantitative results if expression 39 is used in eq 40.
At the same time, the ETC provides a qualitatively correct dependence Ω ∞ (r 0 , T, F) in a practically significant interval, 1 nm < r 0 eff < 5 nm, and over a fairly wide range of field strength and temperature (see Figure 6) after the following modification.Since in the Onsager model the quantum yield depends on the parameter eFr 0 /2kT, 21 the dependence of the exponent n in eq 39 on this parameter is assumed.A good fit to the Onsager model is achieved by increasing the value of n from 1 to 4 with increasing field F. The distance r 0 serves as a spatial scale.Thus, in eq 39 i k j j j y A similar recent modification of the ETC made it possible to obtain an analytical expression for the temperature and field dependence of mobility in materials with correlated disorder. 35n apparently weak field, e.g., 10 6 V/m, which is insignificant at room temperature, significantly increases the quantum yield at low temperatures; see Figure 6a.One takes this circumstance into account when analyzing the temperature dependence of the quantum yield, obtained experimentally or by MC simulation.The concentration of hopping centers and nondiagonal disorder are also of great importance; see Figure 7.The values obtained with the model developed in this work (solid line) significantly exceed the predictions of the Onsager model (with r 0 = r i = 2.4 nm at low temperatures) but remain considerably smaller than the MC data 17 (squares), even after taking into account the uniform field by eqs 39−41 (dashed line).One should note, however, that the MC simulation 17 used a cubic lattice of hopping centers, while our model assumes a broad spatial distribution of hopping centers, arising due to anisotropy of overlapping wave functions. 5As a result, both the typical jump length and the effective initial separation in the case of MC data 17 should be larger than in our model.Hence, one has to replace r 0 ef by r 0 ef + Δr 0 in eq 40.An estimate of the correction to the effective initial separation, Δr 0 , is given in part S3 of the SI.The corrected model results (see the thin lines with small symbols in Figure 7) are in good agreement with the MC data.
The calculations based on the developed model, including the external electric field (see eq 40), also give an excellent agreement with the experimental photogeneration quantum yield, controlled by geminate recombination, 23  In conclusion, a theoretical model describing the time dependence of the relation between the drift shift and diffusion spreading of the carriers for a realistic Gaussian energy distribution of hopping centers showed that the Einstein relation breaks down over a shorter time interval than that assumed in the previous work. 18However, the calculation   results qualitatively confirm the conclusion obtained in this work by analytic modeling (and earlier Monte Carlo simulations 17 ) that the separation probability of a geminate pair at low temperatures is significantly greater than that predicted by the Onsager model, since an anomalously strong diffusion increases the effective initial separation of the pair.A significant contribution to both the increase in the absolute value of the separation probability and the slowdown of its temperature dependence is also made by the external electric field, which cannot be considered weak at low temperatures and early times with a highly nonequilibrium transport.The results of the calculations using simple analytic formulas (see eqs 37 and 40) and their comparison with the MC simulation data 17 also stress the effect of off-diagonal disorder, which is opposite to the effect of diagonal (energy) disorder.For a given concentration of hopping centers, the effective initial separation distance and the probability of subsequent separation are smaller, if off-diagonal disorder is stronger.

■ ASSOCIATED CONTENT
* sı Supporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpclett.4c00662.Appendices with derivations for the expressions for characteristic time dependences, energies and times (S1), solution for the survival probability in the WKB approximation (S2), and correction of the initial separation to the lattice model (S3) (PDF) ■

Figure 1 .
Figure 1.Left: workflow of the derivation.Right: scheme of electron hops in an electric field of strength F. The shaded rectangle shows the energy regions of the states from which hops to a state with a given energy E are possible at a small temperature.

Figure 4 .
Figure 4. Time dependences of demarcation energy, E d (t), in units of σ.DOS is Gaussian with σ = 0.1 eV.Dashed lines show the long-time asymptote, E d (t) = E tr − kT ln(ν 0 t).

Figure 5 .
Figure 5. Temperature dependences of separation probability obtained from the model of this work and from the Onsager model, a 0 = N −1/3 .
Figure 8, providing a realistic value, r 0 ef = 1.75 nm.The experimental data are taken from Figure 10 (open squares) of ref 23.

Figure 6 .
Figure 6.Field and temperature dependences of separation probability, as calculated from the Onsager model (solid lines) and from the modified effective temperature method; see eqs 39−41 (dashed lines).Temperatures in panel (b) are the same as those in panel (a).

Figure 7 .
Figure 7.Comparison of the results of this model (solid line) and its modifications (by the use of effective temperature) for the case of regular lattice and finite electric field, F = 10 6 V/m (thin lines), with the results of the Onsager model (dashed line) and the MC modeling 17 (squares).

Figure 8 .
Figure 8.Comparison of the results according to the model of this work (multiplied by the estimated dissociation probability of the primary excitation, 0.18) with the experimental data 23 for phenylsubstituted polyphenylvinelene (PhPPV) with aluminum electrodes.The field F is the same as in the experiment; see caption of Figure 10 in ref 23.loc = 10, σ = 0.1 eV, r i = 1.4 nm.