Markov jump transport of trapped charge in thin dielectric layers

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 1D transport through dielectric layers, might incorrectly describe the charge flow through the ultra-thin layers with a countable number of traps, taking into account 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-, twoand 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 dynamics of system evolution for special cases. These solutions can be used to test the code and for studying of charge transport properties of thin dielectric films. Introduction New dielectrics, like high-κ HfO2, Ta2O5, ZrO2 and low-κ ones6–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 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 circuits6. High-κ-based MOSFETs, FRAM and flash memories and low-κ interconnect separators need dielectrics with low leakage currents, while RRAM devices exhibit better performance when used dielectric layers are mostly nonstoichiometric, i.e. include a lot of defects and traps10–12. That is why producing high-quality electronic devices requires hard control of deposed dielectric films 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 (bulkor trap-limited currents). Contact-limited current models describe 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 charge carriers tunnelling through a triangle energy barrier to the dielectric conduction band and does not depend on the temperature13. At high temperatures a significant contribution to the transport might give thermally assisted tunnelling14, 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 structures based on ultra-thin (< 5 nm) dielectric films17. Bulk-limited charge transport models can be divided into two large subgroups: isolated trap ionisation (Frenkel model18, 19, multiphonon single trap ionization proposed by S. Makram-Ebeid and M. Lannoo20) and hops between neighbour traps (thermally assisted traps between Coulomb overlapped centres21, 22, phonon-assisted tunnelling of electrons (holes) between neighbouring traps23). In essence all these models give the probabilities of 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 time-dependent direct or Monte-Carlo simulations with hops between traps, and it requires carrying out many simulations to obtain sufficient statistics for averaging in order to identify general patterns. Another approach is discrete Shockley-Read-Hall-like (SRH) 1D transport equations with the boundary conditions23, 24: dn̄α dt = n̄α−1(1− n̄α)Hα−1,α − n̄α(1− n̄α−1)Hα−1,α + n̄α+1(1− n̄α)Hα+1,α − n̄α(1− n̄α+1)Hα,α+1 +Gα −Rα ; dn̄1 dt = (1− n̄1)I1− n̄1E1 + n̄2(1− n̄1)H2,1− n̄1(1− n̄2)H1,2 +G1−R1; dn̄N dt = (1− n̄N)IN− n̄NEN + n̄N−1(1− n̄N)HN−1,N− n̄N(1− n̄N−1)HN,N−1 +GN−RN , (1) where n̄α is the average filling of the α-th trap, t is the time, Hα,β is the probability rate of a hop from the α-th trap to the β -th one, Gα is the trapped carrier generation rate (due to free carriers localization) and Rα is the trapped carrier recombination rate (due to 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 nt n (here nt is the trapped carrier density, n if the free carriers density), the terms Gα and Rα may be omitted due to their smallness. Recently, discrete equation (1) has been extended to a continuous model25: ∂ n̄ ∂ t =−a ∂ ∂x {[ H+−H− ] n̄(1− n̄) } + a2 2 ∂ ∂x {[ H++H− ]∂ n̄ ∂x }


Introduction
New dielectrics, like high-κ HfO 2 , Ta 2 O 5 , ZrO 2 1-5 and low-κ ones [6][7][8][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 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 better performance when used dielectric layers are mostly nonstoichiometric, i.e. include a lot of defects and traps [10][11][12] . That is why producing high-quality electronic devices requires hard control of deposed dielectric films 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 (bulkor trap-limited currents). Contact-limited current models describe 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 charge carriers tunnelling through a triangle energy barrier to the dielectric conduction band and does not depend on the temperature 13 . At high temperatures a significant contribution to the transport might give 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 structures based on ultra-thin (< 5 nm) dielectric films 17 .
In essence all these models give the probabilities of 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 time-dependent direct or Monte-Carlo simulations with hops between traps, and it requires carrying out many simulations to obtain sufficient statistics for 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 : wheren α is the average filling of the α-th trap, t is the time, H α,β is the probability rate of a hop from the α-th trap to the β -th one, G α is the trapped carrier generation rate (due to free carriers localization) and R α is the trapped carrier recombination rate (due to 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 n t n (here n t is the trapped carrier density, n if the free carriers density), the terms G α and R α 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, H ± = H(±F) are the hop rates between traps along and against electric field F. This approach is approximate, but it works fine for a large number of traps in a dielectric medium. Modern electronic devices are based on thin films 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 √ N/N = 1/ √ N law. 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.

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 2 N different variants of charge localization inside traps (states): n n n 1 , n n n 2 , . . . , n n n 2 N . Each state n n n α has a certain probability of transition to state n n n β during a short period of time δt, which we will indicate as P α,β . Obviously, the probability P α,α of not changing state n n n α can be obtained from the probabilities of all transitions from n n n α : P α,α = 1 − β =α P α,β . Thus, the transition probability matrixP sized 2 N × 2 N can be compiled.
At each moment in time, the system can be characterized by the column vector p p p = [p 1 ; p 2 ; . . . ; p 2 N ] of probabilities, where p α is the probability of finding the system in state n n n α . The evolution of average state of the system during one time step δt can now be described by the expression p p p(t + δt) =Pp p p(t), which can be transformed to the form p p p(t + δt) − p p p(t) = P −Î p p p(t), whereÎ is the identity matrix. The differential form of describing the evolution of p p p can be achieved by considering the limit δt → 0: This equation is known as a Kolmogorov backward equation (KBE). The element Q α,β of matrixQ can be treated as the rate of transition from state n n n α to state n n n β (so, Q α,β δt is the probability of such a transition in short time δt). Diagonal elements can again be found from all of the rates of transition from the state n n n α : Q α,α = − β =α Q α,β . The resulting differential equation (3) can be solved analytically: If we represent the state n n n γ as a column vector n n n γ = n whereN is the matrix of states such that N α,β is the occupancy of the α-th trap in the β -th state.
Each option of filling traps n n n γ is characterized by the matrix of currentsî γ , where i γ α,β is the current between the α-th and β -th traps in the state γ. Obviously, the matrix is skew-symmetric: i γ α,β = −i γ β ,α . 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 into this matrix by adding the corresponding rows and columns to the matrix. Indeed, we can treat each electrode as an additional trap with corresponding currents towards/from all of the other traps. For example, one can treat the left and right electrode as (N + 1)-th and (N + 2)-th traps, correspondingly, so that i

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 p p p stops changing. Therefore, from (3) it follows thatQp p p st = 0. Also, by virtue of its definition, the probability vector p p p 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 2 N + 1 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 'unknowns'). 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 2 N equations: The system of equations can be expressed in the matrix form: whereĴ is the matrix of ones, j j j is the vector of ones. Now, using the matrix notation, the characteristics of stationary state 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). Now the matrix form of this system of equations is: whereQ 1 is the matrixQ with the elements of the first row substituted with ones, j j j 1 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 stationary characteristics in this case is: It can also be seen that column vector p p p st can be found as the first column of matrixQ −1 1 , and column vectorn n n st can be found as the first column of matrixNQ −1 1 (due to the associative property of matrix multiplication).

Hop current model
The elementary acts of hopping conduction in the model are: • Hopping of one charge carrier from an filled trap to an empty one. The rate of such event will be indicated as H; • Extraction (ionization and departure) of one carrier from a filled trap towards an electrode. The rate of such event will be indicated as E; • Injection of one carrier from an electrode into an empty trap. The rate of such event will be indicated as I.
The probability of observing an elementary act during the time δt is the product of δt and the rate of this act. Thus, the probability of observing the event consisting of L elementary acts during the time δt will be proportional to δt L . This means that an addition to matrixP from all events consisting of more than one elementary act does not affect the matrixQ according to its definition. Indeed, the limit in (3) gives 0 for any addition with L > 1. So, considering the events taking place during the transition from one state to another, we can limit ourselves just to events consisting of one elementary act. However, some transitions are not achievable by only one elementary act. For example, if the state n n n α differs from the state n n n β in more than 2 places, the rate of such transition is zero (Q α,β = 0). The number of possible non-zero elements and their fraction in matrixQ can be easily calculated: Thus, even the most general form of matrixQ can be considered sparse (contains mostly zeros) for a reasonably large number of traps, as shown in Figure 1. 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 Figure 2. Under such conditions, we can neglect the probability of electrons moving to the left. In this case, only 3 short-ranged 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 considered matrices are as follows: 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: Average stationary filling can be calculated using the trap filling in different states and found stationary probabilities of these states (9)

5/10
The stationary state is given by the following nonlinear system of equations: − En st 2 = 0. 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 (n 1 = 1,n N = 0) Consider a simple case of uniform linear N-trap chain under a strong right-to-left electric field. Let us also consider, that the current is bulk-limited when I left H E left and E right H I right , 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 uniform linear chain of N − 2 traps under a strong electric field when I = E = H (always filled first trap can be treated as the left electrode with I = H and always empty last trap can be treated as the right electrode with E = H).

Case of N = 4
This case is identical to the case of two-trap chain (the middle two traps) with the corresponding matrices: To acquire the dynamic solution, the initial conditions should be given. Let us study the dynamic of filling the initially empty middle traps. In this case p p p(0) = [1; 0; 0; 0] and equation (4) gives the following evolution of state probability values: The evolution of currents between traps i 12 , i 23 and i 34 is shown in Figure 3. One can see that, at the beginning, the current between left traps is quite large, wile the currents between right traps are small. In a time t ≈ 5/H, a stationary state is established, and the currents become the same i = 0.4qH. 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 being 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.
It is important to note that matricesQ andî γ in the case of linear 4-trap chain under the strong right-to-left electric field with bulk limited current are simply proportional to H. Obviously, this will be the case for the N-trap chain also. A couple of significant conclusions can be drawn from this fact: • stationary average trap filling does not depend on H value; • stationary current will have the form C i qH; • time of establishing a near-stationary solution will have the form C τ /H.
Here C i and C τ are the coefficients that depend on the trap number N only. The values of C i and C τ are shown in Figure 4. One can see that C i tends to the value 0.25 as the number of traps increases for both cases SRH and MJ, but SRH model predicts lower current than Markov jump. Also, Markov jump approach predicts longer transients (steady-state establishment) than the Shockley-Read-Hall-like approach.
Comparison of stationary trap fillings and the current in the case of 4-trap chain for MJ and SRH models is shown in Table 1. One can see, that the Shockley-Read-Hall-like model gives deviations up to 5% from exact values that are calculated in terms of Markov-jump approach. On the one hand, 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.
The same procedure can be carried out for larger trap chains. Comparison of stationary trap filling distributions calculated within SRH and MJ model for trap chains with N = 4, 8 and 12 are shown in Figure 5. One can see that the shape of all curves is close and it tends, with the increase on N, to that shown in Figure 3 in Ref. 25 for the continuum case.

Discussion
We have developed a discrete model of the charge transport in dielectrics in terms of Markov jump, and it takes into account hops of localized charge carriers between traps in a dielectric medium in strong electric fields with 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. 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 large number of traps distributed in 1D, 2D or 3D systems.
The analytical stationary solutions, that are close to experimentally achievable conditions, are found. The comparison with other theoretical models shows a good agreement 25 . The developed model, formulated in linear algebraic equations, can be implemented in the computational code using different libraries like BLAS, LAPACK, MKL, SciPy/NumPy etc. Particularly, proposed approach can enhance the current implementation of the leakage current simulations through dielectric layers in Synopsys TCAD software package 28 . It was shown that the stationary current-voltage characteristics of metal-insulatorsemiconductor 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.