The modeling was carried out in the ZrO2 cubic (c-), tetragonal (t-) and monoclinic (m-) phase. Simulations were performed in terms of the hybrid DFT using a B3LYP functional in the Quantum-ESPRESSO code. The 192-atom c-, t-ZrO2, and 96-atom m-ZrO2 supercells were used to simulate the oxygen polyvacancy. Each subsequent atom removal from the supercell is determined from the point of view of the minimum formation energy. This model reflects a gradual increase in the polyvacancy size. The trap thermal Eth and optical Eoptionization energies are estimated by the formulas offered in .
The calculated band gap values of 5.7, 6.1 and 5.6 eV for c-, t- and m-ZrO2are consistent with the corresponding experimental values 5.1, 5.78, and 5.83 eV . The 4- and 3-fold coordinated (4f and 3f) O vacancy formation energies, 6.11 and 6.17 eV, predict the 4f vacancies prevalence in the real structure.
The Kohn–Sham levels of five VO charged states (q = -2, -1, 0, 1, 2) localized in the band gap indicate that any type of oxygen vacancies can capture both electrons and holes. The second added electron leads either to the singlet (V-2(s)) or a triplet (V-2(t)) defect state, and the latter is 0.06 eV lower in energy. The significant depth of negatively charged defect levels indicates a polaronic effect. The defect peaks width lower than 0.1 eV in TDOS spectra indicate the charge localization in space. The charge density spatial distribution of the excess charge also confirms the charge localization. The positive charge is distributed evenly between the nearest Zr atoms. An additional electron leads to a charge density distribution among two closest to each other Zr atoms, and one can say that there is a kind of Zr-Zr bond.
The calculated values of Eth can be interpreted as the energy gain when the conduction band bottom electron capture onto the defect. The positive Eth values indicate that the charge capture to the O vacancy in m-ZrO2 is energetically favorable. Thus, an oxygen vacancy in the ZrO2 is the amphoteric defect, it can act as the electron and hole trap, and participate in charge transport. The calculated values are consisted with Eth = 1.15  for ZrO2, and Eth = 1.25 eV, Eopt=2.5 eV for the HfO2 . For V0 and V+, a good agreement of the calculation and experiment is observed.
The O polyvacancy atomic structure in c-, t - and m-ZrO2 (Fig. 1) shows that each subsequent vacancy forms near the already existing one, and no more than the 2 removed O atoms relate to any Zr atom. One can assume that this result is valid for amorphous ZrO2, and the external electric field determines the preferred direction of the vacancies chains growth. This is consistent with the filament model of resistive switching. For t-ZrO2, the energetically favorable O polyvacancy atomic structure is a chain along the  direction.
Each subsequent removal of oxygen atom adds one filled level to the band gap. These levels are distributed on the band gap preferentially localized close to the monovacancy level. The oxygen polyvacancy growth in size yields a broad defective band in the energy spectrum of ZrO2 filled with electrons. Since, for the single oxygen vacancy, the extra electron forms something like a Zr-Zr bond, one can assume that ZrO2 can acquire conductivity, which is carried out through the channel of the Zr ions separating the vacancies in this polyvacancy chain.
The work was supported by the Russian Science Foundation, grant #16-19-00002. The modeling was performed on a cluster of the Institute of Semiconductor Physics SB RAS.
 X.D.Huang, R.P.Shi, P.T.Lai, Appl. Phys. Lett. 104, 162905 (2014).
 C.-H.Lai, H.-W. Chen, C.-Y.Liu, Materials 9, 551 (2016).
 D.Muñoz Ramo, J.L.Gavartin, A.L.Shluger, et al.Phys. Rev. B 75, 205336 (2007).
 R.H.French, S.J.Glass, F.S.Ohuchi, et al.Phys. Rev. B 49, 5133 (1994).
 G.Jegert, A.Kersch, W.Weinreich, et al.Appl. Phys. Lett. 96, 062113 (2010).