Interpretation of diffusion coefficients in nanostructured materials from random walk numerical simulation

We make use of the numerical simulation random walk (RWNS) method to compute the ‘‘jump’’ diffusion coefficient of electrons in nanostructured materials via mean-square displacement. First, a summary of analytical results is given that relates the diffusion coefficient obtained from RWNS to those in the multiple-trapping (MT) and hopping models. Simulations are performed in a three-dimensional lattice of trap sites with energies distributed according to an exponential distribution and with a step-function distribution centered at the Fermi level. It is observed that once the stationary state is reached, the ensemble of particles follow Fermi–Dirac statistics with a well-defined Fermi level. In this stationary situation the diffusion coefficient obeys the theoretical predictions so that RWNS effectively reproduces the MT model. Mobilities can be also computed when an electrical bias is applied and they are observed to comply with the Einstein relation when compared with steady-state diffusion coefficients. The evolution of the system towards the stationary situation is also studied. When the diffusion coefficients are monitored along simulation time a transition from anomalous to trap-limited transport is observed. The nature of this transition is discussed in terms of the evolution of electron distribution and the Fermi level. All these results will facilitate the use of RW simulation and related methods to interpret steady-state as well as transient experimental techniques.


Introduction
Electron transport is a key factor in the functioning of the new generation of photovoltaic and optoelectronic devices for lowcost applications.These devices are normally based in mesoporous materials and nanocomposites.Anomalous or dispersive transport 1,2 features are usually observed in these materials, like extremely slow transport combined with density-dependent diffusion coefficients 3,4 and power-law instead of exponential decays. 5Diffusion of electrons injected into mesoporous TiO 2 used in dye-sensitized solar cells (DSC) is one of the most studied examples. 6,7][10] These traps contribute to place the effective diffusion coefficient several orders of magnitude below its bulk value.The description of electron transport in these situations is made most often on the basis of two models that incorporate different physical assumptions.First, the multiple-trapping formalism (MT) 11,12 assumes that charge transport occurs by displacement in extended states slowed down by a succession of trapping and detrapping events in a network of sites with a distribution of trap energies. 13The second model is the hopping approach.
Here the displacement takes place via transitions between localized states.A recent overview of the application of these models in disordered semiconductors has been presented. 14n alternative route to study the electron transport in mesoporous photovoltaic systems is to use Monte Carlo simulation. 5,15,16The random walk numerical simulation (RWNS) is a stochastic computational procedure that allows for a flexible description of electron transport in a network of traps without huge computational demands.This is especially useful in the context of nanostructured materials since the existence of spatial disorder coupled with a broad distribution of trap energies is characteristic of these systems.RWNS simulations have been efficiently used to obtain transient currents and electron lifetimes, 5,17,18 steady-state electron mobilities 17 and surface photovoltage transients. 15,16It has also been applied to study the effect of grain morphology in electron transport in DSCs. 19It is known 17 that RW simulation leads to a Fermi-Dirac distribution when the calculation reaches the stationary state and that photoconductivities follow a power-law when plotted versus the overall electron density.RW simulation has been also utilized by van de Lagemaat and coworkers 18,20 to interpret diffusion coefficients extracted from experiments.
Each of the mentioned approaches has its own advantages as well as limitations.MT (and hopping) is based on well tested physical assumptions about the basic electronic transitions in the system.However, in order to extract experimental conclusions the models require assumptions about spatial homogeneity of the distributions of electronic states that are not well suited to all practical situations.On another hand, the RWNS is rather versatile with respect to specific geometries or initial conditions in a given experiment, but implementing the extended states of the conduction band in RWNS is not feasible and is usually not attempted.The three approaches (MT, hopping, RWNS) require a certain distribution of localized states as an input.Other relevant parameters are the total number of localized states or traps, the total electron density and the attempt-to-jump frequency.All these parameters should be obtained from other sources or can be taken as adjustable in practical implementations. 21herefore, while the different models and approaches share common features, there are also significant differences that complicate the interpretation of the results obtained by the different methods.This paper presents a comprehensive study of the ability of the RWNS method to compute electron diffusivities in nanostructured materials and its relationship to the MT and hopping formalisms.In the theory section (section 2), we provide a detailed discussion of the different frameworks of interpretation, both the physical-analytical and numerical models.By explicit calculation we give formulas for the diffusion coefficient as a function of Fermi level in each of the models and this clarifies the connection between the models in steady state conditions.In section 3 we show the calculations of diffusion coefficients in RWNS as a function of Fermi level, which confirms the previous analytical results.Another interesting point is the use of RWNS in transient experiments.This is specially pointed in view of the application of the numerical simulation method to describe short time and space dynamics of electrons in nanostructured TiO 2 measured by surface photovoltage transients. 16In section 4, we present calculations of the transient behaviour and we show that the results are well described by steady-state diffusion coefficient provided that certain equilibration conditions are satisfied.We finish with the main conclusions.

Definition of diffusion coefficient
The random walks of an electronic carrier determine the jump diffusion coefficient, that has the form 22,23 where Dr i is the displacement of the ith particle at time t, and hi denotes a statistical average.More precisely, the jump (or kinetic) diffusion coefficient defined by eqn (1) reflects diffusion of the center of mass of N particles, while the tracer diffusion coefficient, D*, reflects random walks of a particle If on average, there are no cross correlations between displacements Dr i (t) of different particles at different times, D J and D* become equivalent. 22,23Monte Carlo simulations show that jump and tracer diffusion coefficients are practically identical in a broad interval of carrier densities and tempera-tures. 24Hence, we assume hereafter D* = D J .The jump diffusion coefficient can often be expressed as [24][25][26] in terms of a mean effective jump frequency hni, and the square of effective jump length hr 2 i.
On another hand, experimental information on the fundamental jump rates is often derived from the chemical diffusion coefficient, D n , that relates the flux J n to the gradient of the electron density, n, by Fick's law The diffusion coefficients D n and D J differ by the quantity w n 22,25,26 that is called the thermodynamic factor and is defined as follows where m is the chemical potential of electrons.
The definition of the mobility, u n , is given in terms of the average carrier velocity hv(F)i acquired under electrical field F, at low field values In the presence of an electrical field F, at equilibrium, a local electrical field F is compensated by an opposite variation of the gradient of the chemical potential (thermodynamic driving force): qF = qm/qx = (qm/qn)(qn/qx).Taking into account eqn (4), ( 6) and (7), the generalized Einstein relation is obtained 6,14 The Einstein relation can also be expressed in the form of the classical Einstein relationship However, D J is not in general the diffusion coefficient appearing in Fick's law.

Exponential distribution of localized states
There is wide agreement 27 that random nanoparticulate TiO 2 networks used in DSC exhibit an exponential distribution of localized states (DOS) in the bandgap that is expressed as where N L is the total trap density, k B T 0 the width of the distribution, and E is the energy distance with respect to E 0 , the lower (higher) edge of the conduction (valence) band, for electrons and holes, respectively.E 0 indicates the level of extended states in MT model.In nanostructured TiO 2 layers, values of T 0 range usually between 500 28 and 900 K 29 depending on preparation and measurement conditions.Characteristic values of the total trap density are in the range N t = 10 20 -10 21 cm À3 for mesoporous TiO 2 . 30he exponential distribution has the following property 31 w n = T 0 /T (11)   For the typical values of T 0 , w n E 2-5 at room temperature.Note that eqn (11) is valid only for a deep distribution such that T/T 0 o 1.The diffusion-mobility ratio, eqn (8), is here independent of temperature 32

Multiple trapping
The MT model has been widely used in the DSC area to explain the non-linear dependency of the electron diffusion coefficients and photoconductivities on light intensity and electrical bias (see ref. 6 and 14 for recent reviews).It also permits prediction of the activation energies usually observed in DSC, 12 and the observed correlation between electron lifetime and diffusion coefficient. 11The MT model has been also used to provide a quantitative account of the measured chemical capacitances and chemical diffusion coefficients in highefficiency DSCs. 33In the following, we provide a brief summary of this model to facilitate comparison with the results from RWNS.Let n 0 be the carrier density in the transport states at the energy E 0 (with total number density N 0 ), and n L the density in localized states defined by eqn (10).These densities can be given with respect to the Fermi level, E F , as The total carrier density is n = n 0 + n L .The transport in extended states is characterized by an effective jump frequency n 0n and constant jump diffusion coefficient D 0 J = D 0 n = D 0 .The central kinetic relationship in the multiple trapping models is the following 31 nhn n i = n 0 n 0n (15)   where hn n i is the average jump frequency for all the carriers.Eqn (15) simply expresses the average number of transitions in the transport levels either in terms of carriers in the transport levels or in terms of all the carriers in the system.It follows from eqn (3) that the jump diffusion coefficient relates to D 0 as The chemical diffusion coefficient is obtained with eqn ( 5) and ( 11): 31 At low Fermi level n E n L , and in terms of the total carrier density, the jump diffusion coefficient is 2.4 Random walk numerical simulation The RWNS simulation method is a stochastic calculation in which electrons are moved at random in a 3 dimensional network of N x Â N y Â N z traps with variable release times.These traps are separated by a L , the lattice parameter or trap separation length.The trap energies are distributed exponentially after eqn (10).Each trap i with energy E i is given a release time, t i , according to 10,34 where r is a random number uniformly distributed between 0 and 1, n 0 (=1/t 0 ) is the attempt-to-jump frequency, and q and V i are the elementary charge and a constant applied bias at site ''i''.Eqn (19) implies that electron detrapping is thermally activated with the trap energy being the activation barrier.It can be shown that the rates for transfer between two neighboured sites i and j (inverse of the release times) fulfils the condition of microscopic reversibility or detailed balance. 35In eqn ( 19) E 1 is the reference level that represents the transport level in multiple trapping or hopping model.This may be different, in general, from the energy E 0 that defines the zero of energies in the distribution (10) as shown below for the hopping model.During the simulation the electrons adopt the release times of the sites that they visit.We call waiting time the difference between the release time of an electron and the time already spent by that electron in its trap.For each simulation move the electron with the shortest waiting time (t min ) is moved to an empty neighboured site chosen at random.This time is then used to advance the total simulation time and to reduce accordingly the waiting times of the rest of the electrons.The simulation is therefore advanced by time increments that depend on the traps occupied at a current simulation time, i.e., the RWNS method is an adaptive time-step simulation procedure.
The calculation is performed by running a number of simulations for different realizations or samples of the trap energy distribution for the network of traps.The final result is averaged over the total number of samples utilized.The more samples are run, the better the statistics.An illustration of the RWNS method employed is shown in Fig. 1.The simulation starts at time t 0 = n 0 À1 by distributing randomly N electrons in the traps.The electrons are then moved between neighboring traps that are placed in an ideal cubic lattice.A move is forbidden if the neighboring trap is already occupied.By this way, electrons diffuse through the lattice of trap sites.In addition, we apply periodic boundary conditions along the three directions of space.Hence, an electron crossing a simulation box boundary is automatically reinjected through the opposite side of the box.Proceeding this way, a stationary state (signaled by a constant electron flux or a constant diffusion coefficient) is rapidly achieved.We will discuss below how this stationary state is approached.
During the simulation, the visits of electrons to traps of energy between E and E + dE are recorded and stored in an histogram N(E).This is the occupancy of the energy levels and N(E)/g(E) is the probability of a trap of energy E to be occupied.In the next section we will show how this probability tends towards the Fermi-Dirac function f(E) as the system approaches the stationary regime.On the other hand the diffusion coefficient is computed from the mean-squared displacement where x, y, and z represent the absolute coordinates of electrons (without periodic boundary conditions).As it will be shown below the mean squared displacement varies linearly with time at longer times (except in conditions close to full occupancy or at very short times).Therefore, it is possible to extract the jump diffusion coefficient via In RW simulations electron mobilities can be extracted from computations executed in the presence of a constant electric field F = (V i+1 À V i )/a L .The mobility is then obtained via where j is the current density.This is computed from a histogram in which the net flux of electrons per unit area and unit time along the field direction is recorded.

Calculation of the jump diffusion coefficient in RWNS
Here we calculate analytically the diffusion coefficient in terms of the parameters of the model.Since the hopping distance is fixed by the lattice parameter, we expect the jump diffusion coefficient obtained from RWNS to have the form Therefore, we can obtain an analytical description of the jump diffusion coefficient by calculating the average jump frequency. 36From eqn (19), the jump frequency from the energy E to the transport level is To obtain a mean jump frequency, eqn ( 24) should be averaged over the distribution of trap energies as follows 37 hni ¼ In this last equation, the factor (1 À f) accounts for the unoccupied target sites.However, this factor is significant only when E F E E 1 , and in the following we consider the situation in which the Fermi level is well below E 1 .The main contributions arise from carriers between E F and E 1 , therefore we can write The Fermi-Dirac distribution can be approximated by a step function below the Fermi level and a Boltzmann factor above the Fermi level.Hence, we obtain for the carrier density Calculating the integral in eqn (26), we obtain Therefore, the jump diffusion coefficient is Let us emphasize that there are two ways to calculate the diffusion coefficient given by eqn (23).First, eqn (25) contains no approximations, therefore numerical integration of eqn (25) should give very accurate results compared to the RWNS, provided that the statistics is sufficiently large.On another hand, eqn (29) gives a closed analytical formula though by using some approximations, therefore some deviations from RWNS calculations may be expected.These results are analyzed below.

Hopping transport
In the hopping model the carriers move by direct transitions between the localized states of the distribution in eqn (10). 38,39he transition probabilities are given by the upward and Electrons are distributed at random among all available traps at time t 0 , they are then moved according to the waiting and released times extracted from eqn (19).Whenever an electron crosses a boundary, it is reinjected through the opposite side so that the total number of electrons remains constant.
This journal is c the Owner Societies 2008 Phys.Chem.Chem.Phys., 2008, 10, 4478-4485 | 4481 downward jump rates where r is the distance between sites, a is the localization radius, and E j , E i , are the energies of the target and starting sites, respectively.The jump diffusion coefficient D J is obtained by averaging the hopping rates of eqn (30) over disordered spatial and energy configurations, in order to obtain the mean jump frequency hni.Since the hopping rates depend exponentially both on energy difference and distance, such averaging is usually very difficult, but it is partially simplified in a system with a steep distribution of localized states, where the hopping process is well described with the concept of transport energy. 38,40The rationale for such approach is that in equilibrium the transport is governed by the fastest hop of a charge carrier.The most probable upward jump corresponds to an optimized combination of the distance and energy difference.
For an exponential distribution of localized levels, the result 41,42 is that the fastest hops occur in the vicinity of the transport energy, given by where 36 independently of the energy of the starting site.Apart from this shift of the transport level with respect to extended states level, the hopping model behaves in a very similar way to the multiple trapping model, provided that the Fermi level remains well below E tr .Averaging the mean square displacement and the jump frequency for the hopping model gives the following result: 36

Steady-state results
The previous calculations indicate that the RWNS method in steady state gives the same Fermi-level dependence of the jump (and chemical) diffusion coefficient as the multiple trapping model.The hopping model is also equivalent to multiple trapping, provided that the transport energy approximation is valid.Therefore RWNS also serves to describe the diffusion coefficient of the hopping model.In each case, appropriate correspondence of parameter is needed in order to match one model to another one.For example in simulations of the hopping model we must set E 1 = E tr , being E 0 (the reference of the DOS) fixed.
In the following, we describe the results of simulations of the diffusion coefficient using the RWNS method.This will help to illustrate the characteristic outcomes of the method and to check the analytical expression derived above.
The simulations were carried out using cubic systems defined by a label N/X 3 where N is the total number of particles (N = 10-400) and X the number of sites along each direction (X = 10-28).The number of simulations (samples) required to get good statistics ranged between 500 and 5000 samples.The results are presented in reduced units of a L (length) and t 0 = n 0 À1 (time).

Diffusion vs. overall electron density
By keeping records of the absolute positions of the electrons along the simulation it is possible to compute the mean square displacement using eqn (3).Results for different electron densities are shown in Fig. 2. First of all, it can be observed that individual electrons diffuse more rapidly for systems with larger overall density.This is nothing else but the well-known trap-filling effect 17,43 and a simple numerical demonstration of the fact that diffusion coefficients are Fermi level dependent in thes kind of systems.In Fig. 2 the diffusivities are plotted versus overall electron density in a log-log graph.It is observed that the RWNS simulation results fit very well to eqn (18).Thus, slopes of 1.88 and 3.47 are obtained for T 0 = 900 (see Fig. 2) and 1400 K (not shown), respectively.According to (18) these should be 2 and 3.67.Note that in real units, D 0 B a 2 L /t 0 and, for a L = 2 nm (N t = 1.25 Â 10 20 cm À3 ) and t 0 = 5 Â 10 À13 s, D 0 = 8 Â 10 À2 cm 2 /s.

Diffusion coefficients vs. Fermi level
In order to extract jump diffusion coefficients as a function of Fermi level RWNS calculations are carried out for different carrier densities.The carrier densities were varied by running calculations corresponding to labels 10/28 3 , 10/24 3 , 10/18 3 , 10/ 15 3 , 10/12 3 , 10/10 3 , 100/28 3 , 200/28 3 and 300/28 3 .The Fermi level is then extracted from the N(E)/g(E) histograms once the steady-state has been reached.Alternatively we have run simulations at a fixed Fermi level by replacing the original exponential distribution by a new distribution with no traps below the Fermi level.This idea is illustrated in Fig. 3.The simulations in this latter case are performed with just one electron.Results for the diffusion coefficient versus Fermi level obtained from both types of calculations are shown in Fig. 4. We can observe that both alternatives give essentially the same results, which shows that collective diffusion is equivalent to single carrier random walk.A difference occurs when we approach full occupancy.In this regime, the diffusion coefficient departs from the linear trend due to the lack of empty jumping sites in conditions of high occupancy (above a 20-30%).This problem does not occur in the one-electron calculations and the diffusion coefficient increases linearly in the log-linear plot up to the ''transport'' level (E 1 = 0 eV).
In Fig. 5, RWNS results are compared with the theoretical predictions of section 2.5.The diffusion coefficients fit nicely to the theoretical formulas (25) and (29).This confirms the validity of the analytical approximations used above.It also shows that the diffusion coefficient can be calculated easily, in the steady state, for arbitrary trap distributions, by the integration procedure indicated in eqn (25).

Electron mobilities and Einstein relation
It is interesting to investigate whether the RWNS reproduces the Einstein relation.In order to do this we run simulations with a bias (F = 3 Â 10 5 V m À1 with a L = 2 nm) and extract the mobilities using eqn (22).We compare them with those calculated with eqn (9) using the jump diffusion coefficients extracted from the mean square displacements at the same conditions but with no applied bias.The electron densities were varied by running four different calculations corresponding to labels 10/28 3 , 10/24 3 , 10/18 3 and 100/28 3 .Results are shown in Fig. 6.We observe that the simulated data are in nice agreement with the classical Einstein relationship for the densities considered in this work.Diffusion coefficients extracted from simulations with applied bias do not comply with the Einstein relation, since in this case D J includes a contribution due to the drift.

Time-dependent results
Having shown the equivalence of RWNS and MT in the steady-state it is very important to discuss the relationship between both methods in transient situations, since the RWNS is a powerful tool to describe transient experiments with rather inhomogeneous carrier distribution. 5,16One can readily think of experimental situations involving the injection and ultrafast detection of electrons into the metal-oxide conduction band,    where no equivalent RWNS description can be obtained, since the latter method lacks the extended states.Therefore, for comparing both methods an initial equilibration of electrons in conduction band with electrons in traps is required in the given experiment.
To illustrate the requirements for application of RWNS, we consider here a classical experiment 44 involving the injection of electrons in the conduction band and monitoring the subsequent dynamics of the free electrons.This experiment is typically described using MT arguments that consider the evolution of electrons in traps over time. 44After the initial pulse, electrons are rapidly captured in traps.The capture process is independent of the trap depth, and the electron density in the bandgap is therefore proportional to g(E).The shallow electrons establish equilibrium with the conduction band, therefore a demarcation level can be defined that sinks with time, due to the fact that increasingly deep levels are able to thermalize their carriers.Consequently, the free electron density is observed to decrease as a power law in time, n c p t T/T0À1 , and this is the mark of the exponential distribution in the bandgap. 44][46] We consider now the transient behaviour as determined from RWNS.First, electrons are randomly distributed initially in the lattice, therefore, in the energy axis the initial distribution follows the exponential shape of the DOS.In order to determine how the electron distribution evolves with time we compute the occupancy histogram n(E) and the occupancy probabilities f(E) at different times.This information is presented in Fig. 7.It can be observed that the electron distribution relaxes as predicted by the standard reasoning (see Fig. 2 in ref. 46) with a well-defined demarcation level at each stage of the evolution, that decreases with time since the electrons above it have been already thermalized.Along this process the slope of the distribution for deep traps remains constant until it matches the density of states.In this situation the system has reached the stationary state and the occupation probability reproduces the Fermi-Dirac distribution.
As mentioned above, in the real experiment the evolution pattern indicated in Fig. 7 will be attained only after a certain time (after the injection of the electrons) for electron capture by traps and initial equilibration of shallow levels with the conduction band.Such initial time, however, is also required for the MT description to be valid, as discussed in ref. 46.Therefore, the transient dynamics is equally well described by MT and RWNS.
Having obtained in Fig. 7 the description of electron relaxation in the energy axis in quasi-equilibrium dynamics, it is interesting to extend the local behaviour to the long range diffusion process.In Fig. 8, the D J value extracted from the slopes of the meansquare displacements at different times is plotted versus the simulation time.A transition from anomalous diffusion to the steady-state, trap-limited transport is observed.At short-times the mean square displacement of electrons is not linear with time (anomalous transport).The diffusion coefficient decreases with time with a power-law before it reaches its stationary value, as discussed in recent work by others. 47e have used the maxima of the electron distributions of Fig. 7 to compute the diffusion coefficient at that value of the (quasi)-Fermi level via eqn (16).We find that, before arriving to the stationary regime, the time-dependent diffusion coefficient matches the ''equilibrated'' value predicted by the MT formalism.This shows that the condition of local quasi-equilibrium allows one to predict the fast transient dynamics from steady-state values that can be obtained analytically, as discussed in section 3.

Conclusions
The RWNS method is a computationally simple and versatile simulation procedure that can be used to study transport in systems which a broad distribution of trap energies.It has been proved that it reproduces the main predictions of the multipletrapping and hopping transport models in steady state.We have also seen that a quasi-stationary situation is adopted by the electron ensemble when the system is approaching this steady state.The RWNS method can be used to simulate MT and hopping transport in a variety of situations: transient photovoltages and photocurrents, and to analyse the effect of different trap distributions and geometries.

Fig. 1
Fig.1Illustration of a RWNS calculation.Electrons are distributed at random among all available traps at time t 0 , they are then moved according to the waiting and released times extracted from eqn(19).Whenever an electron crosses a boundary, it is reinjected through the opposite side so that the total number of electrons remains constant.

Fig. 2
Fig. 2 Upper panel: mean square displacement as computed via eqn (20) in a RWNS simulation with exponential distribution of trap energies.Calculations were performed for T 0 = 900 K, T = 300 K and comprised 5000 samples.Lower panel: jump diffusion coefficients as obtained from the slopes of the mean square displacement versus overall electron density.

Fig. 3
Fig. 3 Scheme of the trap energy distributions used in the multielectron and single-electron RWNS calculations.The equivalence between both is shown explicitly.

Fig. 4
Fig. 4 Jump diffusion coefficient as a function of Fermi level for both multi-electron and single-electron RWNS simulations.The Fermi level in multi-electron calculations is extracted from fitting the occupancy curves to a Fermi-Dirac function.

Fig. 7
Fig. 7 Electron occupancy distributions as extracted from RWNS simulations at different times.The distribution of trap energies (DOS) and the probability function in the steady-state (f(E)) are also included in the graph.

Fig. 8
Fig. 8 Time evolution of the jump diffusion coefficient as extracted from RWNS calculations (solid line) and prediction of the multipletrapping theoretical formula of eqn (29).