Charge transport mechanism in dielectrics: drift and diffusion of trapped charge carriers

In this study, we developed a continuum theory of the charge transport in dielectrics by trapped electrons and holes, which takes into account two separate contributions of the current of trapped charge carriers: the drift part and the diffusion one. It was shown that drift current is mostly dominant in the bulk, while the diffusion one reaches significant values near contacts. A comparison with other theoretical models and experiments shows a good agreement. The model can be extended to two- and three-dimensional systems. The developed model, formulated in partial differential equations, can be numerically implemented in the finite element method code.

The knowledge about charge transport processes and mechanisms in dielectrics is critical for modern microelectronics, because many dielectrics, like SiO 2 and high-κ (for example HfO 2 ) are used as gate dielectrics 1,2 and low-κ ones-as insulating dielectrics separating the wire interconnects and transistors from each other in highspeed integrated circuits 3 , as well as blocking insulators in silicon-oxide-nitride-oxide-silicon-type (SONOS) 4,5 and TaN-high-κ-nitride-oxide-silicon-type (TANOS) 6,7 flash memory devices. High-κ dielectrics are promising candidates to be used as functional materials in resistive random access memory devices (RRAM) and memristors, that would involve integrating the most favourable properties of both rapid dynamic random access memory (DRAM) and non-volatile solid drive memory (SSD) [8][9][10][11][12] . Therefore, MOSFETs and TANOS need high-κ and low-κ dielectrics featuring low leakage currents, but RRAM cells require a dielectric medium that exhibits reversible resistive switching. Controlling the process of dielectric film synthesis to manage leakage currents can help creating high-end components for different devices.
In general, the charge transport mechanisms in dielectric layers can be conditionally divided into two groups: contact-limited and bulk-limited via traps ones. The first group is the contact-limited models that describe electrons and holes emission from metal or semiconductor contacts to the conduction band or traps in dielectric layers. The field emission, also known as Fowler-Nordheim model, describes charge carriers tunnelling through a triangle energy barrier to the dielectric conduction band 13 . This is the quantum effect that does not depend on the temperature. As far as the Fermi energy level, as well as the electron energy distribution in metals and semiconductor contacts, depends on temperature, thermally assisted tunnelling might give a significant contribution to the current 14,15 . At the first stage, charge carriers are excited to a certain energy due to the phonon absorption, and then they tunnel through the triangular barrier. In case of high temperatures, the field-enhanced thermionic emission (Schottky effect) of hot carriers takes place 16 . A combination of Fowler-Nordheim tunnelling with the Schottky effect gives the so-called Schottky-Nordheim barrier, which is the barrier model used in deriving the standard Fowler-Nordheim-type equation 17 . In case of ultrathin ( < 5 nm ) dielectric films, charge carriers might exhibit their direct tunnelling from one electrode to another through a trapeze barrier 18 . If the structure is based on a thick ( > 5 nm ) dielectric layer with traps, then the electrons and holes might be injected within the trapassisted tunnelling (TAT) mechanism 19 . At the first stage, charge carriers tunnel from the first electrode to the trap in the dielectric film bulk, and then they tunnel from the trap to the conduction (valence) band in the next electrode. Such a multi-stage process may be more profitable than direct tunnelling or the tunnelling according to the Fowler-Nordheim mechanism.
The first theoretical bulk-limited charge transport model was proposed by Frenkel 20,21 . This model is based on a single charged Coulomb centre (trap). In a strong electric field the Coulomb barrier is decreased, and a localized electron (hole) can leave the trap through three channels: by tunnelling through a reduced energy

Results
Direct simulations. For simplicity, the tunnelling current through dielectrics can be reduced to the tunnelling-driven movement of charges along the multiple independent chains of traps. The probability of tunnelling between traps heavily depends on the spacing between them and electric field and may be different for each pair. Considering only the tunnelling between the nearest traps and assuming tunnelling to distant traps to be highly improbable, the evolution of a chain of traps can be described by the following equations: where P + i is the probability of charge trapping on the ith trap in an infinitesimal amount of time δt ; P − i is the probability of charge leaving the ith trap in an infinitesimal amount of time δt ; n i is the occupancy of the ith trap ( n i = 1 if a trap is occupied with a charge carrier and n i = 0 if a trap is empty); Tun ± i,j is the probability rate of tunnelling between the ith trap and jth trap (superscripts ' + ' and '−' represent forward and backward tunnelling, respectively). The evolution of M-trap chain endpoint traps is described by the following equations: where Inj i is the probability rate of charge injection from the electrode into the ith trap, and Ion i is the probability rate of ionization and subsequent removing into an electrode of a charge carrier from the ith trap.
Considering the occupancy of neighbouring traps to be independent random values, the evolution of this value is given by the following equation: where n i is the average occupancy of the ith trap. The evolution of average occupancy of the traps, placed at the chain endpoints is described by the following equations:  Table 1.
One can see that both approaches give the solutions that are in a good agreement with each other. The tunnelling rate is considered to be the same for any pair of neighbouring traps so that Tun + i,i+1 = P . Hereinafter, t 0 = 1/P is the latency time of tunnelling from an occupied trap to an empty neighboring one.
In the case of bulk-limited current, Eq. (4) can be rewritten simply as: n 1 = 1 , n M = 0 . This scheme is applicable as it is and it was used by different authors [29][30][31] , but it has a number of limitations. First of all, it scales badly with the increase in the number of traps (computing the scheme for a large number of traps can require great computational power). Also, the use of analytical considerations and conventional numerical analysis dictates the need for the formulation of the problem in the form of a system of partial differential equations.
Continuum model. Probability rates Inj i and Ion i depend on the temperature and electric field at the ith trap and can be considered a function of the ith trap coordinate x i . By the same logic, Tun ± i,i+1 can be considered www.nature.com/scientificreports/ a function of coordinate (x i+1 + x i )/2 halfway between traps. By considering n to be a continuous function of time and coordinate and using the second order Taylor expansion of n i+1 , n i−1 , , the following partial differential equation can be obtained after factoring resulting members: where x is the coordinate; Tun ± = Tun(±F) , F is the electric field and a is the distance between traps. The boundary conditions are determined by the current continuity at the boundaries and can be described by the following equations: The dependence of functions Inj and Ion on the values of temperature and electric field at a given coordinate. Current density values can be obtained using the charge density continuity equation and written as: where N is the trap density, q is the elementary charge (absolute value). Due to the structure of Eq. (7), it is convenient to treat the first part of the right-hand side of the equation as a 'drift' current j drift and the second part-as a 'diffusion' current j diff . It is important to note that, although the diffusion current is proportional to the higher order of small parameter a, it cannot be omitted. This fact can be easily seen from the following illustrative example. In the case of uniform tunnelling probability for every trap, Eq. (5), obviously, has no nontrivial stationary solution. Contrariwise, the solution of the full form of Eq. (5) can be found analytically. For example, in the case of bulk-limited current ( n(0) = 1 , n(d) = 0 , where d is the chain length), the solution is: where is the dimensionless number that characterizes the trapped charge spatial distribution in the stationary state; x = x/d is the dimensionless coordinate; β(χ) is the implicit function which satisfies the following equation: The solution (8) now can be used to obtain stationary current values (which are, in most cases, the final objective of the study):  1) and (2) for 500 chains and by calculating Eqs. (3) and (4). where k is the transcendental equation root: The constant B and the functions f 0 (k) , f 1 (k) are following: The solution (9) now can be used to obtain stationary current values:   www.nature.com/scientificreports/ Despite its analytical complexity, the obtained solution (9) and (10) can be indispensable for verification of various numerical codes. Since, the solution (8) is more simple, it can be analysed as follows. In Fig. 3 is a comparison of obtained numerical solutions using different approaches at different moments for the 30-trap chain. A scheme of linear differences is presented by blue symbols, the solution of differential equation (5) is shown by orange lines and the steady-state solution (8) is shown by a red line. The case of bulk-limited current is considered. It should be noted, that the non-stationary solution becomes reasonably close to the steady-state one at t 100t 0 .
The current distributions over a dielectric bulk at different moments are shown in Fig. 4. Figure 5 shows the flowchart for stationary case simulations. The current distributions for non-stationary cases were calculated by solving Eqs. (5)-(7) within standard integration techniques. All simulation parameters are shown in Table 2. One can see that the drift current is mostly dominant in the bulk, while the diffusion one reaches significant values near contacts. This means that the diffusion contribution can not be neglected, especially, in case of strongly out-of-balance (transient) states that take place in real processes in modern ultra-fast microelectronic devices 32 . As expected, the total current does not depend on the coordinate in the stationary case.
It should be noted, that it is not possible to measure j drift and j diff separately, only j total can be measured in real experiments. However, it is possible to calculate j drift and j diff . The algorithm is following: 1. measure J-F-T characteristics; 2. fit simulated curves with experimental ones using any microscopic models of injection from contacts and trap capture/emission processes, obtain the trap parameters (trap ionisation energy, trap density, attemptto-escape rate etc.); 3. calculate n distributions over the dielectric bulk at different voltages using Eq. (9) and obtained trap parameters; 4. calculate j drift and j diff distributions over the dielectric bulk at different voltages using Eq. (10) and obtained trap parameters.
Equation (5) has an obvious generalization for the case:  www.nature.com/scientificreports/ Here i is the coordinate index. Equation (11) can be written in a vector form as follows: where P ± tun = (Tun ± (F x ), Tun ± (F y ), Tun ± (F z )) is a vector of tunnelling probability, the symbol • defines the Hadamard product (also known as the Schur product or the entrywise product): where A and B are two matrices (vectors) of the same dimension.

Model of phonon-coupled traps in dielectrics. Phonon-assisted tunnelling of electrons (holes) between
neighbouring traps. Recently, it has been shown, that the PATENT model adequately describes the current in high-κ dielectrics 30,[33][34][35] . An energy diagram of the electron tunnelling from a phonon-coupled trap to the other one at a distance of a in an external electric field is shown in Fig. 6. The energy dependency from configuration coordinate Q of a system trapped-electron-plus-phonons is shown by U b (Q) curve. The U f (Q) curve corresponds to a "free" electron in the conduction band. Solid lines show the initial state (before the tunnelling act), dashed lines represent the final state (after the tunnelling act). Due to the effect of the external electric field, electrons localized on the neighbouring traps have different energies represented by slanted lines ε(Q) in Fig. 6, and the tunnel event must be accompanied by inelastic processes, like phonon emission and/or absorption to compensate the energy difference. All the above is taken into account in the PATENT model.
According to the PATENT model, the tunnel junction rate between neighbouring traps is the following: where ℏ is the Dirac constant, m * is the effective mass of trapped electron (hole), k is the Boltzmann constant, T is temperature, q is the elementary charge (absolute value), W t and W opt are thermal and optical trap energy, respectively, and Q 0 is the configuration coordinate that characterize the electron-phonon interaction. In case of low electric fields qFa ≪ W t , Eq. (12) can be simplified to  www.nature.com/scientificreports/ where, the pre-exponent ratio is the attempt-to-escape rate, the first exponent is the thermal ionization effect (with the activation energy of (W opt − W t )/2 ), the second exponent is the tunnelling factor and the hyperbolic sine represent the activation energy decreasing due to the external electric field with respect to transitions in both co-field and contrafield directions. The analysis shows that, in non-low electric fields ( qFa 1/4W t ), Eq. (13) gives a large deviation from experiment data, and Eq. (12) should be applied instead.
Ionization of phonon-coupled traps to contacts. Ionization of phonon-coupled traps to contacts is described according to Refs. 29,36 by the following equation: Here is difference between the energy of the dielectric conduction band bottom and the Fermi level in the contact, a * is the trap-to-contact distance, V out the free electron velocity in a contact, f F-D (Ẽ) is the Fermi-Dirac distribution function, Ẽ is the electron energy relative to the Fermi level.
Additionally, thermal ionization, which takes place in case of zero electric field F = 0 , should be taken into account: The total ionization rate is composed of terms (14) and (15): Injection from contacts to phonon-coupled traps. The injection rate Inj from the contact to the trap is similar to Eqs. (14)- (16), taking into account that the electron is injected from the occupied state: (16) Ion = Ion F , for F � = 0, Ion T , for F = 0. The measured current-field characteristics J-F-T are shown by symbols in Fig. 7. The calculated J-F-T characteristics using equations (9), (10), (12), (14)- (19) are represented by lines in Fig. 7. One can see that the calculations in terms of a developed model are in a good quantitative agreement with the experimental data. In the calculations, the electron barriers on contacts as adopted parameters were used as breaks in the conduction band bottom in the contacts and dielectric � Si/HfO 2 and � HfO 2 /Ni that are in agreement with the literature data 35,37 , thermal and optic trap energies are corresponds to an oxygen vacancy in HfO 2 30,35 , the mean distance between traps a = a * = 1.8 nm corresponds to the total trap density N = 1.7 × 10 20 cm −3 that is close to value 2.5 × 10 20 cm −3 obtained in the original article 30 , while the effective electron mass m * = 0.59m 0 (here m 0 is the free electron mass) is closer to the values calculated within ab initio simulations 38 than the original value m * = 0.8m 0 that was obtained when ignoring the diffusion contribution 30 . All simulation parameter values are shown in Table 3.

Discussion
We have developed a continuum model of the charge transport in dielectrics, which takes into account the drift and diffusion contributions to the trapped charge carriers current. The developed continuum model is universal, i.e. it can be supplemented with a different injection from contacts and trap capture/emission microscopic processes, including Hill-Adachi, PATENT, TAT and other trap ionisation models. It is shown that neither of contributions can be neglected in case of ultrathin dielectric films. The analytical stationary solutions, that are close to experimentally achievable conditions, are found. A comparison with other theoretical models and experiments shows a good agreement. The model can be extended to two-and three-dimensional systems. The developed model, formulated in partial differential equations, can be implemented in the finite element method code that is compatible with other partial differential equations, e.g. Poison and thermal equations. The joint solution of these equations can significantly advance the search of optimal parameters for electronic devices and conditions for their fabrication, including promising ones: RRAM, FRAM and many others. The found analytical solutions will be useful for the verification of various numerical codes for simulations of physics of electronic devices, including stress induced leakage currents in FETs and flash memories, currents during the resistive switching of RRAM and memristors and leakage currents in FRAM devices.  www.nature.com/scientificreports/

Methods
Preparation of samples for transport measurements. Transport measurements were performed for metal-insulator-semiconductor (MIS). For the MIS Si/HfO 2 /Ni structures, the 20-nm-thick amorphous hafnia was deposited on a n-type Si (100) wafer by using the atomic layer deposition (ALD) system as described according to Ref. 30 . Tetrakis dimethyl amino hafnium (TDMAHf) and water vapour were used as precursors at a chamber temperature of 250 • C for HfO 2 film deposition. The samples for transport measurements were equipped with round 50-nm-thick Ni gates with a radius of 70 µm.
Transport measurements. Transport measurements were performed using a Hewlett Packard 4155B semiconductor parameter analyzer and an Agilent E4980A precision LCR meter.