Diffusion and aggregation of oxygen vacancies in amorphous silica

Using density functional theory (DFT) calculations, we investigated oxygen vacancy diffusion and aggregation in relation to dielectric breakdown in amorphous silicon dioxide (a-SiO2). Our calculations indicate the existence of favourable sites for the formation of vacancy dimers and trimers in the amorphous network with maximum binding energies of approximately 0.13 eV and 0.18 eV, respectively. However, an average energy barrier height for neutral vacancy diffusion is found to be about 4.6 eV, rendering this process unfeasible. At Fermi level positions above 6.4 eV with respect to the top of the valence band, oxygen vacancies can trap up to two extra electrons. Average barriers for the diffusion of negative and double negatively charged vacancies are found to be 2.7 eV and 2.0 eV, respectively. These barriers are higher than or comparable to thermal ionization energies of extra electrons from oxygen vacancies into the conduction band of a-SiO2. In addition, we discuss the competing pathways for electron trapping in oxygen deficient a-SiO2 caused by the existence of intrinsic electron traps and oxygen vacancies. These results provide new insights into the role of oxygen vacancies in degradation and dielectric breakdown in amorphous silicon oxides.


Introduction
The aggregation of oxygen vacancies as a result of electrically stressing amorphous SiO x (x = 1.3-2) films is thought to facilitate the dielectric breakdown of complementarymetal-oxide-semiconductor (CMOS) devices [1] and electroforming in resistive random access memory devices (RRAM) [2][3][4][5][6]. However, the fundamental mechanisms behind these processes are still poorly understood. Both the creation of additional oxygen vacancies [3,6,7] and the clustering of diffusing vacancies [2,5] have been proposed as mechanisms for vacancy aggregation. Here we focus on elucidating the feasibility of the latter process and use theoretical modelling to investigate the structures and binding energies of vacancy dimers and trimers in a-SiO 2 , energy barriers for individual vacancy diffusion, and whether this diffusion can be stimulated by the trapping of injected electrons at vacancies.
The aggregation of oxygen vacancies has been suggested as the mechanism for electroforming and resistive switching in RRAM devices employing several oxides, such as TiO 2 and NiO, as reviewed in [8]. The formation and disruption of oxygen vacancy conducting channels in these materials has been modelled in [9][10][11]. The calculated interaction profiles between vacancies indicated that specific defect arrangements are strongly favoured for oxygen vacancies in both TiO 2 and NiO [9,11]. The resistive switching mechanism of Ti/HfO x /Pt memory devices, studied using x-ray photoelectron spectroscopy and cross-sectional transmission electron S Supplementary material for this article is available online (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. microscopy, also suggested the significant dynamics of oxygen vacancies near both the top and bottom electrodes [12]. The segregation of oxygen vacancies at HfO 2 grain boundaries in TiN/HfO 2 /TiN based RRAM devices was proposed in [13,14] and dynamic flow between two oxygen vacancy reservoirs has been proposed for the set and reset mechanisms in RRAM devices employing amorphous HfO 2 films [15]. Two RRAM-based studies on monoclinic HfO 2 carried out by Bradley et al explored oxygen vacancy aggregation by diffusion [16] and additional vacancy creation by electron injection into existing vacancy sites [17]. The latter mechanism was favoured due to relatively weak vacancy binding energies when compared to higher barriers for diffusion.
In a broader context of radiation damage, catalysis and surface structure, the clustering of oxygen vacancies has been considered in the bulk and at surfaces of various other oxides, both theoretically and experimentally. Theoretical calculations suggested a relatively weak attraction between neutral oxygen vacancies in the bulk of MgO, CaO, and Al 2 O 3 , whereas repulsion was reported for ZnO [18,19]. A kinetic model proposed in [20] explained the formation of nanocavities in MgO by oxygen vacancy aggregation by diffusion in agreement with experimental observations. In addition, oxygen vacancy clustering has been predicted theoretically in SrTiO 3 [21] and BaZrO 3 [22]. The simulations by Ganduglia-Pirovano et al [23] demonstrated clustering of oxygen vacancies on VO 3 surfaces, which was confirmed by STM measurements. However, in the same study, CeO 2 surfaces showed no tendency for oxygen vacancies to aggregate. On the other hand, the experimental study in [24] has demonstrated a direct relationship between the concentration of the larger size oxygen vacancy clusters and the reducibility/reactivity of nanosized ceria nanorods. These results highlight the strong material and morphology dependence of oxygen vacancy aggregation.
Oxygen vacancy clustering is clearly a viable phenomenon in oxides; however, apart from the preliminary studies in [25,26], we are not aware of any systematic theoretical invest igations of this phenomenon in crystalline and amorphous SiO 2 in the context of defect studies or applications in electronic devices. A di-vacancy model has been proposed as the origin of the twofold-coordinated Si defect in silica glass in [27]. In this paper we therefore use DFT calculations to explore the possibility of oxygen vacancy aggregation by diffusion in a-SiO 2 . Assuming standard temperature and pressure device operating conditions, this could occur via two methods: (i) the formation and diffusion of vacancies in the oxide layer when an electric field is applied across the device, such as through the mechanisms described in [7] and [28], and (ii) the diffusion of pre-existing oxygen vacancies in siliconrich silica [2,5]. Our results demonstrate that vacancy aggregation is thermodynamically favourable only at some sites in the a-SiO 2 structure. As we show below, the calculated barriers for neutral vacancy diffusion are around 4.6 eV, which rules out their effective diffusion. Previous invest igations have suggested that electron injection into amorphous silicon oxides can accelerate the dielectric breakdown process [29][30][31][32][33]. It has recently been demonstrated that injected electrons can be trapped at intrinsic sites in the a-SiO 2 structure which results in deep electron traps [34,35]. Such intrinsic sites can trap up to two electrons and this weakens Si-O bonds and leads to the thermally activated creation of oxygen vacancies and interstitial oxygen ions [7]. These intrinsic sites compete for electrons with neutral oxygen vacancies, which can also trap up to two electrons. Our results demonstrate that electron trapping at vacancies facilitates diffusion by reducing barriers to about 2.7 and 2.0 eV for single and double charged vacancies, respectively. However, these barriers are still too high for the effective aggregation required for dielectric breakdown and electroforming in RRAM devices.

Theoretical methods
The ReaxFF [36] force field and the LAMMPS code [37] were used to generate initial models of a-SiO 2 . Classical molecular dynamics and a melt and quench procedure were applied to periodic cells containing 216 atoms, as described in detail in [34,35,38]. Starting from supercells with a cubic β-cristobalite structure, the system was equilibrated at 300 K and a pressure of 1 atm. Maintaining the pressure at 1 atm, the temperature was linearly ramped to 5000 K using a Berendsen thermostat [39]. The temperature was maintained at 5000 K for 40 ps to ensure a complete melting of the initial structure and then brought down to 0 K at a rate of 8 K ps −1 . The densities, distributions of Si-O bonds and Si-O-Si angles, and neutron structure factors of the obtained models are similar to those obtained in previous theoretical [40,41] and experimental [42][43][44][45][46] studies. An example structure is shown in figure 1.
Density functional theory (DFT) was then used to optim ize these starting amorphous geometries at 0 K, calculate their electronic structures, and study the properties of neutral and charged oxygen vacancy defects. In total five amorphous structures were selected for this study, with cell sizes ranging from 14.835 to 15.095 Å and densities ranging from 2.09 to 2.20 g cm −3 . The XYZ coordinate files of DFT optimized structures are included as supplementary information (stacks. iop.org/JPhysCM/29/245701/mmedia). A range of densities was used since small variations in density can arise from differences in fabrication conditions [45,46]. Additionally, the silicon-rich silica used in ReRAM studies [4][5][6] has been observed to exhibit variations in density within the oxide layer [6]. All DFT calculations were performed using the gaussians and plane waves (GPW) [47] method as implemented in the CP2K code [48]. In attempt to more accurately reproduce the band gap of a-SiO 2 , the PBE0-TC-LRC [49] non-local functional was used in conjunction with the auxiliary density matrix method (ADMM) [50] to reduce the computational expense. A long-range DZVP-MOLOPT basis set was used [51] along with the corresponding GTH pseudopotentials [52] and a 400 Ry plane wave cutoff. This resulted in a mean band gap of 8.0 eV for the a-SiO 2 cells as indicated in figure 2(a), which can be compared with approximately 9 eV from experimental measurements [53]. Lany-Zunger charge corrections were used for calculations on charged cells with a dielectric constant of 3.9 [54,55]. Finally, adiabatic barriers for diffusion were calculated using the climbing image nudged elastic band (CI-NEB) method [56] using 7 frames and a spring constant of − V 5.14 e Å 2 . In some cases where two maxima were potentially present in the 7 frame calculation, the calculation was divided into two separate 7 frame calculations in order to accurately characterise the reaction path.
The calculation of formation energies of vacancies, their aggregates, and vacancy diffusion barriers required a selection of oxygen vacancy sites which encompasses the range of local structural variations present in the amorphous material. In this study, pairs of nearest neighbour oxygen atom sites were selected so that they were situated within N mr -member rings where N 3 7 r m ⩽ ⩽ , with equal weighting given to each value of N mr . When the oxygen atoms were removed, these sites would leave a di-vacancy pair within the corresponding ring. Ring sizes were classified according to King's criterion [57] with silicon atoms being used as nodes. Tri-vacancy sites were then formed from a subset of di-vacancy sites with all values of N mr represented. A total of 54 neutral oxygen vacancies, 25 di-vacancy sites, and 12 tri-vacancy sites were obtained. As shown in previous studies, the removal of an oxygen atom leads to the formation of a Si-Si bond at the vacancy site with the Si-Si distance values ranging from 2.20 to 2.80 Å in the literature [38,[58][59][60]. The use of ring sampling was able to produce a comparable range from 2.25 to 2.70 Å.
30 sites from across this range of bond lengths were chosen to study the effects of electron injection and we aimed to calculate 10 diffusion barriers for the migration of vacancies to a nearest neighbour site with two barriers obtained for every value of N mr . The presented barrier heights describe the migration of the vacancy from the higher total energy configuration to the lower energy configuration. We note that, due to the even sampling of vacancy sites from across ring structures, the presented results do not intend to account for a statistically weighted distribution of sites. Figure 3 shows a histogram of the Si-Si bond lengths obtained using this sampling method in comparison to a previous study [61] with a-SiO 2 structures prepared using the same method [34,35,38] and where every oxygen atom site was sampled in a 216 atom cell. Both distributions peak at between 2.4-2.5 Å; however, this peaking is less prominent in the current study, which reflects the use of even ring sampling. Sites between 2.7-2.9 Å were not encountered in the current study although they form 0.06% of the 144 oxygen sites sampled by Wimmer et al [61].
Formation energies, E q form (− q 2 0 ⩽ ⩽ ), were calculated for neutral and charged oxygen vacancies using the following equation: (1) where E pristine is the total energy of the pristine cell, E defect is the total energy of the cell with an oxygen vacancy present, ∆E F is the Fermi energy level position with respect to the valence band maximum (VBM), E O2 is the total energy of an oxygen molecule, and E corr is a correction term to account for the spurious interaction between images of charged defects in periodic cells calculated using the method of Lany and Zunger [54]. E v is equal to the total energy difference between a cell with a hole introduced into the valence band in the dilute hole gas limit, and the corresponding pristine supercell [54]. For an oxygen vacancy cluster composed of n vac neutral vacancies, its corresponding formation energy, E n form vac , was calculated using the following equation: where E pristine is the total energy of the pristine cell, and E n defect vac is the total energy of the defected cell with n vac oxygen vacancies. In these calculations we assume that the Fermi level is suitably high that the neutral vacancies are stable. This is a reasonable assumption as, according to [62], the +/0 transition level for O vacancies in SiO 2 is around 2.5 eV above the top of the valence band, which is much lower than most relevant band offset energies of SiO 2 with Si or metallic electrodes (e.g. see [63]). At higher Fermi level positions (see also below), extra electrons can be trapped by neutral O vacancies and change their mobility. However, the vacancy diffusion process competes with thermal ionization of trapped electrons into the conduction band of a-SiO 2 . Thermal ionization energies were calculated for single electrons trapped at vacancy sites using DFT total energy values via the method described in [64,65].

Aggregation of vacancies
In order to gain an understanding of interactions between neutral O vacancies, we calculated the formation energies of the di-and tri-vacancy clusters per vacancy and compared them to the mean formation energy of the vacancies when present at the very same locations individually. Figure 4 shows a scatter plot comparing these values for each vacancy cluster where ( ) represent di-vacancy sites and (▴) represent tri-vacancy sites. Example electronic density of states plots for cells with one to three vacancies are shown in figures 2(b)-(d), respectively. A linear fit was performed to the scatter plot in figure 4 using least squares regression and resulted in a correlation coefficient of 0.92 with a gradient of 1.00. The corresponding residual plot revealed a random distribution of data points about the line of best fit, indicating a good fit to the data. These calculations suggest that, on average, there is no overwhelming trend for di-and tri-vacancies to cluster. The difference in energy between data points and the line of best fit can be interpreted as a measure of cluster binding energy for vacancies brought together from an initial infinite separation. 15 sites are positioned below the line of best fit, suggesting a gain in energy after vacancy aggregation relative to the sample average with maximum values of 0.13 eV and 0.18 eV for di-and tri-vacancies, respectively. In 8 of the 15 favourable sites, the vacancies are situated within 3-member rings. A typical relaxed configuration of such a cluster is shown in figure 5(a). The arrangement of Si atoms in a 3-member ring closely resembles that of the energetically favourable 'cyclic' Figure 3. Histogram comparing the Si-Si bond lengths at oxygen vacancies using the even ring sampling in this study as compared to a previous study where all oxygen sites in a 216 atom a-SiO 2 cell were sampled [61]. Both distributions peak at between 2.4-2.5 Å; however, this peaking is less prominent in this study, which reflects the even ring sampling used. Sites between 2.7-2.9 Å were not encountered in this study although they form 0.06% of the 144 oxygen sites sampled [61]. geometry of a Si trimer [66,67]. It should be emphasised, however, that 3-member rings represent a relatively small proportion of oxygen vacancy sites in the a-SiO 2 geometries (≈8%), with the proportion in real samples predicted to be even lower (≈0.2%) [68]. A tri-vacancy cluster can also have a 'Y' shape such that three Si atoms emerge outwards from a central Si atom ( figure 5(b)). Alternatively, three vacancies may also be arranged in a chain consisting of four Si atoms ( figure 5(c)). Four out of the 15 favourable sites in figure 4 correspond to a 'Y' arrangement, whereas tri-vacancy chains are all situated above the line of best fit. This agrees qualitatively with previous studies of the optimum geometry of four Si atom clusters which indicated that the 'Y' geometry is energetically favourable when compared to a chain of four Si atoms [66]. These observations could provide a qualitative explanation for the structure of energetically favourable sites for vacancy aggregation.

The migration of neutral oxygen vacancies
We now calculate energy barriers for the migration of neutral vacancies between adjacent sites to study whether they are mobile enough to aggregate at energetically favourable sites. Two diffusion mechanisms were observed and resulted in a mean migration barrier height of 4.6 eV. A prototypical migration mechanism is illustrated in figure 6(a) for N 4 mr ⩾ . The migrating O remains bonded to Si2 while diffusing towards the vacancy site, resulting in an increasingly stretched Si 1 -O bond (figure 6(a)(2)) and the transfer of electron density from the vacancy to between Si1 and Si2 ( figure 6(a)(4)).
An alternative migration mechanism is illustrated in figure 6(b) for = N 3 mr . The migrating O firstly bonds to Si3, with the electron density initially at the vacancy localizing on the under-coordinated Si1 ( figure 6(b)(2)). This O then reforms a bond with Si1 and the electron density occupies an anti-bonding orbital between Si1 and Si2 (figure 6(b) (3)). As the O settles into the vacancy position, the electron density finally occupies a bonding orbital between Si1 and Si2 ( figure 6(b)(4)). The associated adiabatic barriers for diffusion are summarized in table 1. In addition to an alternative migration mechanism, O vacancies in 3-member rings exhibited the lowest diffusion barrier of 3.2 eV. The calculated barrier heights suggest, however, that at standard temperature and pressure, the diffusion of neutral oxygen vacancies is unfeasible.

Electron traps in defect-Free a-SiO 2
In prior ab-initio studies, oxygen vacancies have been observed to behave as electron traps in a-SiO 2 [38,[69][70][71] with the Kohn-Sham (K-S) trap states located between 1.8 and 3.5 eV below the CBM [38,69]. The additional electron occupies an anti-bonding orbital as shown in figure 7(a).
Recent studies have demonstrated that extra electrons can also become trapped at intrinsic precursor sites in a-SiO 2 [7,34,35,72] which consist of wide Si-O-Si angles (⩾132°) [34,35], as indicated in figure 7(b). Trapping at these sites results in occupied K-S states located ± 3.2 0.05 eV below the CBM [34,35]. The number of precursor sites for spontaneous electron trapping in ideal a-SiO 2 has been estimated to be about × 4 10 19 cm −3 with a larger number of sites capable of electron trapping over an energy barrier [35]. The calculated thermal ionization energies of these traps into the conduction band range between 0.7 and 1.7 eV, dependent on the local environment [35]. It has been demonstrated that these electron traps can also explain the wide optical absorption band at 3.7 eV observed in silica glass samples irradiated at 80 K [72]. As shown below, an overlap in the distributions of intrinsic trap levels and those produced by oxygen vacancies leads to competition between these trapping sites for electrons injected into the system.

Negatively charged O vacancies
A bias voltage of several volts during electrical stressing can promote electron injection into the oxide either via electron tunnelling into the existing traps or via Fowler-Nordheim tunnelling and subsequent trapping at defect sites. Previous studies suggest that electrons trapped at oxygen vacancies in a-SiO 2 could accelerate dielectric breakdown [33]. To model the effect of electron injection, we added an extra one or two electrons to the simulation cell taking into account the correction for periodic images of charged defects. In 10 of the 30 sites chosen for the study, the injected electrons localized onto intrinsic trap sites [34,38] present in the cell before the formation of the vacancy. However, for the 20 remaining sites, the electron trapping occurred at the vacancies, despite the presence of intrinsic electron trapping precursor sites. This effect is correlated with the geometry of vacant sites. In particular, trapping at intrinsic sites appears to frequently occur for short respectively. It can be seen that at approximately 6.2 eV, single electron trapping becomes energetically favourable, while at 6.4 eV double electron trapping can become favourable at some vacancies. Furthermore, figure 9(a) shows a negative correlation between the 0/−1 charge transition level and the Si-Si bond length. This indicates that O vacancies with shorter Si-Si bond lengths tend to be the least favourable sites for electron trapping. In contrast, figure 9(b) suggests a positive correlation between the thermal ionization energy of a trapped electron and the Si-Si bond length, with ionization energies ranging from 1.0 to 2.2 eV with a mean of 1.6 eV. These thermal ionization energies strongly overlap with those for intrinsic electron traps discussed above (0.7-1.7 eV), which explains the competition for electron trapping observed in our calculations. This additionally suggests that electrons trapped at O vacancies with short Si-Si bond lengths would tend to ionize the most effectively. It should be emphasized  is not included for q = −2 since an energy minimum was found which was lower in energy than the vacancy trap states.  that the distribution of vacancy sites shown in the scatter plots are a product of the sampling method used in this study, and a more accurate distribution would require a weighted sampling method.
In CMOS devices, silicon oxides are typically interfaced with Si, with a resulting valence band offset of approximately 4.4 eV [63]. Additionally, in a previous RRAM study, silicon oxide was interfaced with TiN [6], which has a 4.2 eV [73] band offset with the SiO 2 CBM. A MOSFET with a Si gate typically operates at voltages close to 1 V, while previous RRAM devices with TiN electrode layers have been observed to switch at approximately 3 V [74]. Moreover, the initial electroforming step required a greater voltage of approximately 6 V. For Si and TiN electrodes, our results indicate that electron trapping at oxygen vacancies may be energetically feasible at voltages distributed around 2 V.

Diffusion of negatively charged O vacancies
To understand the effect of electron trapping on vacancy diffusion, we attempted to calculate adiabatic barriers for the diffusion of negatively charged vacancies with charges q = −1 and −2. Extra electrons were added to the cells containing neutral vacancies as discussed above. A total of 8 barriers were obtained from the 10 initial pairs of vacancy sites because in two cases the injected electrons localized at intrinsic electron traps in the cell. Our calculations indicate that for q = −1 the mean barrier height is reduced by 1.9 eV with respect to the neutral vacancy case. An analysis of the HOMO state occupied by the additional electron revealed that the migration of an oxygen atom towards the vacancy site results in an intermediate state which resembles an intrinsic trap [34,35]  In the three remaining cases, Si2 was observed to backproject [59,60] through the plane of the nearest O atoms (figure 10(f)(2)). The back-projected configuration formed an intermediate energy minimum between the initial and final configurations with 2 electrons localized on Si2 and Si3, respectively. This suggests a two-step diffusion process  which requires the initial back-projection of Si1 followed by the diffusion of O. The barriers corresponding to this reaction pathway formed the upper limit of the calculated values. For the 3-member ring, the back-projected state was found to be the lowest energy configuration, suggesting that diffusion of the vacancy would not be an energetically favourable process. The associated adiabatic barriers for diffusion are summarized in table 1. An average barrier height of 2.0 eV for double negatively charged vacancies is greater than a mean value of 1.6 eV for thermal ionization energy. Furthermore, in several cases the thermal ionization energies into the conduction band for the vacancies used in the barrier calculations are lower in energy than the diffusion barrier heights. These results suggest that thermal ionization of O vacancies with trapped electrons can compete with the diffusion process. Ionization would result in neutral vacancies, which are immobile.

Discussion and conclusions
We used computational modelling to investigate the structures and binding energies of vacancy dimers and trimers in a-SiO 2 , energy barriers for individual vacancy diffusion, and whether this diffusion can be stimulated by the trapping of injected electrons at vacancies. The calculations of di-and trivacancy clusters demonstrate that there are favourable sites for vacancy aggregation in a-SiO 2 with the maximum binding energies of approximately 0.13 eV and 0.18 eV, respectively. However, an average barrier for neutral vacancy diffusion was found to be about 4.6 eV, rendering effective clustering of randomly distributed neutral vacancies unfeasible. At Fermi level positions above 6.4 eV with respect to the top of the valence band, oxygen vacancies can trap up to two extra electrons. This greatly reduces the barrier for vacancy diffusion to about 2.0 eV. However, our results suggest that electrons trapped at vacancies would tend to thermally ionize before effective diffusion could take place, strongly diminishing the efficiency of this diffusion channel. These results provide an important contribution to our understanding of the role of oxygen vacancies in the mechanisms for dielectric breakdown and resistive switching in amorphous silicon oxides.
The result that the clustering of oxygen vacancies in a-SiO 2 via thermally activated diffusion is inefficient suggests that one should be looking for an alternative mechanism for the accumulation of oxygen vacancies under electrical bias. A recent study [7] shows that the trapping of two electrons at intrinsic electron traps can create Frenkel pairs consisting of an oxygen vacancy and O 2− ion with energy barriers of approximately 0.7 eV. This barrier could be further reduced by an electric field across a device. Furthermore, recent in situ electrical biasing measurements in combination with residual gas analysis (RGA) using a secondary ion mass spectro meter (SIMS) revealed the ejection of oxygen molecules from a TiN/a-SiO x /TiN stack [6]. Additional evidence for oxygen ejection was provided in the form of surface deformation of the top electrode observed in atomic force microscopy (AFM) and transmission electron microscopy (TEM) measurements. Such deformations have also previously been attributed to oxygen gas in titanium and silicon oxides [4,75]. The detection of oxygen emission during electrical stressing indicates oxygen vacancy and interstitial oxygen production. The mechanism proposed in [7] provides a plausible explanation for this phenomenon.
Finally, we note that our results demonstrate that trap assisted tunnelling in a-SiO 2 can involve both neutral oxygen vacancies and intrinsic electron traps which have close charge transition levels [35]. Depending on respective concentrations, both traps can contribute to an electron percolation path through the oxide. In several cases, we observed the trapping of electrons at intrinsic electron traps present nearby the vacancy site. These traps are present before the formation of the vacancy and can interact with electrons injected into the system in two ways: (i) as stepping stones for electron tunnelling through the oxide under bias; (ii) if two electrons meet at an intrinsic site, a new neutral oxygen vacancy can be created via the mechanism suggested in [7], again contributing to an electron current through the oxide. Our calculations suggest that longer Si-Si bonds form trap states with energies below the distribution of intrinsic traps. Electrical measurements by Mehonic et al [5] indicate that trap-assisted tunnelling is the major process responsible for the electron conduction in the soft breakdown path formed in RRAM devices. Further experimental studies of SiO 2 films with different concentrations of oxygen vacancies could shed further light on the mechanisms of these processes.