Fundamental characteristic length scale for the field dependence of hopping charge transport in disordered organic semiconductors

TL;DR

This study reveals that the field dependence of hopping charge mobility in disordered organic semiconductors is governed by the localization length, not site concentration, challenging previous assumptions and providing new insights into charge transport mechanisms.

Contribution

It demonstrates that the localization length determines the electric field dependence of mobility, contrasting with the traditional focus on site concentration in the literature.

Findings

Mobility dependence on electric field is governed by localization length a.

Site concentration N^{-1/3} does not influence the field dependence of mobility.

Effective temperature combines electric field and temperature through eFa.

Abstract

Using analytical arguments and computer simulations we show that the dependence of the hopping carrier mobility on the electric field $\mu(F)/\mu(0)$ in a system of random sites is determined by the localization length $a$, and not by the concentration of sites $N$. This result is in drastic contrast to what is usually assumed in the literature for a theoretical description of experimental data and for device modeling, where $N^{-1/3}$ is considered as the decisive length scale for $\mu(F)$. We show that although the limiting value $\mu(F \rightarrow 0)$ is determined by the ratio $N^{-1/3}/a$, the dependence $\mu(F)/\mu(0)$ is sensitive to the magnitude of $a$, and not to $N^{-1/3}$. Furthermore, our numerical and analytical results prove that the effective temperature responsible for the combined effect of the electric field $F$ and the real temperature $T$ on the hopping transport via spatially random sites can contain the electric field only in the combination $eFa$.

Equations

$$g(\epsilon )=\frac{N}{\sigma \sqrt{2\pi }}\exp \left(-\frac{{\epsilon }^2}{2{\sigma }^2}\right).$$

$${\nu }_{ij}={\nu }_0\exp \left(-\frac{2|{\mathbf{r}}_{ij}|}{a}\right)\gamma ({\epsilon }_j-{\epsilon }_i+e\mathbf{F}\cdot {\mathbf{r}}_{ij})$$

$$\gamma (\Delta \epsilon )=\{\begin{array}{ll}\exp (-\Delta \epsilon /kT),& \text{if Δε>0,}\\ 1,& \text{otherwise},\end{array}$$

$$\begin{array}{rl}\mu (F)=& {\mu }_0\exp \left\{-{\left(\frac{2}{3}\frac{\sigma }{kT}\right)}^2\right\}\\ & \times \exp \left\{\tilde{C}\left[{\left(\frac{\sigma }{kT}\right)}^2-B\right]\sqrt{F}\right\},\end{array}$$

$$\mu (T,n,F)\approx \mu (T,n)\varphi (T,F)$$

$$\displaystyle\begin{split}\phi(T,F)=\exp\Bigg{\{}&0.44\left[\left(\frac{\sigma}{kT}\right)^{3/2}-2.2\right]\\ &\times\left[\sqrt{1+0.8\left(\frac{Feb}{\sigma}\right)^{2}}-1\right]\Bigg{\}}\,,\end{split}$$

$$T_{\text{eff}}={\left[T^2+{\left(\gamma \frac{eFa}{k}\right)}^2\right]}^{1/2}$$

$$\mu \propto \exp \left[-C{\left(\frac{\sigma }{kT}\right)}^2\right],$$

$$\begin{gathered} \mathcal{R}({\epsilon }_1,{\epsilon }_2)\mathrm{d}{\epsilon }_1\mathrm{d}{\epsilon }_2= \\ \sum _{\begin{array}{c}i\\ {\epsilon }_i\in [{\epsilon }_1,{\epsilon }_1+\mathrm{d}{\epsilon }_1]\end{array}}\sum _{\begin{array}{c}j\\ {\epsilon }_j\in [{\epsilon }_2,{\epsilon }_2+\mathrm{d}{\epsilon }_2]\end{array}}f({\epsilon }_i)[1-f({\epsilon }_j)]{\nu }_{ij}. \end{gathered}$$

$$\sum _i\sum _j\Rightarrow V\rho ({\epsilon }_1)\mathrm{d}{\epsilon }_1\rho ({\epsilon }_2)\mathrm{d}{\epsilon }_2\underset{0}{\overset{\infty }{\int }}{r}^2\mathrm{d}r\underset{0}{\overset{\pi }{\int }}2\pi \sin \theta \mathrm{d}\theta ,$$

$$\mathcal{R}({\epsilon }_1,{\epsilon }_2)=V\rho ({\epsilon }_1)\rho ({\epsilon }_2)f({\epsilon }_1)[1-f({\epsilon }_2)]\mathcal{F}({\epsilon }_2-{\epsilon }_1),$$

$$\begin{gathered} \mathcal{F}(\Delta \epsilon )=2\pi {\nu }_0 \\ \times {\int }_0^{\infty }{e}^{-2r/a}\left[{\int }_0^{\pi }\gamma (\Delta \epsilon +eFr\cos \theta )\sin \theta \mathrm{d}\theta \right]r^2\mathrm{d}r. \end{gathered}$$

$$\underset{-\infty }{\overset{+\infty }{\int }}\mathcal{R}({\epsilon }_1,{\epsilon }_2)\mathrm{d}{\epsilon }_2=\underset{-\infty }{\overset{+\infty }{\int }}\mathcal{R}({\epsilon }_2,{\epsilon }_1)\mathrm{d}{\epsilon }_2.$$

$$\begin{gathered} f({\epsilon }_1)\underset{-\infty }{\overset{+\infty }{\int }}\rho ({\epsilon }_2)[1-f({\epsilon }_2)]\mathcal{F}({\epsilon }_2-{\epsilon }_1)\mathrm{d}{\epsilon }_2= \\ [1-f({\epsilon }_1)]\underset{-\infty }{\overset{+\infty }{\int }}\rho ({\epsilon }_2)f({\epsilon }_2)\mathcal{F}({\epsilon }_1-{\epsilon }_2)\mathrm{d}{\epsilon }_2 \end{gathered}$$

$$\frac{\mathcal{F}(\Delta \epsilon )}{\mathcal{F}(-\Delta \epsilon )}=\exp \left(-\frac{\Delta \epsilon }{kT}\right)\text{(at F=0)},$$

$$\frac{\mathcal{F}(\Delta \epsilon )}{\mathcal{F}(-\Delta \epsilon )}\approx \exp \left(-\frac{\Delta \epsilon }{k{T}_{\text{eff}}}\right),$$

$$f(\epsilon )\approx {\left[\exp \left(\frac{\epsilon -{\epsilon }_f}{k{T}_{\text{eff}}}\right)+1\right]}^{-1},$$

Full text

I Introduction

Organic semiconductors attract currently much attention in the scientific community as materials desired for applications in modern electronics. The term “organic semiconductors” covers a large class of materials with a broad variety of properties. Organic semiconductors can be fabricated in crystalline form, as for instance, pentacene and ruberene Ostroverkhova (2016). The energy spectrum in such materials has a classical band structure with charge carriers moving as free particles or polarons in the conduction and valence bands. The main focus in research on organic materials is put, however, on organic disordered semiconductors (ODSs), such as polymers and low-molecular-weight systems Silinsh (1970); Bässler (1993, 1981); Borsenberger et al. (1993); van der Auweraer et al. (1994); Pope and Swenberg (1999); Hadziioannou and van Hutten (2000); Brabec et al. (2003); Bruetting (2005); Baranovski (2006); Schwoerer and Wolf (2007); Sun and Dalton (2008); Tessler et al. (2009); Meller and Grasser (2010); Geoghegan and Hadziioannou (2013); Baranovskii (2014); Kuik et al. (2014); Köhler and Bässler (2015); Nenashev et al. (2015); Laquai et al. (2015). The interest in ODSs is caused by their optoelectronic features and by easy manufacturing, as compared to organic crystals. In contrary to crystalline materials, ODSs possess neither structural regularity, nor spatially extended electronic states. Instead, electronic states in ODSs are spatially localized Silinsh (1970); Bässler (1981, 1993); Borsenberger et al. (1993); van der Auweraer et al. (1994); Pope and Swenberg (1999); Hadziioannou and van Hutten (2000); Brabec et al. (2003); Bruetting (2005); Baranovski (2006); Schwoerer and Wolf (2007); Sun and Dalton (2008); Tessler et al. (2009); Meller and Grasser (2010); Geoghegan and Hadziioannou (2013); Baranovskii (2014); Kuik et al. (2014); Köhler and Bässler (2015); Nenashev et al. (2015); Laquai et al. (2015). This happens because the overlap integrals for the weak Van-der-Waals interactions between neighboring structural units (molecules or molecular complexes) in ODSs are much smaller than the energy scale of disorder, which prevents the formation of extended electronic states Silinsh (1970); Bässler (1981, 1993). Therefore, charge transport in ODSs is due to incoherent tunneling (hopping) of charge carriers between localized states that are randomly distributed in space Silinsh (1970); Bässler (1981, 1993); Borsenberger et al. (1993); van der Auweraer et al. (1994); Pope and Swenberg (1999); Hadziioannou and van Hutten (2000); Brabec et al. (2003); Bruetting (2005); Baranovski (2006); Schwoerer and Wolf (2007); Sun and Dalton (2008); Tessler et al. (2009); Meller and Grasser (2010); Geoghegan and Hadziioannou (2013); Baranovskii (2014); Kuik et al. (2014); Köhler and Bässler (2015); Nenashev et al. (2015); Laquai et al. (2015). Our paper deals with the description of charge transport in this hopping regime and the results are valid for ODSs, not for organic crystals.

The most popular theoretical model to describe charge transport in ODSs is the so-called Gaussian Disorder Model (GDM), according to which localized states have a Gaussian energy distributionSilinsh (1970); Bässler (1981, 1993); Oelerich et al. (2012)

[TABLE]

Here, $\sigma$ is the energy scale of the spectrum, usually estimatedBässler (1993) in ODSs to the order of $\sigma\approx 0.1$ eV and $N$ is the concentration of randomly distributed localized states (sites). A typical estimateBaranovski (2006); Baranovskii (2014) for the latter parameter is between $N\simeq 10^{20}$ cm*-3* and $N\simeq 10^{21}$ cm*-3*.

The hopping rates are usually assumed Bässler (1993) to be described by the Miller-Abrahams expressionMiller and Abrahams (1960). For each pair of sites $(i,j)$, the rate $\nu_{ij}$ is determined by their energy difference $\varepsilon_{j}-\varepsilon_{i}$ and position difference $\mathbf{r}_{ij}\equiv\mathbf{r}_{j}-\mathbf{r}_{i}$:

[TABLE]

with

[TABLE]

where $a$ is the localization length of charge carriers, $\mathbf{F}$ is the electric field, and $\nu_{0}$ is a prefactor determined by the tunneling mechanism. The localization length $a$ in ODSs is estimated Gill (1972); Rubel et al. (2004) at the order of $10^{-8}$ cm, which is much smaller than the intersite distance $N^{-1/3}$. Therefore, we follow the usual assumption Silinsh (1970); Bässler (1981, 1993); Borsenberger et al. (1993); van der Auweraer et al. (1994); Pope and Swenberg (1999); Hadziioannou and van Hutten (2000); Brabec et al. (2003); Bruetting (2005); Baranovski (2006); Schwoerer and Wolf (2007); Sun and Dalton (2008); Tessler et al. (2009); Meller and Grasser (2010); Geoghegan and Hadziioannou (2013); Baranovskii (2014); Kuik et al. (2014); Köhler and Bässler (2015); Nenashev et al. (2015); Laquai et al. (2015) that $a$ can be considered to be independent of the concentration of sites $N$.

While powerful and transparent analytical theoretical tools have been developed to describe the dependencies of the hopping mobility $\mu$ on $T$, $N$, $a$, $\sigma$, and on the concentration of carriers $n$, as highlighted in recent reviews Baranovski (2006); Tessler et al. (2009); Baranovskii (2014); Nenashev et al. (2015); Ostroverkhova (2016), theoretical studies of the dependence $\mu(F)$ have mostly been focused on computer simulations. The group of Bässler simulated $\mu(F)$ on a cubic lattice and fitted results in the form of the parameterized equation Borsenberger et al. (1991); Bässler (1993); Köhler and Bässler (2015)

[TABLE]

where $\mu_{0}$ is a field-independent prefactor.

Two parameters, $\widetilde{C}$ and $B$ are involved in this fitting. The parameter $\widetilde{C}$ is assumed to depend on the lattice constant $b$ (distance between localization sites) having the value $\widetilde{C}=2.9\times 10^{-4}$ cm1/2/V1/2 for $b=0.6$ nm Borsenberger et al. (1991); Bässler (1993); Köhler and Bässler (2015). Although simulations were performed on regular cubic grids, a non-diagonal disorder has been introduced into simulations by Bässler et al. Borsenberger et al. (1991); Bässler (1993); Köhler and Bässler (2015) in order to mimic spatial disorder. The exponent $(2|\mathbf{r}_{ij}|/a)$ in Eq. (2) was rewritten in the form $2\lambda b|\mathbf{r}_{ij}|/b$ , where $b$ is the lattice spacing, and the parameter $\lambda$ can be viewed as the inverse localization length. The factor $2\lambda b$ was distributed in a Gaussian manner with the width $\Sigma$ around the value $2\lambda b=10$. The parameter $B$ in Eq. (4) was set equal to $B=2.25$ for $\Sigma<1.5$ and to $B=\Sigma^{2}$ for $\Sigma\geq 1.5$. Equation (4) is one of the most frequently used equations in the context of organic semiconductors Tessler et al. (2009); Baranovskii (2014).

A similar approach to determine $\mu(F)$ was used by Pasveer et al. Pasveer et al. (2005), who reduced the lattice GDM of Bässler et al. to the case $\Sigma=0$ and herewith completely eliminated spatial disorder. Calculating numerically $\mu(F)$ in the framework of this reduced GDM on a cubic lattice, Pasveer et al. fitted results to the analytical formula

[TABLE]

with $\phi(T,F)$ in the form

[TABLE]

where $b$ is the lattice constant. The latter equations are sometimes considered universal and they are the basis Coehoorn and Bobbert (2012) for the commercially available OLED simulation software tools [Simulation software SETFOS3.2, product of Fluxim (www.fluxim.com); Simulation software SimOLED3.x, product of Sim4tec (www.sim4tec.com)].

Pasveer et al. Pasveer et al. (2005) mentioned that Eq.(6) “should merely be considered as a description of the numerical data in a limited parameter range” promising to rationalize this parametrization in future work. We show below that neither Eq. (4) nor Eq. (6) can be rationalized because they do not contain decisive parameters responsible for the field-dependent mobility $\mu(F)$. Equations (4) and (6), which are used by thousands of researchers, were obtained by fitting the numerically simulated data under the assumption that the decisive parameter for the dependence $\mu(F)$ is the intersite distance, parameter $b$ in Eq. (6). We rigorously prove below that this assumption is wrong and the intersite distance is irrelevant for the field-dependent mobility in disordered systems. One should instead use the localization length $a$ as the decisive length scale determining the field dependence of the hopping carrier mobility $\mu(F)$. A theoretical recipe on how to describe the dependence $\mu(F)$ in disordered materials will be formulated below, which should encourage researchers to reanalyze their data on $\mu(F)$ in disordered organic semiconductors.

The paper is organized as follows. In Sec. II, we first stay for simplicity in the framework of the reduced GDM used by Pasveer et al. Pasveer et al. (2005), i.e., on a cubic lattice without spatial disorder. We show that already in this oversimplified case, Eq. (4) and Eq. (6) are incompatible with each other even if the same material parameters in these equations are used. We further show that the results of computer simulations by Pasveer et al. Pasveer et al. (2005) are incompatible with the results of computer simulations by Bässler et al. Borsenberger et al. (1991); Bässler (1993); Köhler and Bässler (2015) carried out in the framework of the same reduced GDM on the cubic lattice (i.e. for $\Sigma=0$). Performing our own computer simulations, we prove that the localization length $a$, not even present in Eqs. (4) and (6), is responsible for this discrepancy in the simulations and that the localization length affects decisively the field dependence of carrier mobility.

In Sec. III, we consider the GDM on spatially random sites, i.e., not anymore on a lattice, and show by computer simulations that the localization length $a$ is the only spatial scale responsible for the field—-dependent hopping mobility. Our computer simulations show herewith that the intersite distance, present in the form of lattice constant $b$ in Eq. (6), is irrelevant for the field—-dependent mobility $\mu(F)$.

In Sec. IV, we further show by computer simulations that the dependence of the carrier mobility on the electric field $F$ can be described by inserting the field-—dependent effective temperature $T_{\text{eff}}(F,T)$, instead of the real temperature $T$, into the temperature dependence of the hopping mobility, which is well understood and described at low electric fields Baranovski (2006); Baranovskii (2014); Nenashev et al. (2015). Herewith our computer simulations on random sites rigorously prove the idea by Shklovskii et al. Shklovskii (1973); Shklovskii et al. (1990); Marianer and Shklovskii (1992); Baranovskii et al. (1993); Cleve et al. (1995); Jansson et al. (2008), who already suggested many years ago that the field—-dependent effective temperature, which contains the localization length $a$ as the only relevant spatial parameter, describes the combined effects of electric field and temperature on the hopping mobility.

In Sec. V, we prove the concept of the effective temperature for spatially random sites by analytical calculations. It is shown that the effective temperature does exist and that it depends on the localization length $a$, and not on the concentration of sites $N$.

Concluding remarks are gathered in Sec. VI.

A short version of this work has been made publicly available in Ref. Nenashev et al. (2017).

II Localization length affects $\mu(F)$ in the lattice model

Before considering a realistic case of a spatially disordered system, let us analyze the simulated data on the field–dependent mobility $\mu(F)$ available in the literature Bässler (1993); Pasveer et al. (2005), which were obtained on regular cubic lattices and served for parametrizations by Eqs. (4) and (6). The concentration of sites $N={b^{-3}}$ is used for the plots in Fig. 1 in order to consist with the data in other figures calculated for random sites.

Let us first check the compatibility of Eqs. (4) and (6) with each other. In order to enable the comparison, we plot the data of Bässler et al. Borsenberger et al. (1991); Bässler (1993) for the case $\Sigma=0$, i.e., with $B=2.25$, since Pasveer et al.Pasveer et al. (2005) simulated for $\Sigma=0$. The value $T=$300\text{,}\mathrm{K} was used in simulations by Bässler *et al*., which gives $\sigma=$0.075\text{\,}\mathrm{eV} for $\sigma/kT=3$. Using the realistic value $\sigma/kT=3$, we plot by a dotted line in Fig. 1 the curve for $\mu(F)/\mu(0)$ given by Eq. (4) and by a dashed line the curve for the function $\phi(T,F)$ given by Eq. (6). The difference in the dependencies $\mu(F)$ given by Eqs. (4) and (6) for the same sets of parameters is striking.

Trying to reveal the reason for such a large discrepancy, we also plot in Fig. 1 the simulated data Bässler (1993); Pasveer et al. (2005) that served as the basis for fittings by Eqs. (4) and (6). The apparent inability of Eq. (4) to fit the simulated data evidences the poor accuracy of this equation but it can hardly be considered as an issue of fundamental importance. However, it is surely an issue of fundamental importance to elucidate the difference in the results of the two simulations Bässler (1993); Pasveer et al. (2005) for the same value $\sigma/kT=3$ because the difference between the data obtained in simulations by the group of BässlerBorsenberger et al. (1991); Bässler (1993) and by Pasveer et al.Pasveer et al. (2005) is comparable to the total effect of $F$ on $\mu$.

The apparent difference in the simulated systems lies in the choice of the parameter $b/a$. While the group of Bässler simulated for $b/a=5$, Pasveer et al. simulated for $b/a=10$. In order to check the validity of those previous simulations, we carried out simulations on a cubic lattice similar to those carried out by Bässler et al. and by Pasveer et al.. Our data, plotted in Fig. 1 for $b/a=10;5;3$, confirm the data by Pasveer et al. with $b/a=10$ and the data by Bässler et al. with $b/a=5$, implying that the computer simulations by both research groups Borsenberger et al. (1991); Bässler (1993); Pasveer et al. (2005) were correct. However, it has not been recognized in previous simulations that the shape of the dependence $\mu(F)/\mu(0)$ is sensitive to the choice of $b/a$.

This result shows that neither Eq. (4) nor Eq. (6) can be considered as universal because these equations do not even contain the localization length $a$. Furthermore this result shows the apparent deficiency of doing physics by computer simulations. Parameterized phenomenological equations, such as Eq. (4) and Eq. (6), do not contain the material parameter $a$, which is decisive for the field–dependent mobility $\mu(F)$, as evidenced in Fig. 1.

Being interested in the dependence $\mu(F)$ for realistic spatially disordered systems rather than for cubic grids, we will consider in the rest of this paper a system of sites distributed in space randomly.

III Localization length determines $\mu(F)$ for random sites.

In order to discern the decisive length scale ($a$, $N^{-1/3}$, or some combination of these parameters) for the field dependence of $\mu$ in a system of random sites, we performed computer simulations using the standard Monte Carlo procedure. A disordered system is created with $140\times 140\times 140$ sites distributed randomly in a box of $L=140$, so that the average inter-site distance $N^{-1/3}$ is unity. The site energies are chosen randomly according to the DOS given in Eq. (1). A single electron is placed onto a random site $i$ and in each simulation step performs a hopping transition to another site $j$ with probabilities weighted by the MA hopping rates given by Eqs. (2) and (3). After each hop, the system time is advanced by $\tau=\nu_{ij}^{-1}$. Initially, the electron is allowed to make $5\times 10^{7}$ relaxation hops to ensure steady-state conditions, after which statistics is collected for $5\times 10^{8}$ hopping transitions. The simulations were repeated and averaged $20$ times. The realistically chosen parameters were $\sigma/kT=4$ and $0.18\leq a/N^{-1/3}\leq 0.30$.

Our simulation results for $\mu(F)/\mu(0)$ are plotted versus $FN^{-1/3}$ in Fig. 2 and versus $Fa$ in Fig. 3. The results look really remarkable. While the plots as a function of $FN^{-1/3}$ differ from each other for different values of the parameter $a/N^{-1/3}$, as they do in the case of the lattice model shown in Fig. 1, the data fall onto a universal curve when plotted as a function of $Fa$. The deviations for $a=0.18N^{-1/3}$ at high $F$ are caused by the effect of the negative differential conductivity discussed elsewhereLien and Shklovskii (1981); Nenashev et al. (2008). The universality of plots $\mu(F)/\mu(0)$ versus $F$ in units $\sigma/ea$ proves that the localization length $a$, and not the intersite distance $N^{-1/3}$ (present in Eq. (6) in the form of the lattice constant $b$), is the decisive length scale for the dependence $\mu(F)/\mu(0)$.

IV Effective temperature dependent on $a$ is responsible for $\mu(F)$

One might be tempted to invent new phenomenological fitting equations in the spirit of Eqs. (4) and (6) for $\mu(F)/\mu(0)$ that would take into account the effect of the parameter $a$. Instead we suggest to recall the idea by Shklovskii, who recognized the importance of the localization length $a$ for the dependence $\mu(F)$ more than 40 years agoShklovskii (1973).

Let us try to understand, why the localization length $a$ and not the intersite distance $N^{-1/3}$ is the decisive length scale for the field dependence of $\mu(F)/\mu(0)$. ShklovskiiShklovskii (1973) considered the case $T=0$ and pointed out that when a charge carrier tunnels in the field direction over some distance $x$, its energy gain due to the applied electric field amounts to $\delta=eFx$. The tunneling probability $\nu(x)\propto\exp(-2x/a)$ can then be rewritten as $\nu(\delta)\propto\exp(-\delta/kT_{\text{eff}})$ with $T_{\text{eff}}\simeq eFa/2$.

For the case of finite temperatures, i.e., for $T\neq 0$, ShklovskiiShklovskii (1973); Shklovskii et al. (1990); Marianer and Shklovskii (1992) and successorsBaranovskii et al. (1993); Cleve et al. (1995); Jansson et al. (2008) have shown that the combined effects of the electric field $F$ and temperature $T$ on the hopping mobility can be expressed in the form of the so-called effective temperature

[TABLE]

with $\gamma\approx 0.67$.Marianer and Shklovskii (1992); Baranovskii et al. (1993)

This result is non-trivial and it looks counterintuitive. The electric field enters the theory only via Eq. (2), i.e., via the combination $e\mathbf{F}\cdot\mathbf{r}_{ij}$. The length of a hop $|\mathbf{r}_{ij}|$ is of the order of the intersite distance $N^{-1/3}$. Therefore, one might expect the combination of parameters $eN^{-1/3}F$ to be essential for the field-dependent mobility. Shklovskii instead argued Shklovskii (1973) that the localization length $a$, i.e. the feature of a single localized state, and not the intersite distance $N^{-1/3}$ is responsible for $\mu(F)$. Taking into account that the Stark effect (determined by the length $a$) is not considered, this proposition sounds revolutionary. Only very recently it has been proven Nenashev et al. (2017) that indeed $a$ and not the intersite distance $N^{-1/3}$ is responsible for the dependence $\mu(F)/\mu(0)$, as described in Sec. III.

The counterintuitive and revolutionary nature of Shklovskii’s idea might be the reason for the fact that it was ignored by the broad scientific community. For instance, in recent review papers Kuik et al. (2014); Laquai et al. (2015), Eqs. (4) and (6) are considered as the main theoretical achievement in the study of charge transport in ODSs. Another possible reason might be the lack of a straightforward proof for this rather counterintuitive concept. Notably, it has never been shown before that only $a$, and not the intersite distance $N^{-1/3}$, is responsible for the dependence $\mu(F)/\mu(0)$, as the parameter $N^{-1/3}/a$ was always fixed and not varied in the simulations confirming the validity of Eq. (7) Marianer and Shklovskii (1992); Baranovskii et al. (1993); Cleve et al. (1995); Jansson et al. (2008). For instance, Marianer and Shklovskii Marianer and Shklovskii (1992) suggested Eq. (7) as the result of computer simulations using the fixed value $N^{-1/3}/a=3$. Their result can be plotted as a function of $eaF$, and, with the same success, as a function of $eN^{-1/3}F/3$. The data in Sec. III and in the previous paper Nenashev et al. (2017) prove, however, that the localization length $a$, and not the intersite distance $N^{-1/3}$ is responsible for the dependence $\mu(F)/\mu(0)$, as assumed in the concept of the effective temperature.

Following this concept Shklovskii (1973); Shklovskii et al. (1990); Marianer and Shklovskii (1992); Baranovskii et al. (1993); Cleve et al. (1995); Jansson et al. (2008), the dependence of the charge carrier mobility $\mu(F)$ on the electric field in hopping conduction can be obtained by inserting the effective temperature $T_{\text{eff}}$ on the place of the laboratory temperature $T$ in the analytical expressions for $\mu(T)$ obtained at low $F$. The temperature dependence of hopping mobility $\mu(T)$ in the GDM at low carrier concentrations is known to have the formBässler (1993)

[TABLE]

where the coefficient $C$ has typically the value $C\approx 0.4$, only slightly depending on the ratio $N^{-1/3}/a$.Baranovskii et al. (2000) In Fig. 4, the mobility $\mu$, obtained in our simulations is plotted as a function of $(\sigma/kT_{\text{eff}})^{2}$, where $T_{\text{eff}}$ is given by Eq. (7) with $\gamma=0.67$. The results perfectly agree with the prediction of Eq. (8) with $T=T_{\text{eff}}$, $C=0.37$, as shown in Fig. 4 by the solid line. Simulations for Fig. 4 were carried out for the parameters sets $\sigma/kT$ between 3 and 4 with the step size 0.25 and $eFN^{-1/3}/\sigma$ between 0.1 and 3.9 with the step size 0.2. The values of $\mu$ in Fig. 4 are normalized by the mobility values at highest $F$ and $T$.

Experimental data for the field-dependent mobility at low carrier concentrations should be compared not with Eq. (4), or Eq. (6), but rather with Eq. (8), in which temperature $T$ is replaced by the field-dependent effective temperature $T_{\text{eff}}$ given by Eq. (7). Such a comparison allows one to determine the value of the localization length $a$ experimentally.

At high carrier concentrations $n$, the temperature dependence of the mobility is described Baranovskii et al. (2002); Baranovski (2006); Baranovskii (2014); Nenashev et al. (2015) by the Arrhenius law instead of Eq. (8). In order to describe the dependence $\mu(F)$ in this regime, one should replace the temperature $T$ in the Arrhenius equation with the effective temperature given by Eq. (7).

V Effective temperature proven analytically

In Sec. IV, we provided analytical arguments for the validity of the effective–temperature concept, where $T_{\text{eff}}$ only depends on the localization length $a$, as suggested by Shklovskii Shklovskii (1973) at $T=0$. Below we provide additional arguments valid also at finite $T$ in favor of the localization length $a$ as the decisive spatial scale responsible for $\mu(F)/\mu(0)$ in a hopping motion of charge carriers via random sites. The results proven in the rest of this report can be formulated as follows.

i) The energy distribution of charge carriers does not depend on the concentration of sites $N$ at a fixed ratio $n/N$.

ii) At non-zero field and temperature the distribution function is the Fermi function with the effective temperature, independent of the site concentration $N$.

The applied electric field changes the rates of carrier transitions between sites, so that site occupation probabilities at nonzero field can differ from their equilibrium values given by the Fermi–Dirac distribution. Here, we will show that one can find the occupation probabilities in the presence of an electric field from a simple integral equation, Eq. (14), assuming that these probabilities are the same for all sites of the same energy, as is granted for zero field. Herewith the occupation probability of some site $i$ is a function $f(\varepsilon_{i})$ solely of the site energy $\varepsilon_{i}$, as in the case at zero field. Note that the energy of a site is considered without the contribution of the potential of the external field. The sites are considered as randomly placed in space, without correlations between their positions and energies. No further assumptions will be involved.

Consider all the carrier transitions from sites within some energy range $[\varepsilon_{1},\,\varepsilon_{1}+\mathrm{d}\varepsilon_{1}]$ to sites within another range $[\varepsilon_{2},\,\varepsilon_{2}+\mathrm{d}\varepsilon_{2}]$, where widths $\mathrm{d}\varepsilon_{1}$ and $\mathrm{d}\varepsilon_{2}$ are small compared to $kT$. Let us denote as $\mathcal{R}(\varepsilon_{1},\varepsilon_{2})\,\mathrm{d}\varepsilon_{1}\,\mathrm{d}\varepsilon_{2}$ the number of such transitions per unit time in the whole sample. Then,

[TABLE]

Since site positions and energies are uncorrelated, the vectors $\mathbf{r}_{ij}$ are uniformly distributed over the three-dimensional vector space with the density $V\rho(\varepsilon_{1})\,\mathrm{d}\varepsilon_{1}\,\rho(\varepsilon_{2})\,\mathrm{d}\varepsilon_{2}$, where $V$ is the volume of the sample. If $V$ is large enough, the vectors $\mathbf{r}_{ij}$ fill the space densely enough to enable integration instead of summation in Eq. (9),

[TABLE]

where polar coordinates $r$ and $\theta$ are introduced in the space of vectors $\mathbf{r}_{ij}$. Directing the polar axis along the field $\mathbf{F}$, and taking into account that $\varepsilon_{i}=\varepsilon_{1}$ and $\varepsilon_{j}=\varepsilon_{2}$ to the accuracy of $\mathrm{d}\varepsilon_{1}$ and $\mathrm{d}\varepsilon_{2}$, one obtains from Eqs. (9)–(10) the following representation for the rate $\mathcal{R}(\varepsilon_{1},\varepsilon_{2})$,

[TABLE]

where

[TABLE]

Now it becomes easy to formulate the carrier balance equation in the steady state. The rate of carrier transitions from the vicinity of energy $\varepsilon_{1}$ to the vicinity of $\varepsilon_{2}$ is proportional to $\mathcal{R}(\varepsilon_{1},\varepsilon_{2})$, and the rate of reverse transitions is proportional to $\mathcal{R}(\varepsilon_{2},\varepsilon_{1})$. Integration of these rates over $\varepsilon_{2}$ provides the total carrier loss from/gain to the energy $\varepsilon_{1}$, and the equality of loss and gain determining the steady state, takes the form

[TABLE]

Inserting Eq. (11), one obtains the following balance equation:

[TABLE]

This is the master equation for calculating the carrier distribution function $f(\varepsilon)$ in the presence of the external electric field. It proves that $f(\varepsilon)$ does not depend on the site concentration $N$. Indeed, since $N$ contributes to this equation only as a factor in the density of states $\rho(\varepsilon_{2})$, it is present in both sides of the equation, and the factors $N$ cancel. Therefore, the electric field $F$ affects the carrier energy distribution only in the combination $eFa$, but not in the combination $eFN^{-1/3}$. This supports our data obtained by Monte Carlo simulations depicted in Fig. 3.

The question might arise on how sensitive this conclusion is with respect to the choice of the expression for transition rates. So far we considered the Miller-Abrahams expression given by Eqs. (2), (3). We would like to emphasize that the electric field $F$ affects the carrier energy distribution only in the combination $eFa$ also for all other shapes of the transition rates, in which the distance $r$ of a hop appears in the combination $r/a$. If the transition rate can be represented in the form of Eq. (2), Eqs. (9)-(14) keep their form. Therefore the conclusion about the decisive role of the localization length $a$ as the only relevant length scale for the dependence $\mu(F)/\mu(0)$ is valid also for Marcus transition rates, which have the form of Eq. (2), though depending on the matrix renormalization energyMarcus (1964).

Furthermore, Eq. (14) provides a basis for justifying the concept of the effective temperature at nonzero temperatures $T$. To see this, let us first note that at zero electric field, the following relation holds for any $\Delta\varepsilon$:

[TABLE]

as evident from Eqs. (3) and (12). In this case, according to detailed balance, the solution of the master equation Eq. (14) is the Fermi–Dirac distribution. Similarly, if (at nonzero electric field) there is such a quantity $T_{\text{eff}}$ that for any $\Delta\varepsilon$

[TABLE]

then the solution $f(\varepsilon)$ of Eq. (14) should have the form of the Fermi–Dirac function with the effective temperature $T_{\text{eff}}$ instead of the real temperature $T$:

[TABLE]

with an appropriate value of the Fermi energy $\varepsilon_{f}$.

Using Eq. (12), we numerically checked the validity of the relation in Eq. (16) for the whole range of electric fields at $\sigma/kT=3$ and $\sigma/kT=4$. This relation is proven to hold in the range of energy differences $\Delta\varepsilon\in[-\sigma^{2}/kT,\;\sigma^{2}/kT]$, which has been provenBässler (1993); Baranovski (2006); Tessler et al. (2009); Baranovskii (2014); Nenashev et al. (2015); Ostroverkhova (2016) as decisive for charge transport in organic semiconductors. Moreover, in the limit of small carrier concentrations we verified that the solution $f(\varepsilon)$ of Eq. (14) follows Eq. (17) in the important energy range $\varepsilon\in[-\sigma^{2}/kT,\;0]$. Herewith it is apparent from the above consideration that the effective temperature introduced in Eqs. (16) and (17) cannot depend on the site concentration $N$, so that the electric field contributes to the effective temperature only in the combination $eFa$.

VI Conclusions

By computer simulations and by analytical calculations we showed that the localization length $a$ of charge carriers in the localized states is the only spatial parameter responsible for the dependence of the hopping mobility $\mu$ on the applied electric field $F$ in a system of random sites. Remarkably, this parameter $a$ is not present in Eqs. (4) and (6), which are often treated as theoretical predictions for $\mu(F)$ and used in device simulations. Results of the current report exclude Eqs. (4) and (6) as candidates for the description of the dependence $\mu(F)/\mu(0)$. Instead, the effective temperature that contains the electric field $F$ in the combination $eFa$ is responsible for the dependencies of the carrier mobility $\mu$ on $T$ and $F$, as illustrated in Fig. 4.

In essence, theories developed for the temperature–dependent hopping mobility $\mu(T)$ at vanishingly small electric fields $F$, as described in recent reviews Baranovskii (2014); Nenashev et al. (2015), are capable to account also for the field–dependent mobility at high $F$ if the temperature $T$ in the low–field theories is replaced by the field–dependent effective temperature $T_{\text{eff}}(F,T)$ given by Eq. (7).

Acknowledgements.

Authors are indebted to Boris Shklovskii for valuable comments. Financial support of the Deutsche Forschungsgemeinschaft (GRK 1782) is gratefully acknowledged.