Complex Oxides under Simulated Electric Field: Determinants of Defect Polarization in ABO3 Perovskites

Abstract Polarization of ionic and electronic defects in response to high electric fields plays an essential role in determining properties of materials in applications such as memristive devices. However, isolating the polarization response of individual defects has been challenging for both models and measurements. Here the authors quantify the nonlinear dielectric response of neutral oxygen vacancies, comprised of strongly localized electrons at an oxygen vacancy site, in perovskite oxides of the form ABO3. Their approach implements a computationally efficient local Hubbard U correction in density functional theory simulations. These calculations indicate that the electric dipole moment of this defect is correlated positively with the lattice volume, which they varied by elastic strain and by A‐site cation species. In addition, the dipole of the neutral oxygen vacancy under electric field increases with increasing reducibility of the B‐site cation. The predicted relationship among point defect polarization, mechanical strain, and transition metal chemistry provides insights for the properties of memristive materials and devices under high electric fields.

In this study, we conducted all density functional theory and Berry phase simulations [1,2] with the Quantum Espresso package version 6.4.1 [3,4] . The plane-wave kinetic energy cutoff was set to 100 Ry and charge density cutoff to 400 Ry. Optimized norm-conserving Vanderbilt pseudopotentials, NCSR (ONCVPSP v0.4), from PseudoDojo [5,6] with Perdew-Burke-Ernzerhof exchange-correlation functional for solids (PBEsol) [7] are used in this study. Lattice constants for ABO 3 (A = Ca, Sr and Ba; B = Ti, Zr and Hf) were obtained using 21 different volumes, fitted with the 3rd order Birch-Murnaghan equation of state [8] . Reciprocal space was sampled using 8x8x8 displaced Monkhorst-Pack k-point grid [9] with no smearing in all calculations. All compounds were fitted under cubic symmetry constraint. Fitted parameters for all compounds are shown in Figure S1.
With the fitted lattice constants combined with our local Hubbard U method (details in Section 2) on ATiO 3 (A = Ca, Sr and Ba), we constructed 2x2x2 supercells for defect calculations. However, CaZrO 3 , CaHfO 3 and SrZrO 3 showed distorted structure by simply having an oxygen vacancy present without applying electric field; CaTiO 3 and SrHfO 3 showed distorted structure after applying very low field strength parallel to the 4+ − − 4+ chain (B = Ti and Hf). Since the cubic phase is the high temperature phase for these compounds [10][11][12] , the cubic symmetry constraint on these materials is unreasonable and made the supercells very unstable. By introducing an oxygen vacancy, the ions obtained more degrees of freedom, resulting in the distorted structures shown in Figure S2. In this work, we would want to conduct our simulations on materials under the same crystallographic symmetry. Therefore, we only selected SrTiO 3 (STO), BaTiO 3 (BTO), BaZrO 3 (BZO), and BzHfO 3 (BHO) in cubic phase to study their defect polarization responses under different electric field and strain. It is worth noting that we did not purposely enforce high structural symmetry to our materials except for BaTiO 3 . The structural ground state for BaZrO 3 , BzHfO 3 and SrTiO 3 are already cubic phase under room temperature [10,13] , and selecting cubic phase BaTiO 3 is necessary to provide reasonable comparison among cubic structures. We have shown the summary of the equilibrium lattice constant, bulk modulus, and band gap for the selected cubic phase perovskites in Table S1, along with the selected experimental values for each property. To ensure that the supercell size is big enough with minimal effect on our results, we have conducted additional simulation using a 3x3x3 SrTiO 3 supercell. As shown in Figure S3, both the relative electric Gibbs free energy of formation for oxygen vacancy and oxygen vacancy dipole moment in unstrained SrTiO 3 showed similar values and trend in different supercell sizes (2x2x2 vs. 3x3x3), indicating that 2x2x2 supercell is big enough for our calculation.
Previous work [14,15] showed that there are four possible configurations of neutral oxygen vacancy in SrTiO 3 . In our work, with spin polarized calculation and different magnetic moment initializations, we found that nonmagnetic state (singlet state) was the most stable configuration for such defect. Even with non-zero starting magnetic moment on the Ti ions adjacent to the defect site, the system would eventually relax to zero magnetic moment state. If we enforced magnetic solution for the defect (triplet state), the energy of ferromagnetic state was higher than the nonmagnetic state under all field strength as shown in Figure S4. Therefore, we focused on the ground state of neutral oxygen vacancy, singlet state with zero magnetic moment, in our current work. Amore recent study has used hybrid functionals [16] and showed that both electrons of neutral oxygen vacancy should be trapped in the in-gap defect state with a zero magnetic moment. With self-consistent DFT+U methods [15] and comprehensive surveys of all defect configurations, Ricca et al. have also shown that ferromagnetic vacancy is not always favored and one had to enforce Ti atoms adjacent to the vacant site at a specific distance to stabilize the ferromagnetic defect  [13] 186 179 [17] 2.46 3.25 [18] BTO 3.9725 4.00 [17] 180 162 [17] 2.35 3.2 [19] BZO 4.1826 4.1955 [10] 163 126 [20] * 3.73 4.9 [10] BHO 4.1491 4.1705 [10] 172 116 [21] * 4.18 5.4 [10] *Polycrystalline samples Figure S1. BM fitting for ABO3 (A = Ca, Sr and Ba; B = Ti, Zr and Hf) with 8x8x8 displaced Monkhorst-Pack kpoint grid [9] . 0 , 0 , 0 , 0 ′ and 0 are the equilibrium energy, equilibrium volume, bulk modulus, first derivative of bulk modulus with respect to volume, and equilibrium lattice constant, respectively.

b) Electric field calculation -Berry phase details
In this study, we conducted electric field simulations on STO, BTO, BZO, and BHO with a local Hubbard U method (for STO and BTO), Berry phase approach (lelfield = .true.) and modern theory of polarization [1,22] . Hubbard U was applied using the DFT+U simplified version by Cococcioni and de Gironcoli [23] , which corresponded to "lda_plus_u_kind = 0" in Quantum ESPRESSO. We used 2x2x2 supercell on both perfect and defective systems, with reciprocal space sampled using 2x2x2 Monkhorst-Pack k-point grid [9] , displaced by (0.5 0.5 0.5 structures were allowed to fully relax, while the shape and lattice constant were constrained. We set the convergence condition for electronic structure relaxation to be 10 -9 Ry, and ionic relaxation stopping criteria to be 4x10 -6 Ry for total energy and 4x10 -5 Ry/bohr for forces. We have shown the electronic density of states (DOS) for each compound under all field strengths in Figure S5. We adopted our previous method [24] , and applied small Gaussian smearing (0.004 Ry) to speed up our calculation which did not show any fractional occupations. Nevertheless, we tested the polarization behavior for unstrained STO with three different settings -2x2x2 k-point with smearing, 2x2x2 k-point without smearing, and 2x3x2 k-point without smearing (for Wannier calculations), which all showed identical polarization response under all fields as shown in Figure S6. We have also validated that all materials in our study did not undergo any phase change/distortion under the strain and electric fields considered in this study. We examined it by performing a backward relaxation (relaxing the structure from finite field to zero field), and observed no finite polarization from the structure at zero field. We limited the tensile strain on SrTiO 3 to be smaller than 2% (up to +1.3%) and on compressive strain on BaTiO 3 to avoid severe phase distortion.
In the main manuscript, we have shown the static permittivity of SrTiO 3 under all fields using local Hubbard U method with the lattice constant determined by DFT. Due to the smaller lattice constant obtained by DFT (3.8882 Å) compared to experiment measurement (3.900 Å), the static permittivity was underestimated. In Figure S7, we have shown the static permittivity of SrTiO 3 using local Hubbard U method with experimental determined lattice constant (green), which were extremely close to the measured values obtained from previous experimental work [25] (orange) under all fields.   (green) and measured values replotted from previous experimental work [25] (orange).

Local Hubbard U method a) Hubbard U values comparison
As discussed in the main text, for STO with oxygen vacancies, Hubbard U is necessary to obtain the correct in-gap defect states; otherwise, a system with vacancies would result in a questionable metallic solution under electric field [26,27] . In Figure S8, we have shown the integrated local DOS of the two electrons which occupied the highest energy states under zero field using traditional DFT method and our local Hubbard U method. It is clear that even without applying electric field, the electrons have spread to the Ti t 2g band, forming the 3d yz orbital shape.
Without Hubbard U on Ti 3d states, STO would form metallic solution under very low field strength. However, previous studies have also shown that the material polarization response under electric field would be greatly underestimated by the Hubbard U [28][29][30] .
In order to resolve this dilemma, we first tried combining Berry phase with hybrid functionals, hoping to achieve correct electronic and ionic structure at the same time. Nevertheless, although there was no theoretical reason that these two approaches should not work together, the current Quantum Espresso package structure was not implemented to use these two methods simultaneously. It appeared that the iteration processes of these two approaches interfere with each other, and would result in questionable results for self-consistent single-shot calculation, and would fail to converge for ionic relaxation [31,32] . Therefore, we returned to the DFT+U based method, and altered it by only introducing Hubbard U on the two Ti ions adjacent to the vacant site.
We first tested the Hubbard U values on Ti for localizing the defect in-gap state with 0.1, 1, 2.5 and 5 eV, where only 5 eV of Hubbard U value was able to provide integer occupation numbers, and lower Hubbard U values showed fractional occupation numbers. We then did a systematic survey on the influence of U values (3eV, 4eV and 5eV) with 2x2x2 k-point gird on different physical properties like polarization, , formation energies, band gaps and E CBM -E VO under different field strengths. As shown in Figure S9, different U values affect the defected system under low fields more obviously. However, the trend of the curves remained consistent, indicating that all our observations and analyses still hold independent of the selected U value.
To select the most adequate Hubbard U value in this study, we further tested the defected STO under electric field with denser k-point grid (from 2x2x2 to 2x3x2). Using Hubbard U values equal to 3eV and 4eV, we obtained metallic solution for defected STO under field strength 9MV/cm and 11MV/cm, respectively, while U value equal to 5eV showed consistent results up to 13MV/cm. Furthermore, recently self-consistent Hubbard U study on STO also suggested the ideal Hubbard U value on the neighboring Ti ions of the neutral oxygen vacancy being 4.8eV [15] .
With the equation describing the relationship between breakdown field and band gap for binary high-k dielectrics [33] , we could estimate the intrinsic breakdown field for STO being above 11MV/cm (although STO is not a binary oxide). Recent experimental study also showed that the breakdown field for BTO thin film can reach up to 8MV/cm in a heterojunction device [34] . From the estimations above, we believe the intrinsic breakdown field for STO should be above 10MV/cm in our simulation, since there was no geometric factor that could result in electric field concentration or tunneling in the simulation, causing the intrinsic breakdown field strength in simulation being much higher compared to experimental values. Moreover, the previous study also did a systematic survey of Hubbard U determination for transition metal binary oxides [35] , where a Hubbard U value equal to 5eV showed the best results on Ti ions. Although the pseudopotentials used in that study, GBRV [36] ultrasoft pseudopotential [37] , were different from the pseudopotentials used in our study, we believe it still provided important insight for Hubbard U value determination. Therefore, we selected 5eV to be our local Hubbard U value for STO and BTO in this study. Note that the local Hubbard U was used in both defected and perfect structures on the corresponding Ti ions to be consistent. Using the selected U values and parameters described in previous sections, we obtained band gaps of 2.46 eV for perfect STO, and 2.54 eV for defected STO. The underestimation of band gaps compared to the experimental value (3.2eV [19] ) was attributed to the lack of Hubbard U corrections on all the Ti and O ions in.
With 2 = 2 + 2 , where 2 being the oxygen chemical potential, 2 being the difference in chemical potential of 2 between the reference condition and the condition of interest, and 2 being the total energy of an isolated oxygen molecule computed by DFT, we set 2 to be zero to reflect oxygen rich condition, and obtained absolute oxygen vacancy formation energy under zero field of 5.41 eV. This is comparable to 6.2 eV under the same oxygen-rich conditions in ref [16] using HSE hybrid functionals.
For wide bandgap semiconductors (> 3.5eV from density functional theory simulation) such as BZO and BHO, local Hubbard U correction is not required to obtain correct in-gap defect state under facilitating subsequent application of the electric field. Nevertheless, we were curious about how the local U correction could affect the oxygen vacancy dipole moment even when such correction is not needed, especially when comparing to BaTiO 3 . In Figure S 10, we have shown the oxygen vacancy dipole moment for BHO under compression, and we found that the influence of local Hubbard U correction on such wide bandgap material is very small.
where U and P are the internal energy and polarization of defect-free (perf) or defected (def) SrTiO3 under electric field, respectively, is the chemical potential of oxygen, V is the volume of the supercell, and ⃑ is the electric field. The first part of the equation represents the relative internal energy of the system under electric field, and the second part of the equation represents the dot product of oxygen vacancy dipole moment, , and electric field.
In Figure S11, we have shown the relative Gibbs free energy along with its two components, relative internal energy and the dot product of and electric field. As internal energy increased relatively under electric field, the increase in , which resulting in the decrease in relative Gibbs free energy, dominated over the internal energy change. Therefore, is the main term determining the Gibbs free energy of neutral oxygen vacancy under electric field.
Figure S11 Relative electric Gibbs free energy for neutral oxygen vacancy in STO with field applied parallel to 4+ − − 4+ chain (red) with two components: relative internal energy (yellow) and • (blue) which increase and decrease relative electric Gibbs free energy, respectively.

Strain effect on dielectric and electronic properties in STO
Strain is known to influence the dielectric constant of materials significantly. In Figure S12, we have shown the strain effect on dielectric constant, polarization, bandgap, and E CBM -E VO with respect to electric field in STO. It is clear that tensile strain could increase the dielectric constant of STO significantly under low field, but the dielectric constant also dropped rapidly with respect to the field strength. Bandgap is known to have an inverse relationship with dielectric constant [38] .
In Figure S12

Electric field direction perpendicular to + − − + chainoctahedral rotation
When electric field was applied perpendicular to 4+ − − 4+ chain, defected STO started to show distorted structure at around 0.4MV/cm as shown in Figure S13. In all our calculations, we enforced cubic symmetry on all compounds. However, STO has tetragonal phase as stable structure below 105K [13] . Therefore, we believe that by introducing an oxygen vacancy, and applying the electric field along the undesirable direction (perpendicular to 4+ − − 4+ chain for STO) on defected STO, we could induce this structure distortion/phase transition. We then used the distorted structure from the 0.4MV/cm calculation, and decreased the field down to zero. We found that the crystal structure remained distorted even when it was fully relaxed under zero field, and assumed a polar phase. Although no such phase change/structure distortion was found in perfect STO under any field strength, we believe it would not be meaningful nor physical to compare the physical properties of the perfect cell with that of the defected cells in different phases. Therefore, we manually constructed a distorted perfect STO structure, by adding one oxygen ion into the distorted defected STO which was fully relaxed under zero field, and allowed the perfect structure to fully relax under zero field again. We then conducted the same process for electric field calculations on this distorted perfect STO. For the data analysis, since defected STO started to show distorted structure when electric field exceeded 0.4MV/cm, we also used the distorted perfect STO to calculate the formation energies and for field strength higher than 0.4MV/cm, and selected the normal perfect STO for field strength smaller than 0.4MV/cm. Therefore, there was a "stitch" at field strength 0.4MV/cm, which indicated the boundary between normal and distorted structures. In Figure S14,

Wannier calculation details and site-decomposed polarization
To trace the site polarization difference between a defected and perfect cell, we utilized maximally localized Wannier functions [39,40] with WANNIER90 version 3.1 [41] to obtain the sitedecomposed polarization from our Berry phase results generated from PWSCF, Quantum Espresso [3,4] . We only conducted Wannier calculations on STO and BHO since we believed these two compounds could already represent the trends for BTO and BZO, respectively. Similar to our previous work [24] , we used s-type Gaussians as our initial guess, and allowed projection centers to be guiding centers during the Wannierisation routine (guiding_centres = T). We also found that by increasing the density of the k-point grid, the accuracy for Wannier centers' position could be significantly improved under low fields, while not affecting the Berry phase results, as shown in Figure S6. Therefore, we used denser k-point grids for our Wannier calculations, 2x3x2 k-point grid for field applied parallel to 4+ − − 4+ chain, and 2x2x3 k-point grid for field applied perpendicular to 4+ − − 4+ chain. We selected the MLWFs numbers by targeting the valence bands in perfect structures, and both the valence and defect states in defected structures. We set our convergence threshold on Ω (spread of Wannier functions) to be 5x10 -8 (conv_tol = 5.0E-8) within two successive iterations (conv_window = 2).
In the following paragraphs, we will use site dipole to represent the dipole moment difference of    We also showed the comparison between the Berry phase approach and Wannier function for BHO below with electric field applied parallel and perpendicular to 4+ − − 4+ chain in Figure S22 and Figure S23, respectively. It is clear that the differences between the two methods are still very small under all fields, and the deviation in the plots seemed more obvious only because of the small in BHO. In Figure S24 and Figure S25, we also showed the site dipole between perfect and defected BHO. In contrast to STO, the site dipole in BHO is very small for each site, and the curves basically followed the trend of vacant site dipole curves.     b) Compressive strain effect on negative under low field As discussed in section (a), several sites in the defected cell became less polarizable than the corresponding sites in the perfect cell, which contributed negatively to (negative sites). This phenomenon is more prominent under low field strength, and compressive strain played an important role in affecting it, which caused the slight difference in between compressed 2.1% BTO (always positive) and unstrained STO (small negative value). Compressive strain is known to decrease the polarization of ions [42] . Therefore, under compressive strain, the ions in both defected and perfect cell would all became less polarizable, and eventually become almost immobile under extreme compressive strain, where we can see them being inert to external electric field. In Figure S26    However, the main source promoting is the Hf1 ion instead of the trapped electrons. As shown in Figure S29     In Figure S31, we have shown the field induced dipole moments for different defect clusters with respect to the B site cation electronegativity. We believe the most adequate defect cluster to directly connect with B site electronegativity is the 4+ − − 4+ chain, as shown in Figure   S31 (a). In Figure S31 (b), we can see the field induced dipole moment for single oxygen vacant site does not perfectly correlate with B site cation electronegativity. This discrepancy mainly resulted from that we should not correlate electronegativity with a "single site", but to connect it to a "bonding state". The difference in electronegativity between two elements have been widely used to described the covalency of the bond between two different elements. Although there was no bonding between the trapped electrons and the neighboring B site cation, they still affect each other's polarization simultaneouslya two-way influence. Therefore, it would be unphysical to isolate the oxygen vacant site alone when correlating to B site cation electronegativity; instead, we should include the whole defect cluster to capture the whole picture. In Figure S31 (c), we have also shown the − 4+ − − 4+ chain with respect to the B site cation electronegativity. The trends in Figure S31 (a) and (c) are almost identical, and − 4+ − − 4+ actually contributed to slightly more than 4+ − − 4+ alone. Nevertheless, the oxygen in between the B site cations does not directly interact with the trapped electrons in the vacant site, where it was indirectly affected by them under electric field. Furthermore, the contribution difference between two different defect clusters was also very small. Therefore, in order to have a simpler and more direct connection between the defect cluster dipole moment and B site cation electronegativity, we believe 4+ − − 4+ chain should be the best defect cluster we should select. Other than the electronegativity of the B site cation, we also tried to correlate the defect polarization to a more material-specific quantity -the effective mass of electron at conduction band minimum (CBM). We believe that the more reducible the B site cation (equivalent to higher electronegativity), the more easily it can "trap" the electron, which will lead to larger electron effective mass. Therefore, defect cluster dipole and oxygen vacancy dipole moment should also correlate positively to electron effective mass.
To validate the hypothesis, we calculated the band structures and the corresponding electron effective mass at CBM along two high symmetry reciprocal space directions (Γ to Χ and Γ to Μ) [43] for BaTiO 3 , BaZrO 3 and BaHfO 3 with all materials having the same lattice constant (lattice constant of unstrained BaTiO 3 ) to exclude the effect of lattice volume. As shown in Figure S32, the trend of defect cluster dipole and oxygen vacancy dipole moment both correlate positively to electron effective mass regardless of band and k-space direction (Γ to Χ or Γ to Μ). Such trend is identical to that being plotted against B site cation electronegativity ( Figure 4 in the main manuscript).
Finally, since the − − of BaTiO 3 and Ba(Zr, Hf)O 3 responded differently to lattice strain, where the former increase with tensile strain, and the latter increase with compressive strain, we did a rough interpolation and found a possible transition region where the strain should have no effect on − − around B site cation electronegativity being around 1.4. Interestingly, there is no element having electronegativity being exactly 1.4, where only uranium (U) with electronegativity being 1.38 [44] is close enough to fall into this transition region. Figure S32. Defect cluster dipole and oxygen vacancy dipole moment with electric field applied parallel to 4+ − − 4+ plotted against electron effective mass at CBM. Electron effective masses were obtained from heavy electron (he) band and light electron (le) band, with two high symmetry directions in reciprocal space, Γ to Χ and Γ to Μ [43] .