Despite extensive experimental and theoretical studies, the atomistic mechanisms responsible for dielectric breakdown (BD) in amorphous (a)-SiO_{2} are still poorly understood. A number of qualitative physical models and mathematical formulations have been proposed over the years to explain experimentally observable statistical trends. However, these models do not provide clear insight into the physical origins of the BD process. Here, we investigate the physical mechanisms responsible for dielectric breakdown in a-SiO_{2} using a multi-scale approach where the energetic parameters derived from a microscopic mechanism are used to predict the macroscopic degradation parameters of BD, i.e., time-dependent dielectric breakdown (TDDB) statistics, and its voltage dependence. Using this modeling framework, we demonstrate that trapping of two electrons at intrinsic structural precursors in a-SiO_{2} is responsible for a significant reduction of the activation energy for Si-O bond breaking. This results in a lower barrier for the formation of O vacancies and allows us to explain quantitatively the TDDB data reported in the literature for relatively thin (3–9 nm) a-SiO_{2} oxide films.

## I. INTRODUCTION

Dielectric breakdown (BD) is one of the most important phenomena determining semiconductor device reliability. It is characterized by an abrupt increase in current flowing through a dielectric layer, which typically happens when electric field exceeds the dielectric strength of the material, i.e., ∼15 MV/cm in a-SiO_{2}. As a result of BD, the silica layer loses its insulating properties; it is usually assumed that an O-deficient highly conductive region is formed, leading to an orders-of-magnitude increase in current and quasi-ohmic (depending on the current compliance) I–V characteristics.^{1}

Although BD in a-SiO_{2} has been investigated for more than 50 years, several aspects related to the atomistic mechanisms responsible for time-dependent dielectric breakdown (TDDB) are still unclear. In particular, the kinetics of the process, the atomic defects assisting in the current increase that occurs during BD, and the related structural modifications (i.e., the presumed creation of a highly oxygen deficient conductive path) are not fully understood.

Several models of BD in a-SiO_{2} have been proposed in the literature: the thermochemical E-model,^{2,3} the anode hole injection (AHI) 1/*E* model,^{4–6} the power low voltage model,^{7–10} and the exponential $ E 1 / 2 $ model.^{11–14} Each of these models explains some of the experimentally observable trends (i.e., field/voltage and temperature dependencies and Weibull statistics) of the time to BD (t_{BD}), which is the amount of time that the dielectric film can sustain a constant voltage stress without losing its insulating properties. However, these models are either empirical or based on over-simplified physical descriptions that do not properly address the microscopic complexity of the bond-breaking process, which can be locally affected by charge carriers, adjacent defects, and statistical variations of the bond properties. For example, in the thermochemical model proposed by McPherson, which is based on oxygen vacancy defect generation due to the Si-O bond breakage, the parameters describing the microscopic quantities (i.e., bond-breakage activation energy, field acceleration factor) are not fully justified and sometimes simply used as fitting parameters in order to reproduce the experimental data.^{15} In addition, the role of electron injection and trapping in the BD process has not been considered.

In order to connect the microscopic understanding of the BD process with quantitative description of TDDB, a more physics-based approach is required. This should allow one to explore the role of pre-existing defects and defect precursors in the BD process. A variety of atomistic models have been reported in the literature for the electrically active defects that support charge trapping and transport, and are responsible for MOSFET reliability mechanisms, such as bias-temperature instability (BTI)^{16} and dielectric breakdown.^{15,17} The role of electron trapping has been studied in a-Si_{3}N_{4} by combining *ab initio* calculations with a continuum-level transport model^{18} and charge trapping at defect sites has been shown to have a strong effect on the performance and reliability of electronic devices that employ a-SiO_{2} films as the gate insulator.^{19,20} These defects have typically been studied using a combination of density functional theory (DFT) and experiments. Hole trapping in silica has been modelled in previous studies^{21–24} and several hole trapping defects are well established.^{25–27} Electron trapping, on the other hand, is less understood and has often been attributed to hydrogen-related network fragments.^{28–31} It has been established that Ge impurities substituting for Si in *α*-quartz and a-SiO_{2} can trap both holes and electrons.^{32–34} Bersuker *et al.* suggested that electrons can be trapped by SiO bonds in the a-SiO_{2} network, and weaken these bonds to facilitate SiO bond dissociation.^{35} It has recently been shown that electrons can indeed trap at intrinsic precursor sites, such as wide O-Si-O angles, in a-SiO_{2} (Refs. 36 and 37) and create occupied states deep in the band gap. However, a single trapped electron is insufficient to break or significantly weaken a Si-O bond. It turns out that when such wide angle intrinsic sites trap two electrons, the Si-O bond dissociation activation energy is reduced significantly, which facilitates the formation of O vacancies and interstitial O ions.^{38} These studies suggest that both pre-existing defects and charge trapping play an active role in a-SiO_{2} film and MOSFET degradation.

Understanding the potential implications of this mechanism for device operation and reliability requires a multi-scale simulation framework that would allow us to use microscopic material and defect characteristics to predict macroscopic electrical device behavior. The purpose of this paper is thus to assess whether microscopic mechanisms of electron-assisted defect generation can be responsible for the BD of a-SiO_{2} by adopting such a novel multi-scale modeling approach. We focus on elucidating the role played by carrier injection and by existing atomic defects in the bond breaking process. Using this multi-scale model, we show that the bi-electron trapping at specific precursor sites in a-SiO_{2} results in a significant reduction of the O vacancy defect formation energy and thus facilitates BD. We demonstrate that based on this microscopic mechanism, one can reproduce not only the experimental TDDB data (along with their statistics and voltage dependence) reported in the literature for a wide range of thin a-SiO_{2} films (3–9 nm), but also the whole kinetics of the BD process.

The paper is organized as follows. The previously proposed models for describing BD in SiO_{2} are briefly reviewed in Section II. The multi-scale modeling framework combining *ab initio* calculations and device simulations is described in Section III. The novel oxygen vacancy generation mechanism proposed in this paper is described in Section IV, while the TDDB simulation results are discussed in Section V followed by Discussion and Conclusions.

## II. A REVIEW OF BD MODELS PROPOSED IN THE LITERATURE

To put this work in the context of previous studies, we briefly review the main existing TDDB models. These can be grouped into two categories: (1) Phenomenological models—which are based on empirical (mathematical/statistical) relations that can account for the TDDB dependence on experimental stress/electrical conditions. However, such models are not capable of identifying the physical mechanisms responsible for BD. (2) Physics-based models that can describe the defect generation rate and reproduce the TDDB data dependence on atomistic properties of materials. One key element that all models have in common is a description of the field (or voltage) acceleration factor

which indicates how fast TDDB changes upon the electric field (or voltage) application. Since TDDB experiments are typically conducted using voltage accelerated tests, correctly modeling the *γ* factor is of utmost importance for reliability projection under real device operating conditions. For this reason, TDDB models are also often identified by the relationship they provide between TDDB and the electric field. Below we provide a brief overview of the four BD models most often used for the interpretation of TDDB data.

### A. The thermochemical model

The thermochemical model attributes the generation of defects within the dielectric to the interaction between an external electric field and the inter-atomic bonds.^{2,3} The bond breakage rate depends on the bond vibration frequency and the probability that the chemical bond will receive enough thermal energy to be broken. Moreover, this model accounts for the role of the external field (*E*), which lowers the energy needed for the bond-breakage. These processes are captured by the following formula:

where *τ*_{0} is a constant related to the bond vibration frequency, $ \Delta H 0 $ is the zero-field activation energy of the bond-breakage process, *k _{B}* is the Boltzmann constant, and

*T*is the temperature. Although $ \Delta H 0 $ is linked to the nature of the chemical bonding in an oxide, this quantity has been often used as a fitting parameter rather than being directly related to the microscopic mechanism of the bond-breakage process.

The thermochemical model is also known as the E-model because it provides a theoretical foundation to the exponential dependence of TDDB on field, which has been observed experimentally at low stress voltages in thin dielectric films. Nevertheless, describing oxide breakdown as a field driven phenomenon does not explain the polarity dependence of BD in contrast to the experimental results that show different TDDB times for the same field in devices of different dielectric thicknesses.^{39}

### B. The anode hole injection (AHI) model

This BD model^{4–6} is based on the idea that holes generated at the anode are responsible for oxide damage. Under high-field electrical stress, electrons are injected into the oxide conduction band by Fowler-Nordheim (FN) tunneling and accelerated by the field towards the anode, where they can create holes through impact ionization. These holes may then tunnel back towards the cathode, and trap at some defect precursors (e.g., oxygen vacancies, forming $ E \u2032 $ centers) to create active defects (i.e., electron traps) responsible for dielectric degradation. Since electrons and holes are thought to be injected by FN tunneling, TDDB is described by an equation, whose field dependence resembles that of the FN current

where B and H are constants associated with electron and hole tunneling, while $ \tau 0 ( T ) $ is a temperature dependent pre-factor. Due to its characteristic dependence on the electric field, this model is usually referred to as the 1/*E* model.

Although this model is able to explain TDDB data in a wide range of experimental conditions, the AHI model has been criticized for some of its assumptions: (i) the defect generation rate due to holes is orders of magnitude higher than that due to electrons; (ii) the AHI model fails to explain the strong temperature dependence of TDDB, since the FN tunneling current is almost temperature independent. This model was mainly applied to relatively thick SiO_{2} layers, and it appears unsuited for thin oxides, where a too high a field (or voltage) would be required to generate holes at the anode.

### C. The anode hole release (AHR) model

The AHR model is based on the idea (experimentally verified) that hydrogen is involved in the generation of defects.^{7–10} Similar to the AHI model, electrons injected by tunneling can reach the anode with enough energy to release hydrogen atoms passivating dangling bonds at the Si/SiO_{2} interface. These positively charged H ions diffuse back through the oxide and can interact with other dangling and weak bonds, creating defects. The AHR model shows a TDDB power law dependence on the voltage

where *B*_{0} and *N* are parameters extracted from experimental data. One of the main weaknesses of the AHR model is that it does not explain the strong temperature dependence of TDDB. As with the AHI model, the AHR model cannot be applied to thin oxides, since a too high field would again be required to generate the high voltage needed for hydrogen release at the anode.

### D. The exponential $ E 1 / 2 $ model

The exponential $ E 1 / 2 $ model was originally developed for low-k silica dielectrics.^{11–14} It is based on the idea that the BD process is strictly connected to the charge transport mechanism and is considered to be assisted by pre-existing defects, whose density is relatively high in these materials (this is true also for other dielectrics, e.g., high-k). TDDB dependence on the field follows that of Poole-Frenkel or Schottky conduction (which is considered to be the dominant charge transport mechanism)

where *λ* is the root-field acceleration parameter and $\varphi $ is the barrier height. It is important to stress that the physics of the current-induced degradation has not yet been clarified. In addition, the $ E 1 / 2 $ model cannot explain the strong temperature dependence of TDDB.

### E. Discussion of BD models

The models described above are either empirical or based on an over-simplified physical description that does not properly address the microscopic complexity of the bond-breaking process that is responsible for BD, which can be locally affected by charge carriers, adjacent defects, and variations of local atomic structures. In order to provide deeper understanding of physical BD processes, an approach that is capable of connecting the atomistic material-dependent mechanisms to the macroscopic (measurable) electrical BD data is required.

The $ E 1 / 2 $ model lacks firm physical foundation; the voltage and temperature dependence of BD are assumed (without any conclusive proof) to be similar to models used for charge transport where these effects are attributed to either Poole-Frenkel or Schottky emission. In addition, the AHI and AHR models cannot be applied to relatively thin dielectrics such as those employed in state-of-the-art technologies, because the voltage drop is too low to generate a significant high energy population of the carriers responsible for either hole generation or hydrogen release.

On the other hand, in the thermochemical model proposed by McPherson^{15} (which assumes the Si-O bond breaking), the parameters connected to microscopic quantities (the bond-breaking activation energy, field acceleration factor) are often used as fitting parameters and the role of the injected charge and the proximity to existing defects (e.g., pre-existing O vacancies) is neglected.

Here, we use a novel simulation framework to understand the role played by electron injection and by existing atomic defects in the bond breaking process. This allows us to connect the microscopic phenomena occurring in the a-SiO_{2} network to the macroscopic electrical behavior observed at device level, i.e., voltage and temperature dependence of TDDB, and its Weibull statistics.

## III. MULTI-SCALE MODELING FRAMEWORK FOR SIMULATING SiO_{2} BREAKDOWN

### A. Simulation framework

Our device simulation framework benefits from an input from DFT simulations, as sketched in Figure 1. This multi-scale model self-consistently describes the main physical mechanisms present in SiO_{2} when the material is subjected to an electrical stress^{40} using the parameters generated by DFT calculations that explicitly consider the SiO_{2}-specific defect characteristics and defect generation mechanisms.

Charge transport is modeled self-consistently by including a variety of conduction mechanisms (relevant in dielectrics), such as direct tunneling, defect assisted contributions modeled in the framework of the multi-phonon Trap-Assisted-Tunneling (TAT),^{41} and carrier drift across either the conduction/valence band and defect sub-bands. The defect properties determined from the DFT analysis are used in the calculation of TAT current contributions accounting for the electron-phonon coupling: in particular, the defect thermal ionization and relaxation energies are the key parameters affecting the current and its temperature dependence.^{41}

The current distribution within the device volume is calculated consistently with the local potential, accounting for defect charge state and occupation. The power dissipation associated with the charge transport across the volume of the device is calculated by summing all contributions due to electron thermalization at defects and electrodes

where *R* is the rate of electron flow and $ \Delta E $ is the energy released by the electrons (either during trapping/detrapping events at defect sites or at electrodes).^{42} The calculated power dissipation is subsequently used to derive the temperature profile across the device by solving the Fourier's heat flow equation

where *k _{TH}* is the thermal conductivity of silicon dioxide.

In order to reproduce microscopic material modifications occurring during BD and the kinetics of the process, the simulation framework allows us to describe atomistic processes (using parameters from DFT calculations) that lead to creation of new defects, i.e., the breaking or distortion of Si-O bonds. The defect generation rates for the most relevant defect creation processes are implemented into the simulation framework accounting for the local field and temperature.^{42} This allows us to consistently model the field- and temperature-driven feedback that occurs during the BD process: the creation of new defects increases the leakage current, which in turn increases power dissipation, temperature, and the local defect creation rate. Once BD conditions (i.e., the presence of a minimal defect cluster^{42}) are triggered, this leads to a very fast acceleration of the defect creation process, culminating in the formation of a highly deficient filament responsible for the abrupt increase in current.^{1}

The stochastic nature of the process is accounted for using the Monte-Carlo method that determines the positions of the new defects as they are created. The field and temperature-driven motion of interstitial ions and vacancies is also statistically modeled consistently alongside charge transport. The complete multi-scale simulation framework which includes current, degradation, and breakdown simulations is included in the Ginestra™ commercial simulation package.^{40}

### B. Atomistic simulations

The defect creation is studied using a combination of classical molecular dynamics (MD) and density functional theory (DFT) calculations. The ReaxFF force-field^{43,44} and a classical MD melt and quench procedure are used with the LAMMPS code^{45} to generate starting structures representing non-defective continuous random networks of a-SiO_{2}.^{46} These structures containing 216 atoms each were evaluated^{46} by comparing the distributions of Si-O bonds, Si-O-Si angles, and neutron structure factors to prior theoretical and experimental studies before being used for DFT calculations. DFT calculations of the electronic and geometric structures and nudged elastic band calculations^{47,48} of adiabatic barriers were performed using the Gaussians and Plane Waves method^{49} implemented in the CP2K code.^{50} The PBE0-TC-LRC nonlocal functional^{51} was used in conjunction with the auxiliary density matrix method (ADMM)^{52} to mitigate the expense of using hybrid functionals. Finally, these calculations were performed using a 400 Ry plane wave cutoff, the Goedecker-Teter-Hutter (GTH) pseudopotentials,^{53} and a basis sets optimized from molecular calculations.^{54} Using this setup produced an average band gap of 8.1 eV for a-SiO_{2} with a range of 7.1–8.4 eV across a library of 500 structures.^{46} 10 periodic models containing wide angle intrinsic electron traps were then selected from this library with an average band gap of 8.1 eV and a range of 7.8–8.3 eV. The scope of DFT calculations in the multi-scale modeling scheme (see the green boxes in the flow chart in Figure 1) is to calculate: (i) electron-injection stimulated mechanism of formation of oxygen vacancies and interstitial oxygen ions and the associated energy barriers; (ii) the defect/precursor trap properties (i.e., thermal ionization and relaxation energies) required to calculate the electron transfer rates; and (iii) the field acceleration factor for the defect generation processes. These characteristics are then used to simulate the kinetics of the BD process and TDDB data (represented by the yellow boxes in the flow chart in Figure 1).

## IV. AN OXYGEN VACANCY GENERATION MECHANISM

The multi-scale modeling framework discussed above was employed to investigate the microscopic mechanism of BD that highlights the role of electron injection in the generation of defects in the SiO_{2} network.^{38} This mechanism stems from the recent discovery that extra electrons can be trapped in the a-SiO_{2} network and form deep electron states in the band gap of a-SiO_{2}.^{36} These trapping sites correspond to wide O-Si-O angles ( $ > 132 \xb0 $) in the otherwise continuum random network a-SiO_{2} structure and can accommodate up to two electrons. As a result, the energy barrier to break the Si-O bonds adjacent to the trapped bi-electron is lowered to around 0.7 eV on average.^{38} The detailed description of these electron traps, their optical absorption and EPR signatures, and electron injection-induced bond dissociation in a-SiO_{2} can be found in Refs. 36 and 38. Importantly, this bond breaking mechanism produces neutral O vacancies and negatively charged $ O 2 \u2212 $ interstitial ions which would not recombine easily. These $ O 2 \u2212 $ ions can diffuse away through the material via a previously studied pivot mechanism^{55} characterized by a low 0.3 eV energy barrier and have been observed experimentally to migrate rapidly towards the positive electrode and release into the gas phase.^{55} The estimated concentration of these intrinsic electron traps is at least 4 × 10^{19} in a-SiO_{2}.^{36} and they are used as a starting point for simulations performed in this work.

In order to incorporate this electron-trapping-assisted mechanism into the device model, one needs to account for the effect of the applied electric field on the energy barrier for Si–O bond breaking. For this purpose, we adopted the thermochemical model^{2,3,15} formalisms where the probability of breaking a Si–O bond is described as

where $ E A , 2 e $ (0.7 eV) is the activation energy required in the absence of the field for the irreversible oxygen ion displacement from its equilibrium position, while *k* is the relative dielectric constant of the SiO_{2} film. The effective dipole moment *p*_{0} can be estimated from the valence of the metal-ion (*Valence*), the distance (*d*) between ions, and (*n*) is the number of bonds aligned favorably with respect to the electric field with the strength *E* (Refs. 2, 3, and 15) as

The doubly occupied intrinsic electron trap can be characterized as a distorted SiO_{4} tetrahedron with a 178° average O-Si-O bond angle and two stretched Si-O bonds with bond lengths of 2.1 Å on average (see Figure 2(a)). When two electrons are trapped at a wide O-Si-O bond angle, only the two bonds participating in this wide angle are weakened. Of these two stretched Si-O bonds, only one can be oriented favorably with respect to the applied electric field at a time, as shown in Figure 2(a). Within the local geometry scheme described by McPherson,^{2,3,15} the dipole moment *p*_{0} can be related to the valence of the metal-ion, the distance between the metal and the O atom, and the number of bonds aligned favorably with respect to the applied field ** E**. For the intrinsic electron trap in a-SiO

_{2}, the Si metal-ion with two extra localized electrons can be assigned a valence of 2, the average Si-O bond length is 2.0 Å, and the maximum number of bonds that can be aligned favorably is 1. Using these parameters, the maximum dipole moment (

*p*

_{0}) is predicted to be around 2.0 eÅ. However, in the amorphous material, each SiO

_{4}is oriented differently with respect to the electric field, as illustrated in Figure 2(b). When the stretched bonds are perpendicular to the field, the number of aligned bonds approaches 0. This results in a distribution of

*p*

_{0}values ranging from 0 to 2.0 eÅ with an average value of 1.0 eÅ which is used for the BD simulations in Section V.

Newly created neutral O vacancies can then support trap assisted tunneling of charge carriers through the oxide. The properties of these traps have previously been discussed,^{41,56} highlighting that neutral O vacancies with a thermal ionization energy between 2.20 and 3.3 eV and a relaxation energy of 0.36 eV are responsible for trap assisted tunneling in a-SiO_{2}. In most situations (70% of 30 total calculated geometries), the neutral O vacancies can trap up to two extra electrons and are stable in the negatively charged state. The typical structure of the doubly negatively charged vacancy is shown in Figure 3.

## V. BREAKDOWN SIMULATIONS

### A. Model

We used the multi-scale simulation framework described in Section III to understand the role played by the charge-assisted defect-creation mechanism described in Section IV in SiO_{2} BD. We simulated the statistical TDDB distributions, their voltage dependence, and the BD kinetics in terms of the evolution of current, temperature, and defect concentration.

The charge-assisted defect creation mechanism was implemented into the device simulation using the thermochemical model formalism, which describes the defect (O vacancies) generation rate *G*

where *R*_{2} is the rate of the double electron capture process occurring at wide O-Si-O angle precursor electron trapping sites in a-SiO_{2} (which sets the maximum frequency of the process) depicted in Figure 2(a). Note that the local temperature *T* rises significantly (with respect to the external one) during the breakdown event,^{42} as discussed later. *R*_{2} is calculated by considering the double electron trapping as a sequential two-step Markov process, as illustrated in Figure 4(a). The Markov chain is comprised by three-states representing the wide O-Si-O bond angle precursor D in three different charge configurations: pristine (state 1), after trapping of an electron (state 2), and after the trapping of the second electron (state 3).

As discussed previously,^{16} the full transition rate from state 1 to state 3 of the double electron trapping process is given by

Here, $ \tau c , 12 $ and $ \tau c , 23 $ are the electron capture times from state 1 to state 2 and from state 2 to state 3, respectively, whereas $ \tau e , 21 $ is the electron emission time from state 2 to state 1.

A schematic presentation of the electron capture and emission processes and some of the corresponding notations is given in Figure 4(b). $ \tau c , 12 , \u2009 \tau c , 23 $, and $ \tau e , 21 $ are calculated for wide O-Si-O angle trapping precursors using the multi-phonon-TAT equations^{41}

$ E j , n $ is either the conduction band, valence band edge, or the energy level (E_{t}) of the jth trap; *N* and *f* are density of states and Fermi-Dirac occupation probability at either the cathode, anode, or the trap, respectively; P_{T} is the electron tunneling probability calculated using the WKB method. The capture, *Ca*, and emission, *Em*, rates accounting for carrier-phonon interactions^{41} are given by relations

Here, *c*_{0} is a constant that depends on the electric field and on the capture cross section of the trap,^{41} *ω*_{0} is the phonon frequency, and *m* and *n* represent the number of phonons exchanged during the trapping and emission processes (see Figure 4(b))

Here, *E _{t}* is the thermal ionization energy of the trap state in the gap with respect to the bottom of the conduction band of a-SiO

_{2};

*x*is the trap distance from the cathode and

_{T}*E*is the height of the energy barrier at the cathode or oxide interface, as illustrated in Figure 4(b)). Note that

_{B}*n*= 0 as emission occurs typically from the ground state of the trap.

^{41}

*L*is the multi-phonon transition probability

^{57}which accounts for the rearrangement of lattice atoms (the so-called relaxation process) required for accommodating (or removing) the electron charge into (or from) the defect during a capture (or emission) event

^{41}

*I _{m}* is the modified Bessel function of the order

*m*,

*f*is the Bose function which provides the phonon occupation number, and

_{B}*S*is the Huang-Rhys factor which represents the number of phonons required for the atomic-scale lattice rearrangement around the defect in order to accommodate the trapped charge. This lattice rearrangement process is described in terms of the relaxation energy,

*E*, that is related to the Huang-Rhys factor in Equation (18) as follows:

_{rel}*E _{rel}* is the intrinsic trap or vacancy relaxation energy due to electron trapping

^{41,58}that ultimately determines capture and emission time constants and sets the TAT temperature dependence. Equations (12) and (13) are used to calculate the time constants involved in both the electron trapping at wide O-Si-O bond angle precursors (Equation (11)) and in the TAT rate of electron flow through the new oxygen vacancies generated as a result of this trapping. As the increasing number of these vacancies is generated, they effectively assist the electron transport determining the observed current increase, which eventually leads to breakdown.

The defect thermal ionization energy, *E _{t}*, is calculated as the total energy difference between the defect state and the state where the defect electron is delocalized in the conduction band and the system geometry is relaxed.

^{59}The defect relaxation energy,

*E*, is calculated as the difference in total energies of the unrelaxed and relaxed defect state with extra electron(s). The parameters used for calculating the capture/emission rates on both wide O-Si-O bond angle precursors and oxygen vacancies are reported in Table I. We note the wide distributions of these parameters characteristic of amorphous structures. These are taken into account in breakdown simulations described below.

_{rel}Parameter . | Value . | Description . |
---|---|---|

Electron injection aided defect generation mechanism | ||

$ E A , 2 e $ | 0.7 ± 0.3 eV (Ref. 38) | Activation energy required to form a Frenkel defect pair |

p_{0} | 1 e Å | Wide O-Si-O angle defect dipole moment |

$ N T , O \u2212 S i \u2212 O $ | 4 x 10^{19 }cm^{−3} (Ref. 46) | Initial wide O-Si-O angle precursor concentration |

$ E t 1 , O \u2212 S i \u2212 O $ | 2.2 ± 0.2 eV | First electron capture thermal ionization energy |

$ E r e l 1 , O \u2212 S i \u2212 O $ | 1.5 ± 0.1 eV | First electron capture relaxation energy |

$ E t 2 , O \u2212 S i \u2212 O $ | 2.7 ± 0.2 eV | Second electron capture thermal ionization energy |

$ E r e l 2 , O \u2212 S i \u2212 O $ | 1.0 ± 0.30 eV | Second electron capture relaxation energy |

TAT charge transport mechanism | ||

Parameter | Value | Description |

E_{t} | 2.75 ± 0.55 eV (Ref. 38) | O vacancy thermal ionization energy |

E_{rel} | 0.36 eV | O vacancy relaxation energy |

$ \u210f \omega 0 $ | 0.06 eV | Phonon energy |

S | 6 | Huang-Rhys factor for O vacancies |

E_{B} | 3.1 eV | Electron tunneling energy barrier at the Si/SiO_{2} interface |

Parameter . | Value . | Description . |
---|---|---|

Electron injection aided defect generation mechanism | ||

$ E A , 2 e $ | 0.7 ± 0.3 eV (Ref. 38) | Activation energy required to form a Frenkel defect pair |

p_{0} | 1 e Å | Wide O-Si-O angle defect dipole moment |

$ N T , O \u2212 S i \u2212 O $ | 4 x 10^{19 }cm^{−3} (Ref. 46) | Initial wide O-Si-O angle precursor concentration |

$ E t 1 , O \u2212 S i \u2212 O $ | 2.2 ± 0.2 eV | First electron capture thermal ionization energy |

$ E r e l 1 , O \u2212 S i \u2212 O $ | 1.5 ± 0.1 eV | First electron capture relaxation energy |

$ E t 2 , O \u2212 S i \u2212 O $ | 2.7 ± 0.2 eV | Second electron capture thermal ionization energy |

$ E r e l 2 , O \u2212 S i \u2212 O $ | 1.0 ± 0.30 eV | Second electron capture relaxation energy |

TAT charge transport mechanism | ||

Parameter | Value | Description |

E_{t} | 2.75 ± 0.55 eV (Ref. 38) | O vacancy thermal ionization energy |

E_{rel} | 0.36 eV | O vacancy relaxation energy |

$ \u210f \omega 0 $ | 0.06 eV | Phonon energy |

S | 6 | Huang-Rhys factor for O vacancies |

E_{B} | 3.1 eV | Electron tunneling energy barrier at the Si/SiO_{2} interface |

### B. Results of simulations

Exploiting the statistical capabilities of the device simulations,^{40} we investigated the TDDB distributions by running 30 simulation trials of a thermally grown SiO_{2} stack in the same stress conditions, i.e., constant stress voltage and temperature. In order to test the accuracy of the simulation results, we considered previously studied p+poly/n-Si capacitors with a 2.7 nm-thick SiO_{2} film.^{17} For the O vacancy creation process, we included only the atomistic mechanism described in Section III. Wide O–Si-O bond angle precursors are randomly generated for every simulated sample with a uniform spatial distribution and with the energy parameters within the energy ranges reported in Table I for the electron capture processes.

The TDDB distributions simulated and measured at three stress voltages are shown in Figure 5. The experimental data are nicely reproduced by our simulations using the parameters for defect precursors and the defect creation processes determined using DFT calculations. Importantly, simulations reproduce correctly the voltage dependence of TDDB distribution, and also the TDDB slope, which is due to the intrinsic stochastic nature of the process, confirming the feasibility of the electron-injection-assisted oxygen vacancy formation model.

In order to further verify the model, we simulated the average of TDDB distributions measured at different voltages on a-SiO_{2} films of different thicknesses, Figure 6. The field dependence of the experimental TDDB data is well reproduced by our simulations.

In order to understand how changing the physical parameters of defect creation mechanism affect the simulation results, we run the simulations by varying $ E A , 2 e $, *E _{t}*, and

*p*

_{0}within reasonable ranges. Results presented in Figure 7 show that the variations of the model parameters from the values estimated from DFT calculations produce the results which deviate significantly from the experimental TDDB data. The high sensitivity of simulation results to model parameters indicates that the O vacancy defect creation process assisted by electron trapping at Si-O precursors is the likely microscopic mechanism responsible for BD in these a-SiO

_{2}films and allows one to accurately describe its statistical distribution and the voltage dependence.

The simulations performed using the multi-scale model allow us also to investigate quantitatively the kinetics of the BD process, understanding the evolution of the defect generation process culminating in the abrupt current increase at the BD. Figure 8 shows the current evolution simulated when applying a constant voltage stress of 3.0 V, along with the 3D map of generated defects monitored at different stages of the SiO_{2} degradation.

During the initial phases of BD, i.e., at time t_{0}, there are no defects in the SiO_{2} film: the current, dominated by direct tunneling, remains constant until about 130 s after the voltage stress application. During this initial phase, electrons are trapped at precursor sites in the film and new defects are generated almost uniformly across the oxide volume with a slightly higher probability of being generated close to pre-existing ones as a consequence of the local perturbation of the electric field induced by their charge state. The new defects assist charge transport through TAT, thus contributing to the current increase and to the associated power dissipation (time t_{1}) and temperature (see Equations (6) and (7)). This enhances the defect generation rate that becomes more localized in the proximity of the higher temperature oxide regions. This process continues until the random formation of a dominant defect cluster (comprised of around 25 vacancies with a mutual distance of no more than 6 Å (Ref. 42)) leads to a substantial increase of the local power dissipation and temperature increase by ∼20 K determining a fully localized defect generation in the surroundings of the hot spot. This event triggers a thermally driven positive feedback between current, temperature, and vacancy generation rates that quickly leads to the creation of a breakdown spot (time t_{2}) formed by an highly O-deficient region.^{42}

The kinetics of the process is dominated by the field- and temperature-induced positive feedback loop responsible for the current runaway,^{42} which can be controlled only by limiting the maximum current flowing through the film, i.e., the current compliance, which interrupts the positive feedback process setting the final size of the BD spot. Once the conductive O-deficient filament is formed, the dominant charge-transport regime changes. The TAT through the isolated O vacancies, dominating the early stages of degradation, is no longer valid when their local density exceeds the critical value of 10^{22 }cm^{−3}, corresponding to an average distance of 3.5 Å between defects^{63,64} (time t_{2} in Figure 8). As the local density of oxygen vacancies approaches the critical one, electrons get increasingly delocalized among adjacent defects eventually leading to the formation of a defect sub-band in which the dominant charge transport is the drift-diffusion mechanism. The effective conductivity is still calculated in the TAT framework by neglecting tunneling and lattice relaxation contributions (whose probabilities are set equal to one) and is consistent with theoretical values estimated according to the Landauer quantum conductance formula.^{63,65}

## VI. CONCLUSIONS

We proposed an atomistic microscopic mechanism responsible for the creation of O vacancies in a-SiO_{2}. Using a novel multi-scale simulation framework which relies on defect precursor properties and the activation barrier for the defect creation process obtained using DFT, we simulated the time-dependent dielectric breakdown (TDDB) distributions at different stress voltages. The good agreement between simulations and experiments confirms that this mechanism could be responsible for the degradation and dielectric breakdown in silica.

The adoption of such a microscopic, physics-based description of the mechanism(s) controlling degradation and breakdown allows a more accurate assessment of stress-induced dielectric degradation (with respect to existing models), thus enabling more reliable predictions of device reliability also at the statistical level. The proposed methodology is general and can be easily applied to other material systems for design-for-reliability applications. In particular, a similar electron-injection-assisted mechanism of Frenkel defect creation recently proposed in monoclinic HfO_{2} (Refs. 66–68) can contribute to BD in HfO_{2} based devices and this process will be considered in separate publication.

## ACKNOWLEDGMENTS

The authors would like to thank M. Morlini for help with simulations, and A.-M. El-Sayed and T. Grasser for useful discussions. We gratefully acknowledge funding provided by EPSRC under Grant No. EP/K01739X/1 Resistive switches (RRAM) and memristive behaviour in silicon rich silicon oxides. The authors acknowledge the use of the ARCHER High Performance Computing Facility via our membership to the UKs HPC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202).

## References

_{2}and Related Dielectrics: Science and Technology