Exact statistical solution for the hopping transport of trapped charge via finite Markov jump processes.

In this study, we developed a discrete theory of the charge transport in thin dielectric films by trapped electrons or holes, that is applicable both for the case of countable and a large number of traps. It was shown that Shockley-Read-Hall-like transport equations, which describe the 1D transport through dielectric layers, might incorrectly describe the charge flow through ultra-thin layers with a countable number of traps, taking into account the injection from and extraction to electrodes (contacts). A comparison with other theoretical models shows a good agreement. The developed model can be applied to one-, two- and three-dimensional systems. The model, formulated in a system of linear algebraic equations, can be implemented in the computational code using different optimized libraries. We demonstrated that analytical solutions can be found for stationary cases for any trap distribution and for the dynamics of system evolution for special cases. These solutions can be used to test the code and for studying the charge transport properties of thin dielectric films.

Introduction

New dielectrics, like high- , , [1-5] and low- ones[6-9], are extremely perspective materials for using them in modern and future electronic devices, including memory and digital ones. That is why the knowledge about the charge transport processes and mechanisms in dielectrics is critical for modern microelectronics. High- dielectrics are promising materials for resistive random access memory (RRAM), ferroelectric random access memory (FRAM) and memristor devices, while low- dielectrics are the insulating dielectrics that separate wire interconnects and transistors from each other in high-speed integrated circuits[6]. High--based MOSFETs, FRAM and flash memories, and low- interconnect separators need dielectrics with low leakage currents, while RRAM devices exhibit a better performance when the used dielectric layers are mostly nonstoichiometric, i.e. include a lot of defects and traps[10-12]. That is why producing high-quality electronic devices requires a strict control of deposited dielectric film quality to monitor leakage currents.

The charge currents in metal-insulator-semiconductor (MIS) and metal-insulator-metal structures (MIM) might be limited by a contact metal-insulator or insulator-semiconductor (contact-limited currents) or by the traps in the insulator bulk (bulk- or trap-limited currents). Contact-limited current models describe the electron and hole emission from semiconductor or metal electrodes to traps (localized states in the bandgap) or to the conduction band in dielectrics. The quantum field emission, first proposed by R. Fowler and L. Nordheim, describes the charge-carriers tunnelling through a triangle energy barrier to the dielectric conduction band and does not depend on temperature[13]. At high temperatures, a significant contribution to the transport might be given by the thermally assisted tunnelling[14,15] and field-enhanced thermionic emission of hot carriers (Schottky effect)[16]. The direct tunnelling from one electrode to another through a trapeze barrier takes place in the structures based on ultra-thin ( nm) dielectric films[17].

Bulk-limited charge transport models can be divided into two large subgroups: isolated trap ionisation (Frenkel model[18,19], multiphonon single trap ionization proposed by Makram-Ebeid and Lannoo[20]) and hops between neighbouring traps (thermally assisted traps between overlapped Coulomb centres[21,22], phonon-assisted tunnelling of electrons (holes) between neighbouring traps[23]).

In essence, all these models give the probabilities of a hop between two traps and between a trap and an electrode. Usually, a set of approaches was used to describe the charge transport between traps in dielectric layers using these probabilities. One on them is the time-dependent direct or Monte-Carlo simulations with hops between traps, and it requires carrying out many simulations to obtain sufficient statistics for data averaging in order to identify general patterns. Another approach is discrete Shockley–Read–Hall-like (SRH) 1D transport equations with the boundary conditions[23,24]:where is the average filling of the th trap, t is time, is the probability rate of a hop from the th trap to the th one, is the trapped carrier generation rate (due to free carriers localization) and is the trapped carrier recombination rate (due to the trap ionization to a conduction/valence band, recombination with free carriers with an opposite sign) and N is the trap number. In case of low free charge carriers density (here is the trapped carrier density, n if the free carriers density), the terms and may be omitted due to their smallness. Recently, discrete equation (1) has been extended to a continuous model[25]:where a is the distance between traps, x is the coordinate, are the hop rates between traps along and against electric field F. This approach is approximate, but it works fine for a relatively large number of traps in a dielectric medium. Modern electronic devices based on thin films are produced with the feature size of several nanometers. In this case, the number of traps in dielectric bulk N is countable, and the filling averaging over them might give quite large deviations due to the law.

All these methods are approximate due to different reasons. Monte-Carlo simulations are approximate by the virtue of them being a numerical calculation method. Both discrete and continuous SRH strongly rely on the fact that fillings of neighboring traps are statistically independent. In fact, they are obviously not independent. That is why the SRH approach is an approximate one. Furthermore, time-dependent SRH equations should be integrated numerically, and it adds an error to the result. The approximate nature of all these methods require using a validation by, for example, some exact analytical solutions as a reference. To obtain such solutions a new model needs to be developed. The aim of the present work is the development of an exact universal model of trapped charge transport via a hopping between traps in strong electric fields applicable both for the case of countable and a large number of traps.

Results

Markov-jump-like (MJ) model formulation

Let us consider the volume of a dielectric containing N traps. From the classical standpoint (considering electrons to be strongly localized inside traps) there is a finite set of different variants of charge localization inside traps (states): . Each state has a certain probability of transition to state during a short period of time , which we will indicate as . Obviously, the probability of not changing state can be obtained from the probabilities of all transitions from : . Thus, the transition probability matrix sized can be compiled.

At each moment in time, the system can be characterized by the column vector of probabilities, where is the probability of finding the system in state . The average state evolution of the system during one time step can now be described by the expression , which can be transformed to the form , where is the identity matrix. The differential form of describing the evolution of can be achieved by considering the limit :This equation is known as a Kolmogorov backward equation (KBE). The element of matrix can be treated as the rate of transition from state to state (so, is the probability of such a transition in short time ). Diagonal elements can again be found from all of the rates of transition from the state : . The resulting differential equation (3) can be solved analytically:If we represent the state as a column vector , where is the filling of the th trap in state (takes the value 0 if the trap is empty and the value 1 if a charge carrier is captured in the trap), then the average ensemble filling of traps can be expressed as an average for all states according to their individual probability values:where is the matrix of states such that is the occupancy of the th trap in the th state.

Each option of filling traps is characterized by the matrix of currents , where is the current between the th and th traps in state . Obviously, the matrix is skew-symmetric: . So, the matrix of currents can be expressed as an average for all states according to their individual probability values:The current between electrodes and traps can also be introduced in this matrix by adding the corresponding rows and columns to the matrix. Indeed, we can treat each electrode as an additional trap with the corresponding currents towards/from all of the other traps. For example, one can treat the left and right electrode as th and th traps, correspondingly, so that and are the currents from the left electrode to the th trap and the current from the th trap towards the left electrode; and are similar currents for the right electrode.

Stationary state

In many cases, the desired characteristic of a system is the stationary state achieved with time and not the dynamics of its development. In a stationary state, the column vector stops changing. Therefore, from (3) it follows that . Also, by virtue of its definition, the probability vector must always satisfy the additional constraint that the sum of its elements is equal to 1. This means that to find a stationary vector of probabilities, it is necessary to solve a system of linear equations:It can be solved as is or can be converted to a more standard form. Obviously, the system is overdetermined (the number of equations is larger than the number of the ‘unknown’). However, the conventional square non-homogeneous system can also be produced from it in a number of ways. One way is to use the result of adding the first equation to the rest equations:The system of equations can be expressed in the matrix form:where is the matrix of ones, is the vector of ones. Now, using the matrix notation, the stationary state characteristics can be expressed explicitly:It is important to note that, despite the mathematical elegance, this way of expressing the solution may not be the most optimal one. Depending on the numerical methods in use, it might be wise to keep the sparsity of matrices. The sparse square non-homogeneous form can be achieved, for example, by simply removing the second equation from the system (5). Indeed, as far as , one could see that the sum of the last equations in (5) is equivalent to the second equation. So, the second equation can be discarded as linearly dependent on other equations. It should be noted that the same is true for any equation in (5) except for the first one. Now the matrix form of this system of equations is:where is the matrix with the elements of the first row substituted with ones, is the vector, such that the first element is equal to 1 and the rest of the elements are equal to 0. The explicit form of the stationary characteristics in this case is:It can also be seen that column vector can be found as the first column of matrix , and column vector can be found as the first column of matrix (due to the associative property of matrix multiplication).

Hop current model

The elementary acts of hopping conduction in the model are:The probability of observing an elementary act during the time is the product of and the rate of this act. Thus, the probability of observing the event consisting of L elementary acts during the time will be proportional to . This means that an addition to matrix from all events consisting of more than one elementary act does not affect the matrix according to its definition. Indeed, the limit in (3) gives 0 for any addition with . So, considering the events taking place during the transition from one state to another, we can limit ourselves just to the events consisting of one elementary act. However, some transitions are not achievable by only one elementary act. For example, if the state differs from the state in more than 2 places, the rate of such transition is zero (). The number of possible non-zero elements and their fraction in matrix can be easily calculated:Thus, even the most general form of matrix can be considered sparse (contains mostly zeros) for a reasonably large number of traps, as shown in Fig. 1.

Figure 1

Illustration of sparsity of matrix for , 5 and 6. White corresponds to zero-valued elements (impossible transitions).

hopping of one charge carrier from an filled trap to an empty one. The rate of such event will be indicated as ;

extraction (ionization and departure) of one carrier from a filled trap towards an electrode. The rate of such event will be indicated as ;

injection of one carrier from an electrode into an empty trap. The rate of such event will be indicated as .

Many of the indicated possible non-zero elements can also be treated as zeros due to the fact that the probability of charge transitions is extremely dependent on the distance and the hopping of a trapped charge carrier into a trap not adjacent to it (far jump) is unlikely, compared to a jump into the neighboring one. The same applies to the carrier injection from the electrode into an empty trap and extraction from an occupied trap towards the electrode. So, for the simplicity of computations, one can limit oneself to considering only short-range events.

Examples

In this section we will apply the model to some simple cases of the hopping conduction to illustrate the general mathematical formulation.

Two traps in a strong electric field

Consider a simple case of a 2-trap chain under the strong right-to-left electric field, as illustrated in Fig. 2. Under such conditions, we can neglect the probability of electrons moving to the left. In this case, only 3 short-range events can happen: injection from the left electrode into the first trap, jumping of the electron trapped in the first trap towards the second trap and extraction of the electron trapped in the second trap towards the right electrode. In this case the considered matrices are as follows:

Figure 2

Illustration of different states and electron transition probability values for under the strong right-to-left electric field.

where q is the signed carrier charge (equal to the elementary charge modulo).

Solving equations (5) using (7) gives us the stationary probability values of different states:The average stationary filling can be calculated using the trap filling in different states, and the found stationary probabilities of these states (9) areThe stationary current matrix can be calculated using the state currents (8), and the found stationary probabilities (9) areHerein, Shockley–Read–Hall-like transport equations (1) describing the dynamics of average occupancy areThe stationary state is given by the following nonlinear system of equations:This nonlinear system of equations can be solved analytically:According to the SRH model, the stationary current can be found using average filling:

Bulk-limited current (, )

Consider a simple case of a uniform linear N-trap chain under a strong right-to-left electric field. Also consider that the current is bulk-limited when and , i.e. the interfaces between electrodes and a dielectric layer are transparent for charge carriers. In this case, the first trap is filled, and the last one is empty. It is easy to see that this case is identical to the case of the uniform linear chain of traps under a strong electric field when (the always filled first trap can be treated as the left electrode with and the always empty last trap can be treated as the right electrode with ).

Case of

This case is identical to the case of a two-trap chain (the two middle traps) with the corresponding matrices:To obtain the dynamic solution, the initial conditions should be given. Let us study the dynamic of filling the initially empty middle traps. In this case and equation (4) gives the following evolution of state probability values:The dynamic solution for the average trap filling can be found using the probability values:The time dependence of currents between traps on time is as follows:The evolution of currents between traps , and is shown in Fig. 3. One can see that, at the beginning, the current between the left traps is quite large, while the currents between the right traps are small. In a time , a stationary state is established, and the currents become the same . It should be noted that we do not take into account the quantum nature of an electron, because all processes under consideration take place at room temperature and above. Also, electron quantum states are destroyed during inelastic (de-)trapping acts and the transport is not ballistic. That is why we do not get the conductivity quantum in the obtained solutions.

Figure 3

Evolution of currents between traps , and for a uniform chain of four traps under strong electric field.

It is important to note that matrices and , in the case of the linear 4-trap chain under the strong right-to-left electric field with a bulk-limited current, are simply proportional to . Obviously, this will be the case for the N-trap chain also. A couple of significant conclusions can be drawn from this fact:Here and are the coefficients that depend on the trap number N only. The values of and are shown in Fig. 4. One can see that tends to the value 0.25 as the number of traps increases for both SRH and MJ cases, but the SRH model predicts a lower current than MJ model. Also, the Markov jump processes approach predicts longer transients (steady-state establishment) than the Shockley–Read–Hall-like approach.

Figure 4

The and coefficients values.

stationary average trap filling does not depend on the value;

stationary current will have the form ;

time of establishing a near-stationary solution will have the form .

The comparison of stationary trap fillings and the current in the case of 4-trap chain for the MJ and the SRH models is shown in Table 1. One can see that the Shockley–Read–Hall-like model gives deviations up to from the exact values that are calculated in terms of the finite Markov jump processes approach. On the one hand, the obtained deviations of the current might be smaller than different noises in real experiments, but they might be too large for simulations of high-quality electronic devices.

Table 1

Comparison of analytically found stationary trap filling distributions and currents between traps for MJ and SRH models.

	Markov jump-like	Shockley–Read–Hall-like	Error (%)
n¯1st	1	1	0
n¯2st	35=0.6	5-12≈0.618	≈3
n¯3st	25=0.4	3-52≈0.382	≈-5
n¯4st	0	0	0
i^st	25qH=0.4qH	qH3-52≈0.382qH	≈-5

The same procedure can be carried out for larger trap chains. The comparison of stationary trap filling distributions calculated within the SRH and the MJ models for trap chains with , 8 and 12 is shown in Fig. 5. One can see that the shape of all curves is close and it tends, with the increase in N, to that shown in Figure 3 in[25] for the continuum case.

Figure 5

Trap filling distributions, calculated in terms of SRH and MJ approaches for , 8 and 12.

Discussion

We have developed a discrete model of the charge transport in dielectrics in terms of finite Markov jump processes, and it takes into account hops of localized charge carriers between traps in a dielectric medium in strong electric fields with the injection from and extraction to electrodes (contacts). The developed model is formulated as a matrix differential equation or a system of linear algebraic equations. The dynamic solutions describing the charge transport are found using a matrix exponential. The stationary solutions describing the charge transport are found by solving a system of linear equations. The major drawback of the presented method for the charge current description is the dramatic rise in the computational complexity with the increasing amount of trap. However, an electronic system, containing a dozen of traps, can be easily simulated using a laptop PC. Moreover, we demonstrated that the solutions can be acquired by manipulations with sparse matrices. This means that the numeric computation of transition processes can be solved quite quickly using optimized algorithms[26,27]. The developed model is universal and can be applied to a relatively large number of traps distributed in 1D, 2D or 3D systems. Also, the matrix sparsity can be increased by ignoring extremely long hops in a specific spatial trap distribution, thus, reducing the total difficulty of the task.

The analytical stationary solutions, that are close to experimentally achievable conditions, are found. The comparison with other theoretical models[25] shows a good agreement, which proves their good applicability, at least, for the stationary problems under examination. Also, the presented model allows one to get exact analytical solutions of the charge transport evolution in time, while Shockley–Read–Hall-like transport equations require their integration to get the same results. This procedure is more difficult: it requires more computational resources and time. So, the advantage of the presented solutions over others lies in the fact that it can be used as an exact reference to test other (approximate) solutions in both dynamic and stationary cases.

The developed model, formulated in linear algebraic equations, can be implemented in the computational code using different libraries, such as BLAS, LAPACK, MKL, SciPy/NumPy etc. Particularly, the proposed approach can enhance the current implementation of the leakage current simulations through dielectric layers in the Synopsys TCAD software package[28]. It was shown that the stationary current-voltage characteristics of metal-insulator-semiconductor and metal-insulator-metal structures can be calculated analytically for any trap distribution. Also, the model implemented as a numeric code can be used in simulations of flash, FRAM, RRAM memories and memristor resistive switching with different spatial designs, including vertical and planar ones.