Abstract
We propose a model of the kinetics of reversible breakdown in metal-insulator-metal structures with afnia based on the growth of fractal patterns of defects when the insulator is subject to an external voltage. The probability that a defect is (or is not) generated and the position where it is generated depend on the electric field distribution. The new defect moves accordingly to fractal rules and attach to another defect in a tree branch. When the two electrodes sandwiching the insulating film are connected a conductive filament is formed and the breakdown takes place. The model is calibrated with experiments inducing metastable soft breakdown events in Pt/HfO2/Pt capacitors.
1. Introduction
The loss of insulating properties of thin oxide films is due to the phenomenon of dielectric breakdown. The irreversible breakdown is one of major sources of failure in integrated circuits based on metal oxide semiconductor field effect transistors (MOSFET) technology. Reversible soft breakdown of gate dielectrics is not fatal with respect to the MOSFET switching but leads to undesired current leakage through the gate which affects the power performance. On the contrary, deep control of the reversible changes of conductivity between high and low resistance states in high-k based metal-insulator-metal (MIM) capacitors is today strongly desired for application in resistive switching memory (RRAM). In fact, among the technology options for nonvolatile memory devices, RRAM is gaining an increasing consideration for its superior scaling opportunities, high speed, and low operating voltages [1–3]. For all these reasons, understanding reversible breakdown in high-k dielectric films is today of major interest for microelectronics. The metastable soft breakdown of the dielectric film is due to the built-up and rupture of conductive filaments (CF) formed of defects which can be achieved by consecutively applying proper voltages at the two electrodes. While the key factors driving the resistive switching have been identified and already widely modeled [4–11], the kinetics of CF formation is still under debate. At an atomistic scale, some insight in the time evolution of filaments is given by Monte Carlo modeling in 2D [12] and 3D [13], which relates the time to breakdown to the applied voltage and the filament evolution to vacancy migration. The model proposed here describes the kinetics of growth of conductive filaments until breakdown within an external driving factor (voltage pulses ) in terms of fractal aggregation of defects.
2. Method
The damage structures observed in dielectric breakdown in many cases have the form of trees [14–18]. A tree is at least at some level of approximation a fractal [15] or self-similar object in the sense that a branch of the tree again looks like a tree and so on. In fractal structures the relation of the total length of all branches inside a circle of radius (the number of all the grid points) and the radius itself is a power law with noninteger exponent FD: , where FD is the fractal dimension. 2D numerical modelling of dielectric breakdown in terms of fractal patterns of defects was already proposed in the past, in the case of solid [19], liquid [20], and gaseous [21] materials. In all those simulations and visual inspections, breakdown in dielectric media was always found to feature fractal patterns with . In the present case the dielectric breakdown in metal-insulator-metal (MIM) capacitors is studied. The defect distribution is filamentary, growing from one electrode to the other. The condition will be respected. The filaments are assumed to be electrically conducting so that the field at the filament tip exceeds the original field at the electrode plate. This leads to an amplification, propagation, and multiplication (branching) of the local field enhancement. So, once a characteristic value of the electric field is exceeded, a rapid and abrupt increase of the charge density near the filament tip is achieved. The avalanche behaviour is a key feature of the phenomenon descending from its fractal or self-similar nature. With respect to previous study of breakdown of dielectrics in MOS capacitors, in the present model there are two fundamental improvements: (i) the probability that a defect is (or is not) generated depends on the maximum value of electric field in the mesh and (ii) the position in the grid where the defect is generated depends on the local electric field. The cross-section of the oxide film is meshed by a 2D grid. The matrix is given by elemental cells with size (representing a point). At a generic time , a distribution of defects is present, as depicted in Figure 1 (with , ): black cells of the matrix represent point defects. Breakdown occurs when defects are piled up in a way to connect the two electrodes. The distance to be covered to complete a filament is (dashed arrows in Figure 1) and the minimum distance in the whole matrix is (bold red arrow). In many situations the local electric field in each row of the matrix is higher than due to the needle shape of neighbor filaments and this accelerates the filament build-up. A rigorous evaluation of the field distribution would require burdensome numerical simulations. In our treatment, it will be posed and the slowest case of breakdown will be considered. This leads to possible underestimation of the times to breakdown. The steps of the modelling procedure are resumed in the flow diagram shown in Figure 2.
At time zero a pseudo-random 2D distribution of native defects in the dielectric is set, resulting in a mix of point defects and portions of filaments (chains of more than one point defect). Such initial distribution recalls the presence of preferential defect lines along grain boundaries already existing in the polycrystalline material as deposited. This is an approximation of the real distribution of native defects which, obviously, extends over the three dimensions. On the other hand, the 2D model will be calibrated with experimental results performed on sets of real samples. Therefore, the variability of the 2D native defect distribution coming out from the final comparison with experiments will actually be overestimated since it is in charge also of the variability in the third dimension. obeys to two conditions. The first condition regards the maximum length of native filaments; the second one regards the maximum area inside which native filaments may exist: (1) the maximum length of native portions of filaments is a fraction of the oxide thickness; (2) the maximum area occupied by native defects is a fraction of the total matrix area. As physical intuition suggests, the distribution of native defects affects the statistics of breakdown. In Figure 1, there was displayed example of , with and . A sensitivity study to parameters and is shown in Figure 3. To the aim, some different distributions were taken into account varying values of in the range between 10% and 50% and in the range between 5% and 15% (, ) at the same voltage ( V). The number of iterations to breakdown is calculated for a number of virtual samples featuring those values of and and results are displayed in Weibull plots. In Figure 3 some of the pseudo-random distributions better attain to the Weibull statistics, and some others do not, but this is not the point. In fact, the native defect distribution which will be finally adopted in our simulations will be obtained by fitting the Weibull distributions of experimental data of time-to-breakdown measured in a wide range of electrical conditions. As one can see, the calculated Weibull slope () increases with increasing (the maximum number of native defects) and decreasing (the maximum length of native filaments).
(a)
(b)
The electric field distribution in the matrix is calculated at the step 2 of Figure 2. At time , a new defect can (or cannot) be generated, with a probability which depends on (step 3 of the flow diagram). The analytical expression of the distribution probability is a trunked Boltzmann-like equation where is the local electric field to breakdown and is a characteristic parameter. The truncation defines the validity range of the model, since when the electric field exceeds the critical value the breakdown occurs certainly. The Boltzmann-like model of probability expressed in (1) assumes a physical meaning when it is compared with the generation probability holding in the thermochemical model (TCM) proposed to explain the field-dependent breakdown in thin SiO2 films [22]. This point will be developed in the next section. When only few defects are present (e.g., at short times), distance is great and is low with respect to . The ratio generations/iterations are very sensitive to the matrix size and applied voltage and may vary in several orders of magnitude in practical conditions. When the new defect is generated the procedure goes to the next step; otherwise another iteration is performed (step 4 in Figure 2). Generation only at metal interfaces is considered since, in general, defectiveness at interfaces is increased with respect to bulk [23–26]. The position of the newly generated defect (i.e., the column of the matrix) follows a discrete probability distribution () which is related to the electric field in that column through the relation [21, 27]: The sum in the denominator is a normalized factor as to have . This corresponds to step 5 in Figure 2. at and at the generic time is sketched in Figure 1 as a function of the matrix column (once the defect distribution is given). In the 2D approach, at each step, the generated defect can move in five possible directions on the cross-section plane (right and left diagonal, right and left lateral, and downward), occupying one of the five adjacent cells. Finally, it aggregates to another existing defect in a tree branch (or reaches the other side) and stops (step 6 in Figure 2). The breakdown fractal structures may be reproduced by the Diffusion Limited Aggregation (DLA) model [16, 17], introduced as a model for irreversible colloidal aggregation. In the original DLA algorithm at only a seed is present in a 2D plane. The first particle is introduced from infinite distance and moves randomly in the plane until it encounters the seed and stops. Subsequently, other particles are added with no correlation, one at time, from infinite distance. They move randomly in the plane until they aggregate to the preexisting structure.
In the case of a structure with parallel plane faces, the fractal tree starts from one of the two interfaces, and then the DLA algorithm cannot be used as it is. We modified it in order to generate fractals structures describing breakdown in metal-oxide-semiconductor or metal-insulator-metal structures. In contrast with the DLA, in the present case particle displacements are not random, but probabilities are ruled by quantities and (whose values range from 0 to 1), which definitely determine the fractal dimension of the tree. Particle movement probabilities are defined as follows: In particular, if the particle lateral movement is inhibited and the straight downward probability (and the fractal dimension) increases with . On the contrary, if the straight downward movement is inhibited and the particle can only move laterally or obliquely. In this case, the fractal structure is as spread (and the fractal dimension decreases) as is higher. For this reason, we call lateral probability and downward probability. The choice of the values of the two parameters and is arbitrary, within the frame of a fractal tree featuring , in accordance with optical inspections reported in literature [16, 17] and literature findings [19–21] for breakdown phenomena in dielectric media, and under the condition , in accordance with physical intuition for the problem under investigation. In our simulations was definitely ensured setting and . Defects attaching to a tree branch concur to the formation of a conductive filament. When a branch of the tree touches the electrode breakdown occurs and the simulation ends (step 7 in the flow diagram). The time to soft breakdown () is the time necessary for the filament forming, which is proportional to the number of iterations to breakdown. It is worth noticing that the tree grows in a volume of the bulk oxide (3D), not on a cross-section plane (2D), and therefore the particle actually has more than the five possible directions described in (3). However, considering the third dimension would not require the introduction of additional equations nor other parameters distinct from and , since the probabilities of the particle movement in the new directions (along a plane parallel to the electrode plates) are the same. Therefore studying the tree growth only in cross-section planes is not limiting. At the same time, simulations performed in 2D are much faster and, definitely, in this work they will be preferred to the 3D approach.
The proposed methodology can be of interest for the study of the electroforming in the resistive RAMs (RRAMs) based on MIM stacks. Electroforming consists in the operation of soft breakdown in fresh devices due to defects (oxygen vacancies) piling up between the two electrodes when the insulator is subject to electrical stress. In the last years, the phenomenon of electroforming in RRAMs has been widely studied experimentally and simulated (see, e.g., [28–31]), since those devices are of outmost interest in microelectronics being candidate to replace floating-gate nonvolatile memories in the future nodes [32]. For this reason, in the next section the proposed model will be applied on MIM capacitors designed for RRAM applications.
3. Results and Discussion
The model has been calibrated using experimental results obtained on MIM capacitors of the type Pt/HfO2/Pt. The HfO2 film (10 nm thick) was ALD deposited on sputtered bottom electrode. Also the top electrode was sputtered and the stack patterned via hard mask etching. Other details on sample preparation, material analysis, and complete characterization of samples in quasi-static condition are reported elsewhere [4, 33]. The values and fix the size of the simulated device being the cell size nm (consistent with defects originated by oxygen vacancies [28]). In order to measure , transient measurements were performed at extremely short times. To the aim, voltage pulses with 10 ns rise time and amplitude . Pulses were displayed on channel 1 (CH1) of an oscilloscope. The voltage drop on a resistor in series to the device under test () was acquired on channel 2 (CH2) of the oscilloscope. kΩ fixed the current compliance to kΩ. The signal displayed on CH2 allowed measuring . As a result, the measured values of are displayed in Figure 4(a) (symbols) as a function of .
(a)
(b)
It is worth noticing that the voltage behavior of is similar to literature data regarding the same insulator and different electrodes [5], but also different insulator, electrodes, and thickness [31]. Our data were compared with simulations performed with the procedure described above. An excellent fit of the experimental trend of versus was achieved using cm/V and MV/cm in the expression of . It is worth noticing that is the value of the local electric field on the edge of the filament which can be much greater than the value of the so-called breakdown electric field calculated by dividing the breakdown voltage by the oxide film thickness (typically lower than 5 MV/cm in HfO2). Once the calculated trend reproduces the experimental one, there is another step regarding times to be completed. In fact, the output of our modeling procedure does not contain any information about absolute times, since the discretization is given by iterations. Therefore, in order to interpolate measured data, a trivial normalization is made, forcing the calculated point at V (obtained after a certain number of iterations) to coincide with the point measured at the same voltage ( s). In this way, a time interval is associated with each iteration: and the variable of discretization becomes the time to breakdown. The calculated curve of with cm/V and MV/cm is drawn in Figure 4(a) with continuous line. For comparison, other two curves with different values of are also drawn in the same figure, which outlines the sensitivity of the procedure to parameter . The fitting curve displayed in Figure 4(a) was calculated using the particular distribution of native defects which obeys to the condition that the experimental Weibull distributions are fitted at any experimental conditions. The Weibull plot in Figure 4(b) includes experimental data got at several different voltages (symbols) and the calculated curves obtained with featuring , (continuous lines).
Recalling the thermochemical model proposed to explain the field-dependent breakdown in thin SiO2 films [22], we have the following expression of the generation rate: where is a rate constant, which depends on intrinsic material properties, is the bond polarization factor, is the effective activation energy of the defect generation. The Boltzmann-like generation probability defined in (1) (for ) and the generation rate defined in (4) can be matched. Then, it comes out that and , with being the local temperature where the defect is generated. Starting from this correspondence, trivial algebra gives proportionality between the local temperature where the defect is created and the activation energy for the creation: Using cm/V and MV/cm, it comes out K·eV−1, which means that with activation energies of a few eVs the local temperature for defect generation is several hundreds of kelvin degrees.
The simulation procedure yields the defect distribution at each iteration. Its evolution until soft breakdown is shown in the 3D plot shown in Figure 5 for the case V. In that plot the time axis is logarithmic, starting at one microsecond since no variation of the defect distribution is appreciable before that time. It is worth noticing that defect filaments begin to grow suddenly after a few microseconds and the first conductive filament is completely formed within a few tens of microseconds. This is the case of an avalanche kinetics rather than a progressive phenomenon. The latter result has a general validity and it is not influenced by the 2D approximation, which probably affects the timing and the entity of the phenomenon.
Further considerations may be done about the rate of generated defects. The (minimum) number of defects () generated in the 2D approach is now plotted in Figure 6(a) as function of time, for different operation voltages. The trend in Figure 6(a) is in agreement with [29]. In that work the number of defects only in the breakdown path is reported and therefore it is much lower than the minimum total number of defects at breakdown reported in Figure 6(a).
(a)
(b)
Those curves may be derived numerically and the calculated rates are reported in Figure 6(b) with symbols. The voltage dependence of the defect rates displayed in Figure 6(b) is not descending from the dimensionality of the approach and therefore is holding in general. Then, the calculated rates are posed equally to the generation rates expressed in (4): . depends exponentially on the electric field, which can be calculated using the fractal growth of the conductive filament. An analytical expression of the electric field as a function of time would be desired, to treat analytically the differential equation . An analytical expression is obtained interpolating the calculated maximum values of electric field. The calculated values of electric field and the analytical interpolating functions are drawn in Figure 7 for several operation voltages, with symbols and continuous lines, respectively. The analytical curves have the following expressions, where , and are constants: We wish to remark here that (6) is not a physically based model but an analytical interpolation of data. Substituting (6) into (4) one obtains In Figure 6(b) the continuous lines represent interpolations of the calculated values obtained with the function reported on the right hand in (7) using the constants values found previously. The procedure allows evaluation of the behavior of the parameter . As a result, is not appreciably dependent on the applied voltage. Quantitatively, we found an average value of around 6 × 108 defects/s, which means that generation of one defect takes approximately 1.6 ns at the maximum electric field. This value should be interpreted as an indication of the timescale of the phenomenon, since it is affected by the 2D approximation. Finally, the voltage and thickness dependence of the number of defects at the moment of soft breakdown () is drawn in Figure 8. As a result, the calculated does not depend on voltage and increases linearly with the oxide thickness. The linearity of can be fairly extended to 3D, while the absolute values of defects at the breakdown should rather be interpreted as an order of magnitude, fruit of the projection on a cross-section plane of 3D trees and branches.
4. Conclusions
It is proposed a model of the formation kinetics of conductive filaments in oxides subject to an external driving factor (voltage pulses). The model is based on fractal aggregation of defects and is inspired by literature works extensively reporting the fractal nature of breakdown events in solids, liquid, and gases, where a fractal dimension around 1.8 was reported. It proposes a modified version of the diffusion limited aggregation model of breakdown fractal trees, in which defect displacements are not random, but probabilities are ruled by quantities which definitely fix to 1.8 the fractal dimension of the tree. It takes into account the presence of native defects and a defect generation probability depending on the local electric field. With respect to other models, the proposed approach predicts the time rate of defect formation, the avalanche onset of breakdown, and the dependence of the average number of defects at breakdown on the insulator thickness and the applied voltage. The proposed methodology can be of interest for the study of the electroforming in the resistive RAMs (RRAMs) based on MIM stacks, which are of outmost interest in microelectronics being candidate to replace floating-gate nonvolatile memories in the future nodes. The model is calibrated comparing simulation results with experiments performed on Pt/HfO2/Pt MIM cells designed for RRAM applications. Excellent simulations of the time to soft breakdown versus voltage curves and of the Weibull distributions are obtained meshing the structure. A 2D approximation has been assumed for brevity, since it allows getting unexplored insights on the breakdown phenomenon with reasonable simulation times. The study of the time to breakdown reveals a sort of avalanche mechanism, with a sudden onset and a rapid conclusion, as predicted in a fractal approach. The correspondence between the physically based model of thermochemical breakdown in dielectrics and the present one yields a proportionality between the local temperature where the defect is created and the activation energy for the creation of the new defect. The average number of defects building-up in a filament can be predicted as function of time and is not dependent on voltage. The timescale for generation of one defect is estimated in the order of a nanosecond. Similarly, the number of defects to breakdown can be predicted as function of voltage and film thickness. As a result, the number of defects in a conductive filament does not depend on voltage and increases linearly with the film thickness.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.