Fractional Derivative Phenomenology of Percolative Phonon-Assisted Hopping in Two-Dimensional Disordered Systems

Anomalous advection-diffusion in two-dimensional semiconductor systems with coexisting energetic and structural disorder is described in the framework of a generalized model of multiple trapping on a comb-like structure. The basic equations of the model contain fractional-order derivatives. To validate the model, we compare analytical solutions with results of a Monte Carlo simulation of phonon-assisted tunneling in two-dimensional patterns of a porous nanoparticle agglomerate and a phase-separated bulk heterojunction. To elucidate the role of directed percolation, we calculate transient current curves of the time-of-flight experiment and the evolution of the mean squared displacement averaged over medium realizations. The variations of the anomalous advection-diffusion parameters as functions of electric field intensity, levels of energetic, and structural disorder are presented.


Introduction
The problem of transport in disordered semiconducting systems is important from fundamental as well as from technological point of view.Anomalous (non-Fickian) diffusive behavior of charge carriers is often observed in these materials [1][2][3][4] and can be attributed to various mechanisms, e.g., multiple trapping of charge carriers into localized states with random energies, phonon-assisted tunneling or hopping, percolation over the conduction states, and others.In most cases, disordered semiconducting systems are characterized by the coexistence of energetic and structural disorder.There remains the need for joint account of material morphology and energy disorder to describe charge carrier dynamics, particularly in models of dye-sensitized solar cells (DSSC, the Grätzel cell) [5][6][7][8] or polymer bulk heterojunction (BHJ) [9][10][11][12].
DSSCs traditionally include a highly porous metal oxide nanocrystalline semiconductors, such as TiO 2 , on a transparent conductive electrode.The efficient electronic transport is responsible for the overall performance of DSSC and disordered structure of porous anode affects essentially on transport characteristics [7,8,13,14].Electron transport in mesoporous TiO 2 exhibits main features of dispersive kinetics in disordered semiconductors [15].The models of phonon-assisted hopping (PAH) and multiple trapping (MT) are often used for the description of charge transport in DSSC.In the MT model (see details in [16]), carriers are divided into mobile or free ( f ) and trapped (t) ones.First group, being delocalized, is responsible for the charge transport.In the PAH model, carriers tunnel between localized states with absorption or generation of phonons.Organic photovoltaic (PV) devices based on polymer blends (such as P3HT:PCBM) forming BHJ [9][10][11][12] are also characterized by strong structural and energetic disorders.The idea of BHJ is based on the extension of the donor-acceptor interface on the entire volume of the working (active) layer.The distributed border of large area between pand n-type domains significantly increases dissociation efficiency of photogenerated excitons.
Stochastic modeling plays a key role in studying charge carrier transport in these systems [17].Random walk based simulations are often performed [1,[18][19][20] to elucidate different transport features.The seminal works on simulation of dispersive transport by Bässler [18] for organic semiconductors and Nelson [15] for nanocrystalline TiO 2 electrodes were followed by a great number of applications of this approach to the modeling of solar cells (see review [21]).
To study combined the effect of energy and structure disorder, simulations of random walks in mesoporous TiO 2 were employed in [5,6,22].Benkstein et al. [5] demonstrated that network topology has a strong influence on the electron transport in mesoporous TiO 2 samples.Their measurements and simulations are described well by percolation theory.They studied DC conduction in mesoporous films and obtained the power-law dependence of the electron diffusion coefficient D on the film porosity P, D ∝ |P − P c | µ , where P c is the critical porosity (percolation threshold).They also noticed that the fraction of terminating particles (dead ends) in the mesoporous TiO 2 film increases essentially when P increases.Cass et al. [6] studied the transient response of photoelectrons by means of a multitime-scale Monte Carlo simulation, and also found that the morphology of the TiO 2 samples has a considerable impact on the transport characteristics.
There is a lack of analytical approaches considering the combined effect of energy disorder and morphology.Apparently, such an approach can be developed as a generalization of dispersive transport models (MT, PAH, CTRW), accounting for the geometry of the structure.As is known, there exist several geometrical models to describe percolation in disordered media: the single-conductor network of Skal and Shklovsky [23], the fractal model of Mandelbrot and Given [24], the Arcangelis-Redner-Coniglio model [25], Stanley's droplet model [26], and others.The popular model to study anomalous diffusive transport in low dimensional percolation clusters is the so-called comb model introduced in [27,28].The dynamic equation to describe transport on the continuous comb has been proposed by Arkhincheev and Baskin [29] and studied in [30][31][32][33][34][35].The advection-diffusion in the comb can be described in the framework of the continuous time random walk (CTRW) concept [36,37].Diffusion in a teeth could be considered as a trapping event characterizing by a random waiting time; effective transport is provided by motion along the backbone [35].The CTRW model for TiO 2 films was proposed in [15].Generalizations of the comb model consider trap-limited diffusion [35], multiple trapping [33] or fractional Brownian motion [38] on combs, fractal comb-like geometry [32,34], multidimensional generalization [39], and others.Assuming the coexistence of two subdiffusive processes: diffusion on a fractal structure and heavy-tailed CTRW, Meroz et al. [40] have noticed the subordination of an ergodic anomalous process to a nonergodic one.The resulting process is nonergodic-time averaged and ensemble averaged mean squared displacements (MSDs) do not coincide.
In the present paper, we try to describe the combined effect of energetic and structural disorder for phonon-assisted hopping in two-dimensional disordered systems in the framework of a generalized model of multiple trapping on a comb-like structure.The basic equations of the model contain fractional-order derivatives.To verify analytical predictions, we perform a series of kinetic Monte Carlo (KMC) simulations of hopping in two types of inhomogeneous media: nanoparticle agglomerate and two-component polymer blend.To elucidate the role of directed percolation in joint influence of structural and energetic disorder, the evolution of mean squared displacement and transient current curves in the time-of-flight method are calculated.The variations of the anomalous advection-diffusion parameters as functions of electric field intensity, levels of energetic, and structural disorder are presented.

Multiple Trapping on a Comb-Like Structure
Here, we continue to develop the model of multiple trapping on a comb structure (MT-comb model) proposed in [33] for dispersive transport in disordered nanostructured samples.The classic comb model consists of 'backbone' coinciding with the x axis and regularly distributed perpendicular 'teeth' directed along the y axis [27,28].The effective transport on the comb is possible only along the x axis, and 'teeth' play role of traps in the CTRW interpretation.The comb geometry describes approximately the effect of 'dead ends' of a percolation cluster [29,32,37].To describe hopping on such a structure, we consider the MT model, although applicability of this replacement is not ensured a priori.It will be validated by comparison with Monte Carlo simulation results.
The most important features of PAH have been understood on the basis of the numerical simulations for different densities of states (DoS) and spatial distributions of traps.Considering the equilibrium conduction, authors of Reference [41] came to the conclusion that the main contribution to conductivity is made by electrons, performing jumps near the characteristic energy level.An analytical formula was obtained for the energy showing that it does not depend on the position of the Fermi level.Independently of them, Monroe [42] considered the energy relaxation of electrons by hopping over the band tail states.The author showed that the electrons starting from the mobility edge make a few jumps and quickly relax in energy, but the type of relaxation changes when the characteristic energy transport level ε tr is reached.The hopping process near and below the transport level is similar to the MT relaxation.The transport level can be considered as an analog to the mobility edge.Only states with energies less than ε tr should be interpreted as traps [3,43,44].We extend this analogy for PAH in a comb-like structure.
Consider the non-Markovian generalization of the Fokker-Planck equation (FP-equation) describing dispersive advection-diffusion in disordered semiconductor [45][46][47][48], where Φ(t) is a memory kernel defined by transport mechanism, DoS, distribution of distances between hopping centers, E is an electric field intensity, µ and D are effective mobility and diffusion coefficient, respectively, and G(x) = p(x, 0) is the initial concentration of generated carriers.Here, we consider constant and homogeneous µ and D.
The Laplace transformation of the FP-Equation (1) leads to In the continuous time random walk (CTRW) model [1], the Laplace transform Φ(s) of the kernel is associated with the transform of survival probability P{θ > t} = Ψ(t) and density ψ(t) of waiting times θ in traps by formula: As was mentioned in Introduction, the diffusion in the comb can be formulated in terms of the CTRW model [36,37].Similarly, MT-comb model can be also represented as CTRW.If l is a random number of visited sites of a dead branch during one sojourn of the branch, localized states are in each site and the corresponding waiting times are independent and identically distributed random variables with pdf ψ(t), then the distribution of the sojourn time in dead ends (averaged over l) is Here, p l is a distribution of l.After the Laplace transformation, the convolution transforms into the product of densities For the transform of the integral kernel of the generalized FP-equation, we have In the case of a tempered power-law distribution of the number of visited states in dead ends, we have p(s) and taking expression − ln ψ(s) = − ln 1 − s Ψ(s) ∼ s Ψ(s) into account at s → 0, we arrive at the following expression For asymptotically large times, according to the Tauberian theorem [49], we can use asymptotic transform of the FP-equation for s → 0, In the MT model [16,33] Here, ρ(ε) is DoS, ω 0 is a trapping rate.For exponential DoS ρ(ε) = ε −1 0 exp(−ε/ε 0 ), we have the power law distribution of waiting times in localized states [1,16]: For the non-tempered case (γ = 0) and for exponential DoS, we have As a result, we obtain the advection-diffusion equation containing the fractional derivative of order αβ where is the redefined FP-operator, is the Riemann-Liouville derivative (RL-derivative) (see, e.g., [50]), α = kT/ε 0 is related to energy disorder (DoS width), and β is determined by the level of structural disorder.
If α = 1 and γ = 0, we have the tempered fractional Fokker-Planck equation where is, so called, the tempered fractional RL-derivative [51].In the case γ = 0, we have the standard RL-derivative 0 D β,0 The generalized FP-Equation ( 5) is obtained under the assumption about negligibly small sojourn times in the backbone between subsequent trapping into 'dead ends'.In the more general case, the transport is characterized by the composition of memory kernels, one of which accounts for walking in 'dead ends' and another arises due to sojourn time in the backbone, The set of considered kernels for different situations is represented in Table 1.In case I, the advection-diffusion equation contains the fractional RL-derivative of order αβ.In case II, we meet the tempered fractional RL-derivative of order β with truncation factor γ. Case III corresponds to a particular type of distributed order fractional equation with two-component fractional operator that leads to transient dispersive kinetics characterized by different dispersion parameters αβ and α at different time scales.

Transient Current Curves
The main features of the dispersive transport have been established by means of time-of-flight (ToF) experiments for amorphous semiconductors.In the ToF method, non-equilibrium charge carriers are generated by a light pulse of short duration on the surface layer of a sample situated between two electrodes.These carriers propagate under influence of the applied electric field, that results in a bias current.The applied voltage is high enough to disregard the field of nonequilibrium charge carriers.Non-Gaussian (e.g., dispersive) transport behavior could be identified by observation of the power law decay of photocurrent I(t).Universal TC-curves are characterized by two power law sections, I(t) ∝ t −1+α for t < t T and I(t) ∝ t −1−α for t > t T (α is the dispersion parameter), and power law dependence of transient time on the sample thickness t T ∝ L 1/α [1].In [52], it has been shown that such behavior of transient current indicates the specific self-similarity of charge carrier propagator, which satisfies the advection equation containing the time derivative of fractional order α.
The ToF photocurrent can be found by integration of carrier concentration p(x, t), Function p(x, t) in its turn could be found by solving the corresponding transport equation.The solution to Equation ( 9), under the ToF conditions, can be written in the subordinated form [53] p(x, s) = ∞ 0 dτ f (x, τ) q(τ, s), q(τ, s) = Φ(s) exp(−τs Φ(s)).(12) Here, f (x, τ) is the solution of the classic FP equation (i.e., Equation (1) with first derivative on operational time τ instead of integral operator), and q(τ, s) is the Laplace transform of the operational time density.Substituting (12) into the Laplace transformation of (11), after integration by parts and changing the order of integration, we arrive at equation where )dx is the transient current at normal transport, which can be found by solution of classic FP-equation.As a result, we have If E is a homogeneous electric field, µ = const, D ≈ 0: Even in case, when waiting times in sites are characterized by the exponential distribution, Ψ(s) = 1/(µ + s), the asymptotic (s → 0) term of integral kernel Φ(s) corresponds to a power law memory kernel with exponent β < 1, Φ(s) ∼ c β µ −β s β−1 .So, the integral operator in the generalized FP-Equation (1) becomes the RL fractional derivative in this case and dispersive transport with parameter β < 1 should be observed.
In the case of kernel (7), we have the following Laplace transform of transient current for pure advection in strong electric field Inverse Laplace transformation of this expression leads to Here, g + (τ) is the one-sided Lévy stable density (subordinator), 0 < η ≤ 1; t T is a transient time usually determined through the break on TC curve.

Kinetic Monte Carlo Simulation
To verify the model, we compare analytical predictions with results of a Monte Carlo simulation of phonon-assisted hopping in porous nanoparticle agglomerate and phase-separated bulk heterojunction.

Simulation of Phonon-Assisted Hopping
Considered samples are rectangular (S × L), where L is a sample length along the field direction and S is a transverse size.Hopping is simulated using the Miller-Abrahams relation for the transition rate between localized states with Gaussian or exponential DoS and arranged in square lattice points.The DoS parameters depend on the local composition of the sample.For porous system, the lattice points belonging to voids are excluded.
An electron (or hole) makes a hop with the Miller-Abrahams rates where r jk is a distance between states and a is a localization radius.Amplitude of hopping rates depends on the difference between energies of localized states ∆U jk = ε j − ε k + eEr jk and E is an electric field intensity.During the computation, hopping rates are estimated according to (17).For a carrier in the j-th localized state, γ jk is calculated for all k in the square with center in r j and side l = 5a.For each j → k transition, the random time τ j →k = −γ −1 jk ln U is generated.Here, U is the random variable uniformly distributed in (0, 1).Transition occurs into the k-th state with minimal τ j →k .Under conditions of the ToF-method, carriers are initially generated near the left electrode.The carriers are reflected from the left electrode, and absorbed by the opposite electrode.

Generation of Particle Agglomerate
We consider a simplified geometry of a mesoporous semiconductor to separate effects of percolation and distribution in energy states on transient photocurrents.A two-dimensional pattern of nanoparticles is generated as a set of N n circles with centers randomly distributed over a square sample.Each circle has a random radius from Gaussian distribution with the mean value R and standard deviation δ R .Circles can overlap each other.Porosity of the obtained agglomerate is calculated as a fraction of empty area.Localized states with a given DoS are distributed inside the particles.Hopping rates between localized states depend on spatial separation and the energies of these states according to (17).We consider samples of size 5 × 5 µm 2 and particles with mean radius R = 120 nm and relative fluctuations of their size δR/ R = 0.5.
Note, that electron transport in the mesoporous TiO 2 films in DSSC occurs by diffusion since the macroscopic electric field across the sample is screened by an electrolyte [54].Here, we consider the mesoporous system without electrolyte and calculate transient currents under conditions of the ToF method.Obtained here, results could not be directly applied to DSSC but provide important understanding of the joint influence of structural and energetic disorder.
Electron trajectories of phonon assisted hopping in nanoparticle agglomerate with N n = 800 for different electric field intensities are demonstrated in Figure 1.Localization of charge carriers in topological traps is enhanced at increased field intensities.Irregularity of spatial distribution of these traps is evident for higher electric fields.

Generation of Spinodal Pattern of Two-Component Blend
Another considered system is BHJ.In theory, BHJ is often implemented as a result of spinodal decomposition of a donor-acceptor semiconductor blend.Simulation of spinodal decomposition and related phenomena is based on the solution of the Cahn-Hilliard equation [55,56] Here D is a bounded domain in R n with sufficiently smooth boundaries, n = 2 in our case, function − f is a derivative of double-well potential W, which is often chosen in the form W(u) = (u 2 − 1) 2 /4, and ε is a small positive parameter, characterizing length of interaction.Function u is a concentration of one of the blend components.Extreme values −1 and 1 of function u correspond to fractions 0 and 1 of chosen component.The solution to the Cahn-Hilliard equation satisfies the mass conservation condition Ω udV = µ = const.
When we consider hopping on a two-dimensional (n = 2) spinodal, patterns of abstract donor-acceptor polymer blend.Donor area is considered to be a set of domains in which u < 0.5, other domains are considered to belong to the acceptor area.A charge carrier performs hopping over one of these areas depending on carrier's sign.We consider patterns corresponding to different µ-values (ratio of components) and different annealing times, determining the morphology of the area accessible for hopping.

Calculation of Transient Current, Mobility and Diffusion Coefficient
Transient current is calculated as a density of the generalized renewal process.The time axis in log-scale is divided into a number of intervals.Each hop i → j, occurring at time t ij , belonging to the k-th time interval contributes to histogram of displacement current I k .Mobility is defined by µ = v x /E, where drift velocity can be determined as v x = L/t T .Here L is a sample thickness and t T is a transient time, or 'time of flight'.The latter can be found through the break on the TC-curves, or from the passage of coordinate L by the mean position of the packet.The coefficient of parallel diffusion can be calculated as In this case, the sample is considered as unbounded that is realized by application of periodic boundary conditions.

Transient Currents in 2D Porous Films
The calculated transient current curves for nanoparticle agglomerates were averaged over 100 realizations.We have chosen exponential DoS and energetic dispersion parameter α = kT/ε 0 = 0.75.For each realization we simulated the hopping of 1000 charge carriers.In all cases, we observed universal transient current curves I(t) with slopes in log-log scale close to −1 + η and −1 − η, before and after the break, respectively.The slopes was determined by means of the least square method in the log-log scale for 10 points before and after the break.Approximation error is of order 0.02.For the homogeneous sample, we observed universal curves with slopes −0.25 and −1.75, that correspond to the chosen α.The dispersion parameter of the structural disorder was determined from the simple relation αβ = η according to (16).Results of calculations for different porosity and two values of electric field intensity are presented in Figures 2 and 3, summarized in Table 2.
At increased electric field, we observe decreasing of the transient time t T , as expected, but the dispersion parameter η also decreases that corresponds to amplification of the structural disorder (β decreases).Trajectories and heat charts demonstrate that additional containment of carriers can be induced by topological traps.This containment depends on electric field and porosity and can be explained by peculiarities of directed percolation.Mean square displacements for hopping diffusion (without electric field) in nanoparticle agglomerate with different porosity are presented in Figure 4. We can see that structural dispersion parameter β d (at diffusion) determined from relation ∆ ∝ t αβ d /2 is sufficiently closer to 1 than β for directed percolation under applied electric bias.

Conclusions
We have described the combined effect of energetic and structural disorder for PAH in two-dimensional disordered systems in the framework of the MT-mechanism on a comb-like structure.Analytical predictions have been confirmed by a series of Monte Carlo simulations of hopping in two types of inhomogeneous media: nanoparticle agglomerate and two-component polymer blend.
The equations derived in Section 2 can be used in the generalized one-dimensional bilayer model of Barker et al. [57] accounting for charge photogeneration, injection, drift, diffusion, recombination, and the effect of space charge on the electric field within the device.Their model utilizes the ordinary FP equations and well describes the steady-state characteristics of BHJ solar cells.On the other hand, dispersive kinetics usually observed in OPV materials cannot be adequately described by the standard FP equation [3,16].Fractional differential models have shown themselves effective in the theory of dispersive transport [45,53].The mentioned generalization requires the development of numerical tools to solve the fractional FP-equations.Some steps in this direction have already been taken.Choo et al. [58] proposed the numerical scheme for the system of fractional FP-equation with RL-derivative and Poisson equation to account for the fluctuation of localized electric field due to the evolution of charge carriers.Rekhviashvili and Alikhanov [59] numerically solve the fractional FP-equation with an alternating electric field.These schemes could be generalized for the presented dispersive transport equations for systems with coexisting energetic and structural disorder.
Monte Carlo simulation of charge carrier hopping in 2D mesoporous samples and spinodal patterns have shown that topological traps play a crucial role in dispersion of charge propagation for all values of electric field.These traps act independently of thermalization process.Due to this structural disorder, dispersive transport can be observed even at small values of energy disorder, when the DoS width is smaller than kT.Theoretical estimations based on MT-comb model confirm these results.Topological traps can be responsible for the separation of particles into two groups of 'fast' and 'slow' carriers.This separation is deepened at higher fields.
Transient current curves are sensitive to the realization of disordered samples.Some of them demonstrate several breaks in a log-log-scale, the latter stage of TC-decay can slow down.The latter behavior in measurements in some works (see, e.g., [60]) is interpreted as bipolar transport, but in our case it arises due to trapping into topological traps.Ensemble averaged TC-curves look like typical universal TC curves of dispersive transport.For the case of exponential DoS, the observed dispersion parameter is equal to product of energy dispersion parameter α = kT/ε 0 and structural dispersion parameter β.For high levels of porosity, parameter β is very sensitive to electric field intensity.

Figure 2 .
Figure 2. Transient currents for different porosity of nanoparticle agglomerate.

Figure 4 .
Figure 4. Mean square displacements for hopping diffusion in nanoparticle agglomerate with different porosity.

Figure 5
Figure 5 presents trajectories of hopping of 15 particles on the same spinodal pattern (µ = 0.1) for three values of the electric field.With increasing E, trajectories are drawn.At the same time, there are branches of carrier localization in topological traps not directly related to energy disorder.Transient current curves of hopping in spinodal patterns for two values of mass ratio µ = 0.1 and µ = 0.2 and Gaussian DoS ρ(ε) ∝ exp(−ε 2 /σ) (σ/kT = 3) are presented in Figure 6.Two samples for each value of µ are obtained for different durations of spinodal decomposition (later realization is denoted as 'after annealing').Correlation functions of two patterns are also presented.The first minimum of these functions is a correlation length λ that determines the typical scale of p-n-domain separation.Annealing of samples increases the separation of pand n-domains.Results are summarized in Table3.
Figure 5 presents trajectories of hopping of 15 particles on the same spinodal pattern (µ = 0.1) for three values of the electric field.With increasing E, trajectories are drawn.At the same time, there are branches of carrier localization in topological traps not directly related to energy disorder.Transient current curves of hopping in spinodal patterns for two values of mass ratio µ = 0.1 and µ = 0.2 and Gaussian DoS ρ(ε) ∝ exp(−ε 2 /σ) (σ/kT = 3) are presented in Figure 6.Two samples for each value of µ are obtained for different durations of spinodal decomposition (later realization is denoted as 'after annealing').Correlation functions of two patterns are also presented.The first minimum of these functions is a correlation length λ that determines the typical scale of p-n-domain separation.Annealing of samples increases the separation of pand n-domains.Results are summarized in Table3.