Reversing A Decades‐Long Scaling Law of Dielectric Breakdown using Hydrogen‐Plasma‐Treated HfO2 ReRAM Devices

Dielectric breakdown (BD) is known to cause component failure in electronic devices and high‐voltage power lines over many decades. In recent years, this failure mechanism has been exploited to intentionally form nanoscale filaments in resistive random‐access‐memory (ReRAM) devices for artificial intelligence (AI). The statistical nature of this failure mechanism, known as the inverse size scaling law based on the weakest‐link theory, dictates the ever‐higher forming voltages for ReRAM devices, in turn, requiring additional peripheral transistors. This law also causes diminishing lifetimes of electronic components for future generations of very‐large‐integrated circuits (VLSI) technology. Currently, the semiconductor industry faces these scaling barriers limiting the miniaturization of both VLSI and AI hardware technology. Here, experimental evidence of a reverse area dependence is presented by introducing an innovative fabrication process with a hydrogen‐plasma‐treated layer in a bilayer HfO2 structure. Using joint order statistics rather than traditional extreme‐value statistics, a physics‐based statistical model is developed in agreement with the experimental data, thus demonstrating that there is no fundamental reason preventing this law from being reversed or altered. These findings will have a significant impact on both technology scaling and fabrication innovation for electronic and/or bioinspired nanomaterials; moreover, stimulate much research in physics, statistics, and reliability.


Introduction
Since dielectric BD phenomena such as fracture (mechanical) or lightning (electrical) have become known to humankind, DOI: 10.1002/aelm.202300296 the size effect of dielectric BD [1][2][3][4][5][6][7][8][9] was already recognized as one of two concerns of "new sciences" by Galileo Galilei in 1638. [10]Induced either by a stress load or by an electric field, the size effect has manifested itself as a geometrical (either area or volume) dependence of dielectric failure.In the electrical BD of dielectrics, this failure is quantified by the measurements called time-to-BD (T BD ) or voltage (field)-to-BD (V BD or F BD , also called dielectric strength).In reliability assessments of polymers used in highvoltage power cables for long-distance transmission, this geometric dependence is known as the enlargement law. [6]In modern VLSI nanoscale semiconductor devices, ultrathin dielectrics are used in a number of different applications.In the front-end-of-line (FEOL), they insulate gate electrodes from the conducting channel.In middle-of-line (MOL) and back-end-of-line (BEOL), they separate various metal wires from each other and serve as the dielectric of capacitors.As such, a single BD event can lead to the failure of an entire chip. [7,8]While dielectric BD is usually regarded as a failure process due to loss of the insulating property in most applications, in recent years controlled dielectric BD has been used to create a conducting filament in ReRAM devices. [8,9,11]For this device type, the term "forming" is frequently used instead of "breakdown", leading to analogous parameters called forming voltage (V FORM ) and time-toform (T FORM ).Similarly, dielectric BD can also be used to fabricate single in-situ solid-state nanopores for DNA sequencing although avoidance of the occurrence of multiple pores in largerarea membranes can be challenging. [12]Regardless of whether dielectric BD is considered a failure mechanism or a fabrication tool, all the experimental data from BD or forming voltage measurements in these different applications reveal this wellestablished inverse area scaling law [4][5][6][7][8][9]11] with a negative slope as shown in Figure 1a which uses data from ReRAM devices as an example.
This scaling law now presents a bottleneck limiting technology miniaturization for ultradense microelectronic circuits, thus preventing the realization of Moore's law for the growth of semiconductor technology.Because of time limitations, measurements of T BD_FORM or V BD_FORM can only be performed on Figure 1.Area-dependence of T BD_FORM or V BD_FORM and their impact on the projection of forming voltage and lifetimes at the specified yield or reliability targets.a) Experimental V BD_FORM or T BD_FORM data for conventional area-dependence of HfO 2 (Sample D for 6 nm and Sample E for 5 nm) obtained by RVS and CVS techniques with a ramp rate of 0.4 V s −1 , respectively.For comparison, the reverse area-dependence data is also included using bilayer HfO 2 stack (Sample A).The detailed discussion of sample A is given later in the text.b) Projection procedures according to conventional area-dependence of T BD_FORM or V BD_FORM for forming-voltage of functional ReRAM devices and the lifetime of an electronic chip from the test device areas, respectively.Note this involves only a horizontal projection for area extrapolation.In addition, a projection from the limited range of the distributions of test devices with limited samples to an extremely small percentage specification (≈10 ppm = 0.001%) is required for an electronic chip 7 .In contrast, for ReRAM devices, a projection to 99.9999-100% is required to ensure high yield in production.Therefore, the knowledge of the functional form of the CDF becomes essential in the second projection step.Note that for a full reliability assessment, a projection of stress BD data from the stress voltage to the use voltage is also required. [7]This step is omitted in (b) for clarity since area-dependent BD is the focus of this work.
devices of certain areas, and on limited numbers of samples.Thus, an extrapolation from the test device area to the total chip area at a cumulative failure specification (an extremely low value of 10 ppm, parts per million) is required as shown in Figure 1b.In the case of reliability assessment, the lifetimes of a chip with a much larger area can be several decades shorter than those of the testing devices, failing to meet product reliability specifications. [7]t would be beneficial to eliminate this area dependence so that devices with different areas can meet the same reliability specification.In contrast, to successfully form ReRAM devices with a yield of 99.9999%, a much higher V FORM is required.The voltage for this one-time forming process is much higher than normal switching voltages, [8,9,11] thus posing a design showstopper because peripheral circuitry needs to include several transistors for each ReRAM device to support this high forming voltage, limiting chip-size downscaling.Therefore, a reversal of this areadependence law would be technologically advantageous.
[3][4][5][6][7][8][9] Defects are generated in an insulator when a voltage is applied across it.These defects begin to percolate and string together to form filaments between the electrodes; the filament which takes the shortest time to span the distance wins the competition of filament formation.Nearly a century ago (in 1926), Peirce [1] was probably the first to realize the connection between dielectric strength and the weakest link statistics theory.Specifically, the weakest-link theory, also called the minimum value theory, [13,14] can be expressed as follows where F u (t) is the cumulative density function (CDF) of BD or filament formation for a unit element of a larger sample.F(t) is the CDF of the entire sample with N number of individual elements.Its corresponding probability density function (PDF), f(t) or f u (t), is given in Equation (2).F u (t) can be considered from either the statistical formulation [1][2][3][4][6][7][8]11] or a physics-based model [5,9,15] for the description of experimental data. Mathemaically, the well-known Weibull distribution, F u (t BD_Form ) = 1−exp−(t BD_Form / BD )  with  BD and  being the respective characteristic scale and shape parameters, is known to preserve its functional form [13,14] after applying Equation (1).If a large area (A) device is composed of N number of equal smaller area (a 0 ), that is A = Na 0 , the application of Equation (1) using the Weibull distribution yields the area dependence of T BD_FORM in an inverse power-law form as follows [7,8] The same power-law form of area-dependence as Equation (3) can be derived [11] for the measurement of V BD _ FORM as shown in Figure 1a whereas an exponential law was also used to fit the BD data. [9]The experimental data supporting this law was reported in the applications of microelectronics technology as earlier as 1979. [4]ver the last several decades, several physics-based models were advanced to explain the statistical nature of dielectric BD based on either numerical simulation or analytic solution of percolation models. [5,9,15,16]One of the pioneering works is to introduce the time dependence of defect generation into their percolation process in a geometrical model. [5]A common feature of these models is based on the concept of a critical defect number, N C , which provides a criterion to trigger a BD event as a filament is formed across an insulator.This cell-based percolation model [5] was the first to correctly capture the experimental T BD_FORM area-dependence ascribed by Equation (3).To our best knowledge, most, if not all, experimental data on dielectric BD reported in the literature; including observations in mechanically induced material failure phenomena [2,3,17] also follow similar areal or volumetric dependence with a negative slope ascribed by the minimum value theory.
Although both experimental observations and physics-based statistical models have established this inverse area-scaling law of dielectric BD over many decades, there is no fundamental physical basis to prevent this law from being changed or reversed.Moreover, the strict application of the minimum value theory of Equation ( 1) is merely a special case of the general framework of order statistics theory. [14]In this work, we will show this scaling law can be reversed by introducing a hydrogen plasma process in the fabrication of our ReRAM devices.An example of the experimental evidence of the reverse area-scaling is shown in Figure 1a with more results to be presented in Figure 4.Moreover, we develop a physics-based statistical model to demonstrate that a defect annihilation process in competition with a defect generation process (Figure 3a) is at the core to lead to such a reversal process while filament forming due to dielectric BD can still take place with satisfactory ReRAM switching performance.

Results and Discussion
The HfO 2 ReRAM devices used in this work are fabricated using the process flow described in Figure 2a. Figure 2b provides a description of device structures of various layers used in the fabrication whereas a transmission electron microscopy (TEM) image of the active device area is shown in Figure 2c.Additional process information is summarized in the experimental section.
The key aspect of the process steps involves a hydrogen-based plasma process [18][19][20] which is intended to introduce oxygen vacancies and interstitials in HfO 2 stacks.Figure 2d provides an overall comparison of forming voltage data collected for a wide range of hydrogen exposure times and different combinations of bilayer thickness values.As seen in Figure 2d, a continuous decrease in forming voltages with increasing hydrogen exposure times.In particular, samples A, B, and C show the lowest forming voltages among the HfO 2 stacks we investigated.More importantly, these samples also demonstrate the reverse area scaling of dielectric BD.To our best knowledge, this unique reverse area scaling characteristic has not been reported in the past.
Our physics-based statistical model is based on general order statistics theory [14] by incorporating both defect generation and annihilation processes as schematically described in Figure 3a.The geometrical cell-based structure used in our model is the same as previously considered for conventional area-dependent BD. [5] Unlike a single random sequence adopted in the minimum value theory, [5,13] we consider the joint probability density function of two independent random sequences, [21] t 1 and t 2 , respectively as given below: Specifically in our model, t 1 and t 2 represent defect generation time, t G , and annihilation time, t A , respectively.More importantly, f(t 1 ) and f(t 2 ) should be based on order statistics theory [14] for the description of dielectric BD.To ensure the weakest link with other splits which are reported previously in Ref. [18].We included these results here for a broad view of the general trend and impact of H incorporation in forming voltage reduction.9).They are used for the model calculation to compare with the experimental data at different voltages in Figure 4a,e.d) The CDFs, F u (t G ) and F u (t A ) of Equations ( 11) and ( 12) for the cell area (a 0 ) using the respective  G (t) and  A (t), with N G = 6.2 × 10 17 cm −3 and N A = 42 × 10 17 cm −3 at 1 s with  G = 0.08 and  A = 0.15, respectively.The corresponding PDFs, f u (t G ) and f u (t A ), are also included for comparison.The value of N C is 8 whereas t diel is 5nm with a 0 = 25 nm 2 as listed in Table 1.The grey area represents the time-window over which the experimental data collected from larger device areas.e) The calculated joint PDF distribution of Equation ( 7) for three different areas using the cell distributions (d).f) The marginal PDF distribution of these three areas obtained by integrating Equation (7).The calculation results shown in (c),(f) eventually lead to the resultant CDF of Equation ( 13) to compare with the experimental data in Figure 4a.Note that the dominant role of the defect annihilation process in conjunction with the defect generation process in (c),(d) in producing reverse area-dependence of T BD_FORM .
property of dielectric BD, the PDF, f(t 1 ) or f(t G ), of the entire sample corresponds to Equation (2) for the shortest time leading to filament formation.For the defect annihilation process, we consider the PDF, f(t 2 ) or f(t A ), of the entire sample follows the highest order statistics to maximize the impact of the annihilation process.According to general order statistics, [13,14] the CDF and PDF distribution of the highest orders, or the maximum values, of a random sample given below With Equation ( 2) and Equation (6), Equation (4) takes on the form as follows In addition, a condition of t G ≤ t A is introduced to ensure that defect generation precedes the annihilation process.It is worthwhile to note that the formulation of Equation ( 7) is entirely based on order statistics [13,14] which is generally valid and applicable for dielectric BD and it is also independent of choice of f u (t G ) and f u (t A ).To obtain the PDFs or CDFs of a unit cell, we assume both generation and annihilation processes follow a Poisson process where N C is the critical number of defects required to trigger BD or filament formation. is the average number of defects, and it is related to the average defect density, , as  = t diel a 0 , where t diel is the dielectric thickness and a 0 is the cell area as shown in Figure 3b.For time dependence of defect generation, we adopt a power-law dependence, G as previously proposed in modeling of gate dielectric BD. [15] Note that we use the subscript, G, specifically for the defect generation process.As a result, the time dependence of the average defect number can be summarized below For the defect annihilation process, we also consider the Poisson process of Equation ( 8), but  corresponds to the average number of annihilation centers so that we specify ≡ A to distinguish it from the average number of defects,  G .Similarly,  A is related to the average density of annihilation centers,  A , as  A =  A t diel a 0 .The critical number N C in Equation ( 8) should represent the critical number of annihilation centers involved in BD or forming process, as its physical meaning will be discussed later.In the case of the annihilation process, we propose a power law for the time dependence of its average density but with a decaying function as  A = N A t − A A .Thus, the time-dependence of the average number of annihilation centers becomes as follows where N G and N A are the density of generated defects and annihilation centers, respectively. G and  A are the respective generation and annihilation rates.Then, the respective CDFs of defect generation and annihilation process can be expressed as follows Their corresponding PDFs, f u (t G ) and f u (t A ), can be readily derived by differentiation of Equations ( 11) and (12).Finally, the CDF of this T BD_FORM involving both defect generation and annihilation processes can be obtained from the marginal PDF of Equation ( 7) by integration over defect annihilation times as follows where A is the normalization constant.In the following, we will present our model calculation in comparison with the experimental data of ReRAM HfO 2 devices as an example for the model validity demonstration.The time-dependent density of generated defects and annihilation centers are shown in Figure 3c for some selected parameters.As expected, defect density continuously increases with time due to the generation process responsible for filament formation.It is generally accepted that oxygen vacancy and interstitial oxygen are responsible for HfO 2 ReRAM switching operations.Thus, in our model (Figure 3b), these vacancies are identified as defect traps for carrier conduction where interstitial oxygens would carry on the role of annihilation centers as discussed so far.The CDFs, F u (t G ) and F u (t A ), of Equations ( 11) and ( 12) for the cell area of a 0 are shown in Figure 3d whereas their PDFs are also shown for comparison.This case illustrates the predominant role of the annihilation process with respect to the defect generation process.A full calculation of the joint PDF distributions is presented in Figure 3e, revealing a shift to longer times with increasing areas.Because both generated defects and annihilation centers increase with increasing area, the presence of annihilation centers with a higher concentration of annihilation centers (Figure 3c) diminishes the impact of defect generation for larger area devices while enhancing filament formation for smaller area devices.This gives rise to a reduced marginal PDF with increasing areas (Figure 3f), ultimately leading to a reverse area-dependence of CDF, consistent with experimental observations such as those in Figure 4a,e as we will discuss further below.
We are not aware of any publications adopting the decay function of Equation (10) for the defect annihilation process as proposed here in dealing with the formulation of dielectric BD or filament formation process.The previous calculation based on density functional theory (DFT) has shown that the vacancyoxygen recombination process can occur almost spontaneously in sub-picoseconds. [22]Our results suggest that this intrinsic spontaneous recombination rate may be modulated or reduced  7) to (13).The parameter values used in each case from (a) and (e) are given in Table 1.f) Polarity dependent forming-voltage (V BD_FORM ) showing similar reverse area dependence for sample A at as more defects are generated in filament regions.This decaying function naturally arises from a physical consideration that interstitial oxygen density should decrease with time as they annihilate oxygen vacancies.It has been shown in the case of amorphous Si, traps can recombine with released interstitial hydrogen. [23,24]On the other hand, several DFT calculations [25][26][27] have been carried out for HfO 2 doped with hydrogen.It is shown that oxygen vacancy (V O ) positions can be filled by stable substitutional hydrogens (H O ). [25] Nearby interstitial oxygens can kick out these hydrogens (H O ) and recombine with oxygen vacancies to form hydroxyl groups as O o − H i defects. [25,26]These results suggest additional channels assisted by hydrogens for defect annihilation.Nevertheless, much more work is required to unveil the microscopic origin of the defect annihilation process as we propose in this work.
In our model, N CG is the critical defect number which was used in the previous work [5] for the modeling of conventional area-dependent BD data.Nevertheless, both N CG and N CA are not experimentally accessible quantities.Here we assume N CA to be equal to N CG for simplicity, so that we retain a single parameter N C in our model.This should not be mistaken to mean that there are equal concentrations of defects and annihilation centers.Note the introduction of the annihilation process is not intended to replace the generation process.Rather, the function of the annihilation process is to modulate the overall defect generation in a manner leading to a variety of T BD or V BD area-dependences, as is observed experimentally.In the case of reverse area-dependent BD, this time dependence of decaying annihilation centers in competition with defect generation plays a pronounced role in the final formation of the resulting filaments.Overall, our model contains six parameters, a 0 , N C , N G ,  G , N A , and  A .
Figures 4a-e exhibit the calculated CDF distributions of our model (Equation ( 13)) in comparison to experimental T BD_FORM data collected from the HfO 2 ReRAM devices of sample A, B, and C, showing excellent agreement.In the model calculation, we assume a value of 25 nm 2 for a 0 based on the work of conductingatomic force microscopy (c-AFM), which shows a range of 9-100 nm 2 for filament size. [28,29]This parameter is relatively insensitive to affect overall fitting results since the device area is much larger than the cell area, a 0 .Moreover, by sampling the parameter space, we arrived at an N C value of 8 as the best choice, (a value of 7 for a thinner bilayer film in Figure 4c).At a first glance, a value of 7 or 8 for N C might seem to be unusually small to be physically reasonable.First, this quantity represents an effective criterion from a statistics perspective to model filament formation, and this value is consistent with the previous modeling work. [5,15,16]Moreover, based on the published DFT calculation, [30] multiple vacancy complexes can be generated in the neighborhood of an original oxygen vacancy due to the energetic correlation effect.Therefore, N C in our model represents a reasonable value of the critical defect number for BD initiation.
Using the reliable experimental technique of isotope exchange annealing, the range of vacancy density is determined to be be-tween 2.8 and 40 × 10 17 cm −3 by a secondary ion mass spectrometry (SIMS) study [31] and further supported by the results of both DFT calculations [32] and molecular dynamics (MD). [33]This is in good agreement with the range of N G values in Table 1 we used to match the experimental data (Figure 4).The value of  G from our work for defect generation is smaller than the value of 0.4 previously estimated using current transient measurements in SiO 2 showing the conventional area-dependent BD. [15] Thus, this value corresponds to a situation where the annihilation process is absent or extremely weak in its time-dependence.These discrepancies can contribute to the smaller value obtained here which involves defect generation in bilayer HfO 2 .In addition, the values of N A and  A are also obtained by the best fitting of the experimental data.In general, we found that the N A values are larger than those of N G in the case of reverse area-dependent BD as shown in Figure 4c and listed in Table 1.
For HfO 2 dielectrics used in ReRAM devices, it is widely reported that the charge states of oxygen vacancies and interstitial oxygen are identified to be V 2+ o and O 2− I , respectively. [22,30]igure 4f provides a comparison of experimental data of V BD_FORM under the negative bias, suggesting that both V 2+ o and O 2− I participate in the annihilation process under the influence of an oppositely applied field.Recent work using the c-AFM technique to reconstruct 3D conductive filaments post-forming also show the movement of oxygen ions among multiple filaments in a ReRAM device. [29]These results are consistent with our model proposed in this work.In addition, the presence of multiple filaments found in Ref. [29] is also supportive of the notion of a defect competition process outlined in this statistical model (Figure 3a).Now, we will further explore the capability of this model below.Figure 5a.displays the calculated F(t) of Equation ( 13) using a set of parameters by only varying  A from 0.04 to 0.15 for three areas but with a fixed defect generation process.It is seen that the model first nicely yields a conventional area dependence ( A = 0.04), then a cross-over area-dependence ( A = 0.08) appears; and eventually, a reverse area-dependence ( A = 0.15) emerges, demonstrating the capability of our model to capture a variety of T BD_FORM area-dependences.In the case of conventional area-dependence, the appreciable F(t)'s appears at a much longer time because the defect generation process dominates over the defect annihilation process.On the other hand, as the defect annihilation process dominates, reverse T BD_FORM area-dependence emerges in a shorter time.These results clearly demonstrate that, while maintaining the weakest link character of filament formation as a result of defect generation, various features of T BD_FORM area-dependence can be realized through fabrication process innovation as we already demonstrated in this work.Figure 5b displays a series of the area-dependence curves of the normalized T BD_FORM at 50%.The available experimental data from Figures (1) and (4) are also included for comparison.These results demonstrate the capability of this model to produce a rich variety of area-dependences of T BD_FORM values.
a ramp rate of 0.4 V s −1 .More detailed information of samples used in this study and the measurement procedure for T BD_FORM and V BD_FORM data is provided in the experimental section.g) Some typical current transients measured using different areas for the bilayer dielectrics of 2 nm/2 nm shown in (c).The top four panels correspond to the currents measured at a stress voltage of 2.2 V whereas the bottom four panels for the currents measured at 0.1 V by interrupting the CVS stress for the sensitive detection of BD at low voltages.It is evident that the BD detections at 0.1V coincide with those at 2.2 V. Table 1.The parameters used in the model calculation for comparison with the experimental data in Figure 4 and the display of the model characteristics presented in Figures (3) and (5).The number in the first column for the corresponding figure numbers.N G and N A are the density of generated defects and annihilation centers at 1 s, respectively.Historically, the minimum value statistics (or the weakest-link theory) has been demonstrated to correctly capture the competition process among generated defects during filament formation and its manifestation in conventional BD area scaling (sometimes, called Poisson area scaling).The competition process not only captures the nature of dielectric BD or forming process but also defines its order statistical characteristics, unlike other statistical descriptions.From a physics perspective, an increasing defect density with increasing time (Equation ( 9)) is a natural outcome of the dominant role of defect generation in the competition process as demonstrated. [5]On the other hand, the reverse BD area-scaling is statistically consistent with a different statistical principle, namely the maximum-value statistics theory of Equations ( 5) and (6). [11,13,14]The rapid rise of leakage currents in Figure 4g for the samples showing the reverse areascaling effect suggests that filament formation still takes place in a competitive manner.In other words, the competition process is the driving mechanism for area-dependence of BD or forming because, otherwise, there would be no area-dependence of BD at all.Therefore, it is only logical to conclude there must exist an additional process that causes defects to be annihilated so that the BD area-scaling is reversed.This is physically consistent with the consideration that both concentrations of generated defects and annihilation centers increase with the increasing area, but the effect of annihilation centers reduces the effects of defect generation without suppressing filament formation due to the decaying time dependence of annihilation centers.The outcome of this competition gives rise to the final filament formation with a reverse area scaling.This further illustrates the important role of the competition process in understanding dielectric BD or forming process as well as using multiple areas rather than only the fixed areas for T BD_FORM or V BD_FORM measurements.In theoretical modeling of dielectric BD or forming process including kinetic Monte Carlo simulation, it is also critical to include defect competition mechanism which is generally neglected.

Conclusion
In summary, by carefully balancing between defect generation and annihilation via defect precursor introduction in the fabrication, our experimental findings show reverse area dependence of T BD_FORM and V BD_FORM using both CVS and RVS tests.A physicsbased statistical model involving both defect generation and annihilation processes is developed, showing good agreement with experimental data.This model is capable of producing a rich variety of area-dependences of T BD_FORM or V BD_FORM from conventional to reverse area-dependence.Specifically, a defect annihilation process with a decay function in time is shown to be essential for the modulation of defect generation process.As demonstrated in this work, the old inverse area-scaling law of dielectric BD is not inherently fundamental.Our work provides a key solution to semiconductor technology scaling and stimulates material innovation in computing technology as well as fundamental research in physics, statistics, and reliability with a wide range of potential applications.

Experimental Section
The detailed information for the fabrication process of the ReRAM stacks used in this study is previously presented in Refs.[18-20] as discussed its key aspects in Section 2.Moreover, the schematics in Figure 2a also compares the hydrogen-plasma-driven approach with the conventional reservoir-driven approach.An atomic layer deposition of 305 °C was used for the growth of HfO 2 thin films on 35-nm-thick TiN bottom electrode films deposited by a reactive sputtering.The hydrogen plasma treatment was used for both single-layer and bilayer structures as shown in Figure 2a,b.Note that this is fundamentally different from a conventional bilayer structure involving two different dielectrics.Particularly, for the bilayer structures, this H-plasma driven process [18,19] was adopted to achieve a defect-dominant layer of HfO 2 and followed by an in situ HfO 2 cap layer to protect the reoxidation of incorporated defects as shown in Figure 2a.Different times of H-plasma treatment was used to achieve various defect levels.As discussed above, since hydrogen can play a role in both defect generation and annihilation processes, the H-plasma treatment time was extended to 40 and 50 s for the samples used in this study as shown in Figure 2d beyond the previous study. [18,19]After the formation of the HfO 2 layers, 20-nm-thick TiN top electrode films were deposited by a reactive sputtering to complete ReRAM stacks as shown in Figure 2b.
The average values of forming voltages or V BD_FORM shown in Figure 2d were collected by using a ramped voltage stress (RVS) method with a current compliance of 100 μA.However, to study the interaction between defect generation and annihilation processes, an area dependence of T BD_FORM or V BD_FORM must be carried out rather than using devices with one area as commonly adopted in the literature.Based on the study reported in this work, it is clear that H incorporation not only reduces the forming voltage but also causes its reverse area-dependence, pointing to the important role of the defect annihilation process.
For electrical analysis of dielectric BD or forming process, the experimental measurements of T BD_FORM using the constant voltage stress (CVS) method were focused for better resolution to compare with the model, although both V BD_FORM and T BD_FORM measurements can be equivalently used to characterize dielectric BD. [11] This is because V BD_FORM measurements using RVS method are compressed in a small voltage window [11] despite its speedy data collection.Nevertheless, the RVS measurements for V BD_FORM data are also carried out to corroborate with the CVS measurements of T BD_FORM as shown in Figures 1a and 4f.
The electrical measurements of BD or forming data are collected using a source-meter-unit (SMU) controlled by an automated program in a personal computer.The data collection procedure for electrical measurements of T BD_FORM or V BD_FORM is well known and can be found elsewhere. [34]All the electrical measurements such as T BD_FORM or V BD_FORM as well as switching IV data are collected at a temperature of 30 °C in this work.Current compliance of 100 μA is always applied during the measurements to avoid catastrophic BD events.In addition, a low sensing voltage of 0.1 V by interrupting the stress is applied for more sensitive BD detection to avoid the situation where BD occurrence may be masked by higher voltage currents.This ensures the accuracy of our BD detection as evidenced in Figure 4g.The BD data in Figures 4a-f is defined where the stress currents reach 100 μA.It is well known that the statistical accuracy of BD measurements is susceptible to sample size variation. [34]o ensure the accuracy of experimental data, we use a large sample size from 63 to 68 devices in the most cases of CVS tests except for three sets of T BD_FORM data presented in Figure 4a,b,e, we used the sample size of 34, 50, and 33, respectively.
Although the main focus of this study is to study reverse area-scaling, we also carried out device switching measurements.Some typical reset and set curves for samples A, B, and C are provided in Figure S1 in the Supporting Information, showing these ReRAM devices with promising switching behaviors at the same time demonstrating reverse area-scaling characteristics of BD or forming.Figure S2 in the Supporting Information displays the cumulative probabilities of ReRAM resistance values measured at 0.1 V of low-resistance state (LRS) and high-resistance state (HRS) for sample A, B, and C, respectively, revealing a large window for memory applications.Again, a very large sample size of more one hundred was used for each sample type shown in Figure S2 in the Supporting Information.

Figure 2 .
Figure 2. a) Process flow and schematics for a reservoir-driven approach using reactive metal for reference process and hydrogen-plasma-driven approach in a single-layer and bilayer structures.The thickness values are given in nanometers whereas the numbers are in seconds for hydrogen plasma exposure times.The process details can be found in Ref. [18,19].b) Schematics of ReRAM device structures with the dielectric layer and top and bottom electrodes.The corresponding thickness values are also provided.The region indicated by HfO 2 is representative of either single-layer HfO 2 or bilayer HfO 2 structures discussed in the text and in the process flow above.c) The cross-sectional TEM image for the active device area.d) The V BD_FORM values with error bars for a wide range of process splits are compared using 10 × 10 μm 2 ReRAM devices with a ramp rate of 0.4 V s −1 at 30 °C.Each data point with the error bars in (d) represents the average value of ≈130 devices.Samples A, B, and C correspond to those showing the reverse area-scaling behavior as discussed in the text.The values of Samples A, B, and C represent the lowest values in comparisonwith other splits which are reported previously in Ref.[18].We included these results here for a broad view of the general trend and impact of H incorporation in forming voltage reduction.

Figure 3 .
Figure 3. Physics-based defect generation-annihilation statistical model.a) Schematics of competition involving defect generation and annihilation processes.b) The cell-based percolation statistical model incorporating both defect generation and annihilation processes.The yellow balls represent defects (O vacancy) and blue balls for annihilation centers (oxygen interstitials) as discussed in the text.The column enclosed by the bold red lines represents the filament formation.c) Time-dependence of annihilation center density,  A (t), in Equation (10) are plotted with three different N A values along with the time-dependent density of defect density,  G (t) of Equation (9).They are used for the model calculation to compare with the experimental data at different voltages in Figure4a,e.d) The CDFs, F u (t G ) and F u (t A ) of Equations (11) and (12) for the cell area (a 0 ) using the respective  G (t) and  A (t), with N G = 6.2 × 10 17 cm −3 and N A = 42 × 10 17 cm −3 at 1 s with  G = 0.08 and  A = 0.15, respectively.The corresponding PDFs, f u (t G ) and f u (t A ), are also included for comparison.The value of N C is 8 whereas t diel is 5nm with a 0 = 25 nm 2 as listed in Table1.The grey area represents the time-window over which the experimental data collected from larger device areas.e) The calculated joint PDF distribution of Equation (7) for three different areas using the cell distributions (d).f) The marginal PDF distribution of these three areas obtained by integrating Equation(7).The calculation results shown in (c),(f) eventually lead to the resultant CDF of Equation (13) to compare with the experimental data in Figure4a.Note that the dominant role of the defect annihilation process in conjunction with the defect generation process in (c),(d) in producing reverse area-dependence of T BD_FORM .

Figure 4 .
Figure 4.The experimental T BD_FORM data, showing reverse area-dependence, in comparison with defect generation-annihilation model.a-c) The CDFs of T BD _ FORM data for sample A at 2.4 V and B at 2.3 V at three areas, and for sample C using four areas at 2.2 V, respectively.d,e) The CDFs T BD_FORM data for sample A at lower stress voltages of 2.3 and 2.2 V, respectively.Note from (a) to €, the symbol represents the experimental data whereas the solid lines represent the calculation results of our defect generation-annihilation model from Equations (7) to(13).The parameter values used in each case from (a) and (e) are given in Table1.f) Polarity dependent forming-voltage (V BD_FORM ) showing similar reverse area dependence for sample A at

Figure 5 .
Figure 5.The calculation results from defect generation-annihilation model to demonstrate the transition from the conventional to the reverse areascaling trend of T BD_FORM .a) Theoretical calculation of the CDFs showing a variety of area-dependence behaviors varying from conventional areadependence to cross-over area-dependence, eventually to reverse area-dependence when the corresponding  A value varying from 0.15 to 0.08, and then to 0.04.The other parameters are given in Table 1.The arrows point to the right for conventional area-dependence while the arrows point to the left for reverse area-dependence.b) The normalized T BD_FORM values at 50% from the model calculation to demonstrate transition from the conventional area dependence of T BD_FORM to area-independence and then to reverse area dependence.The experimental data at 50% of different samples are included for comparison.The normalization is made to t = 10 s at the area of 10 4 μm 2 .