Unveiling Temperature-Induced Structural Domains and Movement of Oxygen Vacancies in SrTiO3 with Graphene.

Heterointerfaces coupling complex oxides exhibit coexisting functional properties such as magnetism, superconductivity, and ferroelectricity, often absent in their individual constituent. SrTiO3 (STO), a canonical band insulator, is an active constituent of such heterointerfaces. Temperature-, strain-, or mechanical stress-induced ferroelastic transition leads to the formation of narrow domains and domain walls in STO. Such ferroelastic domain walls have been studied using imaging or transport techniques and, often, the findings are influenced by the choice and interaction of the electrodes with STO. In this work, we use graphene as a unique platform to unveil the movement of oxygen vacancies and ferroelastic domain walls near the STO surface by studying the temperature and gate bias dependence of charge transport in graphene. By sweeping the back gate voltage, we observe antihysteresis in graphene typically observed in conventional ferroelectric oxides. Interestingly, we find features in antihysteresis that are related to the movement of domain walls and of oxygen vacancies in STO. We ascertain this by analyzing the time dependence of the graphene square resistance at different temperatures and gate bias. Density functional calculations estimate the surface polarization and formation energies of layer-dependent oxygen vacancies in STO. This corroborates quantitatively with the activation energies determined from the temperature dependence of the graphene square resistance. Introduction of a hexagonal boron nitride (hBN) layer, of varying thicknesses, between graphene and STO leads to a gradual disappearance of the observed features, implying the influence of the domain walls onto the potential landscape in graphene.

Introduction

Transition-pan class="Chemical">metal oxides are versatile material systems that offer diverse functionalities, such as ferromagnetism, superconductivity, ferroelectricity, multiferroics, etc.[1] A prototypical transition-metal oxide is SrTiO3 (STO), which is a canonical band insulator and has a large temperature-dependent dielectric permittivity.[2] It is known to possess a temperature-driven ferroelastic transition at 105 K[3] and a quantum paraelectric state persisting to lower temperatures,[4] accompanied by differently oriented structural domains that move under the application of an electric field.[5−7] The moving domain walls can lead to large local polarization at the STO surface.[8−10] Further, oxygen vacancies tend to accumulate at STO domain walls and act as channels for the movement of oxygen vacancies that can dramatically influence the dielectric environment both with temperature and with electric fields.[11,12] This intricate interplay between structural transition and electronic transport has been studied by several groups using different macroscopic and transport probes involving different levels of complexity. The domain walls were found to host interesting properties, such as field-induced ferroelectricity, high mobilities up to low temperature, and enhanced conductivity.[5−7,13]

In proximity with pan class="Chemical">STO, severpan class="Chemical">al electronic properties in graphene can be influenced, predominantly due to the large dielectric permittivity and strong correlation effects in STO.[14−16] STO can influence electron–electron interaction in graphene as manifested in the temperature dependence of both charge transport[17−19] and spin transport and is ascribed to the structural phase transition and presence of surface dipoles[20] in STO. Further, antihysteresis in graphene square resistance, when back gated through STO has been reported by several groups and is an active area of research.[17,19,21]

In this work, we interface pan class="Chemical">STO with pan class="Chemical">graphene to unveil the movement of oxygen vacancies and ferroelastic domain walls near the STO surface by systematically studying the temperature and gate bias dependence of charge transport in graphene. For this purpose, we track the variation of the charge neutrality point (CNP) with temperature and direction of the applied gate bias using a four-probe geometry. By sweeping the back gate voltage, we observe antihysteresis in graphene that is typical of conventional ferroelectric oxides, unlike STO, which is an incipient ferroelectric. Our study on two different fabricated stacks of graphene on STO and graphene on hexagonal boron nitride (hBN)/STO shows additional features in the (anti)hysteresis traces of the square resistance in graphene. By inserting hBN layers, we disentangle the bulk effect from the interface effect.

From our anpan class="Chemical">alypan class="Chemical">sis, we infer that these features are related to the movement of domain walls and their interplay with oxygen vacancies close to the STO surface and depend on the direction of gate voltage sweep. We also find that the separation between the CNP peaks, for the trace and retrace cycles, is proportional to the sweep range of the gate voltage and inversely proportional to the temperature, imprinting the structural phase transitions in STO. Further, density functional theory (DFT) studies of the layer-resolved STO surface calculate the formation energies of the oxygen vacancies and find it to be lower at the surface than in the bulk. Quantitative estimates of the activation energies of the mobile oxygen vacancies were determined from the temperature dependence studies of the graphene square resistance and were found to agree with those obtained from DFT studies. Our study establishes graphene to be a pristine platform to unveil the unusual temperature dependence of domain wall motion and oxygen vacancies that are coupled to the ferroelastic transition in STO.

Results and Discussion

To fabricate the devices for our study, we transferred a monolayer of pan class="Chemical">graphene with or without pan class="Chemical">hBN of different thicknesses onto a 0.5 mm thick bulk STO using the dry-transfer technique.[22] This was followed by contact electrode deposition. A schematic diagram of the device is shown in Figure a. The thickness of the hBN was characterized by atomic force microscopy (AFM) and is shown in the Supporting Information (SI) (Figure S1). The graphene channel resistance was monitored using four-probe measurements so as to exclude the contact resistance. The 0.5 mm thick STO with or without hBN served as the gate dielectric. We compare the charge transport in graphene in direct contact with STO with that when hBN is inserted between graphene and STO. This allows us to examine the contribution of graphene/STO interface on the charge transport. In total, three devices were fabricated: graphene/STO (device 1), graphene/hBN (8 nm)/STO (device 2), and graphene/hBN (23 nm)/STO (device 3). Figure b–d shows the back gate dependence of the square resistance for different sweep ranges and temperatures for device 1 at the sweeping rate of 1.7 s/V. At 150 K, as shown in Figure d, the graphene square resistance shows a typical Dirac curve without any antihysteresis. However, with decreasing temperature, the dielectric constant of STO increases significantly,[2] leading to a narrowing of the Dirac curves and a large antihysteresis, as shown in Figure b,c. At positive gate voltages, graphene is initially populated with electrons. During retracing (positive to negative gate voltage sweep), the CNP in graphene is reached even with a small gate sweep range, leading to a positive shift of CNP. Likewise, during tracing, the CNP shows a negative shift. The separation between the peaks at the CNP, for the trace and retrace cycles, is found to be proportional to the sweep range of the back gate. Interestingly, we observe features in the Dirac peaks for all sweep ranges and for temperatures below 105 K. Such features in the antihysteresis vanishes for temperatures above 105 K.

Figure 1

Antihysteresis in graphene square resistance on STO (device 1). (a) Measurement schematic: a four-probe measurement scheme is used, which excludes the contact resistance. The back gate is applied through a 0.5 mm STO single crystal substrate. The antihysteresis behavior for device 1 was studied at different temperatures for different gate sweep ranges with a constant sweeping rate of 1.7 s/V. The gate sweeping range was systematically changed from ±20 V (black) and ±50 V (red) to ±80 V (blue) for the graphene/STO device at (b) 4 K, well below the phase-transition temperatures, (c) at 105 K, the ferroelastic phase-transition temperature, and (d) at 150 K, well above the phase-transition temperatures.

Antihysterepan class="Chemical">sis in pan class="Chemical">graphene square resistance on STO (device 1). (a) Measurement schematic: a four-probe measurement scheme is used, which excludes the contact resistance. The back gate is applied through a 0.5 mm STO single crystal substrate. The antihysteresis behavior for device 1 was studied at different temperatures for different gate sweep ranges with a constant sweeping rate of 1.7 s/V. The gate sweeping range was systematically changed from ±20 V (black) and ±50 V (red) to ±80 V (blue) for the graphene/STO device at (b) 4 K, well below the phase-transition temperatures, (c) at 105 K, the ferroelastic phase-transition temperature, and (d) at 150 K, well above the phase-transition temperatures.

Both hysteresis and antihysteresis have been observed in graphene with different ferroelectric media.[23−29,30,31]Antihysteresis, specifically, has been attributed to either the presence of ferroelectric-like surface dipoles or to dynamical trapping in the oxygen vacancy sites as well as to band bending at the STO surface due to the surface dipoles that act as the trapping sites.[17,19,21] However, the cause of such antihysteresis is still under debate. Additionally, features in the Dirac peaks were associated with different CNPs in graphene, originating from contact doping,[32] local electrostatic doping from charge traps,[33] p–n junction formed by local gating,[34] or from local dopants from adsorbates.

In this context, our observations on the three typn>es of devices studied are distinctly different from earlier works. n>an class="Chemical">Although we observe antihysteresis for all of the three devices, the multiple features in the Dirac peak appear only for device 1 and 2, as shown in Figures b and 2a, and for temperatures lower than 105 K, while with a thick hBN layer (device 3), such features do not occur at any temperature Figure d–f.

Figure 2

Antihysteresis in graphene square resistance on hBN/STO for different gate sweep ranges with the sweeping rate of 1.7 s/V. (a–c) The antihysteresis curves for the device with 8 nm hBN (device 2) at different temperatures: (a) 4 K, well below the phase-transition temperatures, (b) 105 K, the ferroelastic phase-transition temperature, and (c) 150 K, well above the phase-transition temperatures. The gate sweeping range was systematically changed from ± 20 V (black) and ± 50 V (red) to ± 80 V (blue). (d–f) Curves for the device with 23 nm hBN (device 3) at 4, 105, and 150 K, respectively.

Antihysteresis in n>an class="Chemical">graphene square resistance on hBN/STO for different gate sweep ranges with the sweeping rate of 1.7 s/V. (a–c) The antihysteresis curves for the device with 8 nm hBN (device 2) at different temperatures: (a) 4 K, well below the phase-transition temperatures, (b) 105 K, the ferroelastic phase-transition temperature, and (c) 150 K, well above the phase-transition temperatures. The gate sweeping range was systematically changed from ± 20 V (black) and ± 50 V (red) to ± 80 V (blue). (d–f) Curves for the device with 23 nm hBN (device 3) at 4, 105, and 150 K, respectively.

We carried out several n>an class="Chemical">cooling cycles and traced out similar curves also for other contact electrodes on graphene. Such features were distinctly observed at 4 K, while they disappeared between 60–105 K. Moreover, the exact shape of these features in the Dirac peaks are found to be different in the different regions probed, as well as for different thermal cycles through 105 K, as shown in SI, Figure S3.

This temperature-dependent behavior of the multiple features in the pan class="Chemical">Dirac peak and for the three different devices can be understood from a concerted interplay of associated factors in STO. Besides the enhancement of the dielectric permittivity in STO at low temperatures, the temperature-induced phase transition (around 105 K) in STO triggers the formation of a network of structural domains and domain walls that are mobile and polar down to low temperatures,[5−10,35] and act as trap centers for the oxygen vacancies in STO.[12] These not only change the electrostatic environment in graphene but also generate local gate fields, which assists or opposes the formation of the observed features in the Dirac peak for different polarities and range of voltage sweeps.[36,37] The addition of an hBN layer (device 2) reduces the effective capacitance due to a reduced permittivity of the stack. The total capacitance Ctot per unit area of graphene/hBN/STO can be roughly estimated from the geometrical capacitance of hBN (ChBN) and STO (CSTO) in series:Ctot–1 = ChBN–1 + CSTO–1, where the geometrical capacitance can be estimated by C = εrε0/d when the quantum capacitances are neglected, with εr being the dielectric constant, ε0 being the vacuum permittivity, and d being the thickness of the dielectric. Assuming that εr for STO is 10 000, and εr for hBN is 4, the geometrical capacitances for devices 1 (dhBN = 0 nm), 2 (dhBN = 8 nm), and 3 (dhBN = 23 nm) are estimated to be 2.00 × 107ε0, 1.92 × 107ε0, and 1.79 × 107ε0. This leads to two effects: it reduces the overall electric field from the external gate bias and it reduces the influence of these local gate fields on the electrostatic landscape in graphene. For the thicker hBN layer (device 3), the effect of the local gating field is reduced even further for both positive and negative bias. Hence, no observable influence of the movement of oxygen vacancies is discernible in the Dirac peaks for this device. The graphene channel is found to be electron doped.

To understand the evolution of the antihysterepan class="Chemical">sis with the sweep range of the externpan class="Chemical">al gate bias and temperature, we analyze the charge transport behavior for devices 2 and 3, as shown in Figure , with the same sweeping rate of 1.7 s/V as device 1. When STO is cooled beyond the phase-transition temperature of 105 K, the positively charged oxygen vacancies accumulate at the domain wall network. We note that the sweeping electric field and its polarity has different effects on the movement of the oxygen vacancies. At positive gate voltages, the screening by positively charged oxygen vacancies is less effective, resulting in a higher effective field acting on graphene. The oxygen vacancies thus move away from the STO surface and are manifested in a positive shift of the CNP in graphene. A negative gate voltage, however, enhances the screening efficiency and thereby leads to a lower effective field acting on graphene. In this case, the oxygen vacancies move to the STO surface, resulting in a negative shift of the CNP. The distribution of the mobile oxygen vacancies across the STO surface also depends on the choice of the sweeping range. With increasing sweep range of the gate voltage, the movement of oxygen vacancies enhances, leading to a larger antihysteresis. When the temperature is increased, the sharp decrease in the dielectric permittivity and the disorderliness of the domain wall network in STO result in a lower effective electric field, causing a reduction in the movement of the oxygen vacancies and a reduction in the antihysteresis of the sheet resistance in graphene, which is manifested in Figures b,c and 2b,c,e,f.

The evolution of the antihysterepan class="Chemical">sis with the sweep rate is pan class="Chemical">also investigated, as shown in Figure . The sweep rate does not affect the antihysteresis curves when graphene is directly on STO (Figure a), whereas with the hBN intermediate layer, the antihysteresis decreases with an increasing sweep rate (Figure b). This can be explained in light of the movement of the oxygen vacancies when gating with STO. The movement of the oxygen vacancies is much faster than the fastest sweep rate, which we can achieve with our current setup, and thus the change of the sweep rate does not affect the antihysteresis for device 1. However, with hBN as an intermediate layer, we observed a different scenario: the CNPs overlap for the retracing, while the negative shift of the CNP is less for tracing with a faster sweep rate. Although the movement of the oxygen vacancies is still much faster than the sweeping rate, the screening effect from the oxygen vacancies at the negative gate voltage takes time to be manifested in the square resistance of graphene when hBN is the intermediate layer. Thus, the negative shift is less with a faster sweep rate.

Figure 3

Sweeping rate dependence of the antihysteresis in graphene square resistance for devices 1 and 2. (a) Antihysteresis curves for device 1 do not depend on the sweeping rate. (b) Antihysteresis for device 2 decreases with increasing sweeping rates.

Sweeping rate dependence of the antihysteresis in n>an class="Chemical">graphene square resistance for devices 1 and 2. (a) Antihysteresis curves for device 1 do not depend on the sweeping rate. (b) Antihysteresis for device 2 decreases with increasing sweeping rates.

To further elucidate the dynamics of the pan class="Chemical">oxygen vacancies, their time dependence responses were investigated. Figure a,b shows the time dependence of the pan class="Chemical">graphene square resistance for device 1 at 4 K at Vg = 0 V and Vg = −80 V, respectively. For Vg = 0 V, the time dependence for both trace and retrace cycles is monitored. During the forward trace (from −80 to 0 V), graphene is populated with electrons at Vg = 0 V (t = 0 s), and the square resistance is stable over time, while during the retrace, graphene is populated with holes at Vg = 0 V (t = 0 s). The time dependence, interestingly, shows that the charge carrier type changes from holes to electrons after 190 s and is explained by the movement of oxygen vacancies. When the back gate is swept from a positive gate voltage to 0 V and held at 0 V, the static charges respond immediately, while the mobile oxygen vacancies exhibit a delay in response. Thus, while at t = 0 s, more holes are induced in graphene; over time, an enhancement in the movement of oxygen vacancies changes the charge carrier type to electrons.

Figure 4

Temporal behavior of the graphene square resistance at static gate voltages. (a) Time dependence of the graphene square resistance at Vg = 0 V at 4 K for device 1 (graphene/STO). The red curve shows the time-dependent graphene square resistance after the gate voltage is swept from Vg = −80 to 0 V (tracing), while the black curve shows the time dependence after the gate is swept back to 0 V after 80 V (retracing). (b) Time dependence of the graphene square resistance at Vg = −80 V at 4 K for device 1 (graphene/STO). (c) Time dependence of the CNP positions at different temperatures at −80 V for device 1 (graphene/STO). (d) Relaxation time versus temperature for device 1 (graphene/STO: black dots) and device 2 (graphene/8 nm hBN/STO: red dots). From these curves, the activation energy can be extracted. The activation energy for both devices 1 and 2 is similar. (e) Schematic cross sections of the evolution of electrostatic charge distribution when holding the gate voltage at −80 V. (I)–(III) are defined in (b).

Temporal behavior of the n>an class="Chemical">graphene square resistance at static gate voltages. (a) Time dependence of the graphene square resistance at Vg = 0 V at 4 K for device 1 (graphene/STO). The red curve shows the time-dependent graphene square resistance after the gate voltage is swept from Vg = −80 to 0 V (tracing), while the black curve shows the time dependence after the gate is swept back to 0 V after 80 V (retracing). (b) Time dependence of the graphene square resistance at Vg = −80 V at 4 K for device 1 (graphene/STO). (c) Time dependence of the CNP positions at different temperatures at −80 V for device 1 (graphene/STO). (d) Relaxation time versus temperature for device 1 (graphene/STO: black dots) and device 2 (graphene/8 nm hBN/STO: red dots). From these curves, the activation energy can be extracted. The activation energy for both devices 1 and 2 is similar. (e) Schematic cross sections of the evolution of electrostatic charge distribution when holding the gate voltage at −80 V. (I)–(III) are defined in (b).

We further find that for Vg = −80 V, the carrier type also changes from hole to electron after apn>proximately 330 s, which is pan class="Chemical">confirmed by the gate sweep afterward (as shown in Figure S4). While holding the gate bias at −80 V for 330 s, we find that screening effects due to the movement of oxygen vacancies gradually increase and stabilize over time due to a transition of carrier type to electrons from holes. This process is illustrated in Figure e. At a positive gate voltage, as discussed earlier, the screening is less effective to induce a change in carrier type in graphene (as discussed in SI Figure S5). To elucidate this, we look at the calculations of oxygen vacancy formation energies close to the STO surface layers. The formation energy of oxygen vacancies at layer n is defined by , where Evac,, EO, and E0 are the total energies of the STO slab with vacancies at layer n, oxygen molecule, and pristine STO slab, respectively. The layer-resolved evolution of the formation energies from the surfaces of STO to the bulk are shown for the TiO2 termination of the STO substrate in Figure c. One can easily see that oxygen vacancies have much lower formation energy at the surface TiO2 layer than the subsurface SrO layer. The evolution of the total energies from the surface to the bulk region for both TiO2 and SrO atomic layers are computed using the nudged elastic band (NEB) method[38,39] and are also shown in Figure c.

Figure 5

(a) Top view of STO and the tangent plane in the black rectangle are shown in (b). From the surface layer to the bottom layer, the layers are labeled as 1–10. A schematic representation of the oxygen vacancy diffusion path is shown in (b), and its energy profile obtained from NEB calculations (solid green line) is shown in (c) along with vacancy formation energy in each layer for two types of surface termination (dashed blue and red lines are for TiO2 and SrO layers, respectively). Etot is the total energy of each NEB image, and that of the first image is set to zero.

(a) Top view of STO and the tangent plane in the black rectangle are shown in (b). From the surface layer to the bottom layer, the layers are labeled as 1–10. A schematic repn>resentation of the n>an class="Chemical">oxygen vacancy diffusion path is shown in (b), and its energy profile obtained from NEB calculations (solid green line) is shown in (c) along with vacancy formation energy in each layer for two types of surface termination (dashed blue and red lines are for TiO2 and SrO layers, respectively). Etot is the total energy of each NEB image, and that of the first image is set to zero.

We further analyzed the time dependence of the square repan class="Chemical">sistance at different temperatures. The redistribution of oxygen vacancies at the STO surface is reflected as the change in the CNP in graphene. We can translate the time-dependent square resistance in graphene into shifts of the positions in the CNP, as shown in Figure c, from which the relaxation time can be deduced (Figure d). From this, we evaluated the diffusion energy of the oxygen vacancies. From the Arrhenius equation, , the activation energy has been determined to be 15 meV, as shown in Figure d. The same analysis is performed for both graphene/STO and graphene/hBN (8 nm)/STO. Both show similar values of the activation energies, which again indicates that the mechanism for the antihysteresis is the movement of the oxygen vacancies and is not an interface-driven effect. We further studied this using DFT calculations. As shown in Figure a,b, the diffusion process of oxygen vacancies from the bulk to the surface in STO can be divided into two steps: the first step involves the oxygen vacancy movement from TiO2 to the SrO layer, and in the second step, from SrO to the TiO2 layer. The energy profiles of the two steps shown in Figure c show that the energy barrier is about 50 meV for one oxygen vacancy moving into bulk STO. A more interesting behavior is observed near the surface region in STO. There is a relatively large energy barrier from layer 3 (TiO2) to layer 2 (SrO), impeding the movement of oxygen vacancies, thus effectively confining them up to the third layer. At the same time, the formation energy of the oxygen vacancy is quite low at the surface TiO2 layer (layer 1), indicating the ease of formation of these vacancies. These will be trapped in this layer and face a very strong activation barrier to go to the bulk.

Furthermore, from DFT cpan class="Chemical">alculations and Bader charge anpan class="Chemical">alysis,[40] as discussed in the SI, of fully geometry-optimized structures, we demonstrate that in a graphene/STO system without an external electric field, graphene receives electrons from STO 3.6 × 1013 e/cm2 for a SrO-terminated surface, and loses 1.5 × 1012 e/cm2 for a TiO2-terminated surface. Moreover, according to the density of states (DOS) calculations (Fig. S9 in the SI), graphene receives more electrons when STO has oxygen vacancies. The finding that local differences in the surface-terminating planes in STO, due to the presence of both SrO and TiO2 layers, can bring about differences in the local electrostatic landscape in the graphene channel, is important and also explains the finite differences that are observed in the features in the Dirac peak while using different contacting electrodes in the graphene channel (as discussed in SI Figure S3).

Conclusions

In summary, we have used a monolayer pan class="Chemical">graphene flake to unveil the movement of pan class="Chemical">oxygen vacancies and their interaction with the ferroelastic domain walls in the vicinity of an STO surface. The unique physical and chemical properties in graphene, specifically the van der Waals bonding nature, prevent unwanted electrode–substrate (STO) interaction, facilitating the identification of the exact mechanism responsible for the antihysteresis in the square resistance in graphene. The tunability with an electric field and temperature, of the features observed in the Dirac peak, establishes the interplay of oxygen vacancies and domain walls in STO and rules out the role of any interfacial effects, unlike that reported earlier. We demonstrate that the antihysteresis is caused by the delayed movement of the oxygen vacancies with an estimated energy barrier of 50 meV. Our experimental findings are further corroborated by studies of the formation energies, energy barriers for diffusion of the oxygen vacancy toward the surface of STO, and charge transfer at the graphene/STO interface, calculated from DFT studies. Our approach can be extended to the study of the dynamics of electronic and ionic charge transport and their retention characteristics, using two-dimensional (2D) materials onto oxide substrates, paving the way for electric field control of memory functionalities in electronic devices.

Methods Section

Device Fabrication

The one side polished pan class="Chemical">STO(001) substrates (Crystec GmbH) were treated with a standard protocol with the intention to attain surface-terminating planes of TiO2, and the AFM image of the TiO2-terminated surface is shown in Figure S1. Graphene flakes and hBN flakes (from HQ graphene) were mechanically exfoliated on 300 and 90 nm SiO2/Si substrates, respectively. The thickness of the graphene flake was selected based on the optical contrast and verified by Raman spectroscopy (Figure S7), while the thickness of the hBN flake was determined by atomic force microscopy (AFM) (details in SI, Figure S1). The flakes were then transferred onto STO substrates using the dry-transfer technique,[22] resulting in graphene/STO stack (device 1), graphene/hBN (8 nm)/STO stack (device 2), and graphene/hBN (23 nm)/STO stack (device 3), respectively. To reduce polycarbonate residues on the flakes from the transfer steps, we annealed the devices in an Ar/H2 atmosphere for several hours at 200 °C. Standard electron-beam lithography and electron-beam deposition were followed to make the contact electrodes consisting of Al (5 nm)/Co (35 nm)/AlO (1.8 nm), where AlO was achieved by first depositing Al and in situ oxidization. The final devices are shown in Figure and the optical images of the devices are shown in Figure S2.

Electric Measurements

The graphene channel resistance was monitored using low-frequency (3–15 Hz) lock-in (Stanford Research Systems SR830 lock-in amplifier) technique, and the gate voltage was sourced by a Keithley 2410 Source Meter. The sample was put in a Microstat He2 flow cryostat (Oxford Instruments), and temperature varied between 4 to 300 K for our studies.

Theoretical Methods

We have performed density functionpan class="Chemical">al calculations to estimate the surface polarization and the formation energies of oxygen vacancies in different layers in STO. All of the calculations are performed using the projector augmented wave (PAW) method[41,42]-based density functional code Vienna Ab initio simulation package (VASP).[43,44] The exchange–correlation potential was treated within the generalized gradient approximation in the form proposed by Perdew, Burke, and Ernzerhof (PBE).[45] The wave function was expanded in a plane wave basis with an energy cutoff of 500 eV, and a 5 × 5 × 1 γ-centered k-point sampling was used. Perpendicular to the surface, a vacuum of more than 15 Å was used. To avoid the long-range electrostatic interactions between layers, we have used dipole correction in all of the calculations. The structures were optimized using the conjugate gradient (CG) and RMM-DIIS quasi-Newton algorithms[46] until the Hellman–Feynman force on each atom was less than 0.01 eV/Å. For Bader charge analysis, the wave functions were expanded on a 70 × 70 × 196 grid, and the localized charges were mapped on a 140 × 140 × 392 grid. For the charge partitioning, we employed the algorithm proposed by Henkelman et al.[47] Postprocessing of some calculations were performed using VASPKIT[48] and VESTA.[49]