Controlled manipulation of oxygen vacancies using nanoscale flexoelectricity

Oxygen vacancies, especially their distribution, are directly coupled to the electromagnetic properties of oxides and related emergent functionalities that have implications for device applications. Here using a homoepitaxial strontium titanate thin film, we demonstrate a controlled manipulation of the oxygen vacancy distribution using the mechanical force from a scanning probe microscope tip. By combining Kelvin probe force microscopy imaging and phase-field simulations, we show that oxygen vacancies can move under a stress-gradient-induced depolarisation field. When tailored, this nanoscale flexoelectric effect enables a controlled spatial modulation. In motion, the scanning probe tip thereby deterministically reconfigures the spatial distribution of vacancies. The ability to locally manipulate oxygen vacancies on-demand provides a tool for the exploration of mesoscale quantum phenomena and engineering multifunctional oxide devices.

O xygen vacancies (V ÁÁ o ) are elemental point defects in oxides, and they generally function as mobile electron donors. At sufficiently high concentrations, they can disturb the ground state and promote emergent functional phenomena, such as superconductivity 1 , ferromagnetism 2, 3 , metalto-insulator transition 4 and interface conductivity 5 . In addition, the distribution and dynamics of V ÁÁ o play a central role in many oxide-based energy and memory applications 6,7 . Hence, the ability to manipulate V ÁÁ o provides an opportunity to study the evolution of emergent phenomena and control numerous functionalities, which are essential for developing next generation oxide devices 8, 9 . Traditionally, modification of the vacancy concentration has been carried out by annealing oxides at high temperature in reducing/oxidising atmosphere. Recent works have shown that at room temperature, the vacancy concentration can instead be locally changed by employing an electrical bias from a scanning probe microscope (SPM) tip 4,10 . Such local manipulation enables reversible nanoscale control of metal-insulator transitions 4 and the modulation of interface conductivity 11 . However, in addition to the change in the vacancy concentration, the application of bias through an SPM tip is often accompanied by charge injection 12,13 and the formation of protons/ hydroxyls 14 , which complicates the practical implementation of this approach.
Recently, it has been shown that the mechanical force from the SPM tip can also deplete V ÁÁ o from the contact region 15,16 . However, because this mechanical depletion is less understood, the feasibility of using the force from an SPM tip as an active means of manipulating V ÁÁ o has not been explored. An explanation for the mechanical depletion of V ÁÁ o has been proposed through the so-called piezochemical coupling mechanism [15][16][17] , which involves the converse Vegard effect and the flexoelectric effect. The converse Vegard effect accounts for the decrease in vacancy concentration due to force-induced lattice compression. Meanwhile, the flexoelectric effect considers the electromigration of positively charged V ÁÁ o under the stress-gradient-induced flexoelectric field. Notably, by definition, the flexoelectric effect refers to the generation of an electrical polarisation by stressgradients 18,19 . In this context, Yudin and Taganstev suggested that the flexoelectric field is not a macroscopic electric field 19 ; instead, it should be understood as a pseudo-internal field that polarises a medium in the presence of a stress-gradient [20][21][22][23] . Hence, it remains unclear how this flexoelectric field acts on V ÁÁ o . To be able to mechanically manipulate V ÁÁ o , we must, therefore, thoroughly understand how they respond to the mechanical force from the SPM tip. Subsequently, we have to devise a strategy to employ this force to increase as well as decrease the vacancy concentration in a controlled manner.
Here, we demonstrate a controlled manipulation of V ÁÁ o using the mechanical force from an SPM tip. As a model system, we used a homoepitaxial thin film of SrTiO 3 , which is an archetypal quantum paraelectric oxide with well-known flexoelectric coefficients 24,25 . Using a Kelvin probe force microscopy (KPFM)based imaging scheme, we show that besides pushing V ÁÁ o away in the vertical direction, the force promotes the lateral transport of vacancies during the movement of the tip. Using phase-field simulations, we argue that this mechanical redistribution of V ÁÁ o is driven by the depolarisation field associated with the stressgradient-induced flexoelectric polarisation. The depolarisation field underneath the tip promotes the vertical migration of vacancies, whereas around the contact edge, it traps V ÁÁ o that can move with the tip. Furthermore, we demonstrate that altering the tip geometry can tailor this depolarisation field to preferentially allow the lateral transport of V ÁÁ o . This approach enables a controlled spatial modulation of vacancies.

Results
Probing oxygen vacancies with Kelvin probe force microscopy. We start by discussing the concept of probing distribution of V ÁÁ o with the KPFM technique ( Fig. 1a-c). KPFM is a surface sensitive technique that measures the contact potential difference (CPD) between the tip and SrTiO 3 (STO) surface. The CPD is defined as the offset in the respective vacuum energy levels (in units of V) 26 . However, as schematically illustrated in Fig. 1b, this CPD would change if the STO surface contains V ÁÁ o , which can be locally accumulated by an electrical poling. This change in the CPD can originate from the workfunction of STO 27 , chemical dipoles (including V ÁÁ o -electron pairs), and their orientation 28 . Moreover, the tip bias used during the KPFM measurement has been argued to influence the measured CPD 29 . The contribution of each factor cannot be individually separated, which inhibits a quantitative determination of the vacancy concentration with the KPFM technique. Nonetheless, as shown in Fig. 1c, the KPFM contrast across this V ÁÁ o -rich region can be argued to scale proportionally with the vacancy concentration 27 . To establish a proof of concept, in the following we elaborate on the application of KPFM technique to study the diffusion of V ÁÁ o using a 14 and 120-unit cell (uc)-thick STO film (details of the film are provided in the Methods section and in Supplementary Fig. 1). Figure 1d, e show KPFM images after poling the pristine surface of STO films with a tip bias of −5 V. The images of the 14-uc (120-uc)-thick STO film were taken 20 (16) minutes and 380 (360) minutes after poling. An image contrast is clearly visible across the poled region, which is appearing as a dark rectangular patch. This implies surface charging, either by injected charges, protons/hydroxyls, and/or V ÁÁ o 30, 31 . The injected charges and protons/hydroxyls are expected to decay with a timescale, which is independent of the STO thickness 31 . In contrast, because of diffusion or surface reaction that enables the recombination of V ÁÁ o with oxygen from the ambient, the vacancy concentration is expected to decrease with a pronounced thickness-dependent timescale 32 . Figure 1d, e indeed show that the contrast across the poled region diminishes over time, and the decrease is largest in the 14-uc-thick STO film. Therefore, we argue that the surface charging is caused by V ÁÁ o , which either undergo diffusion or recombine with oxygen from the ambient.
We can distinguish the dominating mechanism that causes the vacancy concentration to decay over time by analysing the timedependency of the degree of equilibrium, S(t) that can be calculated from the KPFM images (see Supplementary Note 1). The time-dependency of S(t) describes how the surface equilibrates after the electrical poling. S(t) will exhibit a semiparabolic (linear) time-dependency if diffusion (surface reaction) is the dominating mechanism 27, 32 . Figure 1f plots S(t) as a function of time. Evidently, the time-dependencies of S(t) are semi-parabolic for both films-implying that the diffusion of V ÁÁ o causes the vacancy concentration to decrease over time.
Arguably, during the poling, the applied tip bias perturbs the equilibrium V ÁÁ o -distribution of the entire film including the surface, which equilibrates through the diffusion of V ÁÁ o . This diffusion occurs both along the out-of-plane and in-plane directions. However, due to the high surface sensitivity of the KPFM technique, the surface-bulk diffusion of V ÁÁ o predominantly affects the time evolution of the KPFM contrast, and thus the time-dependency of S(t). Notably, during this surface-bulk diffusion, the repulsive vacancy-vacancy interaction inhibit V ÁÁ o to migrate independently along the out-of-plane direction 33 . Thus, the time-dependency of S(t) effectively describes how the perturbed volume under the poled area (in Fig. 1d, e) equilibrates (see Supplementary Note 1 for a detailed discussion). Naturally, this volume would be smaller in the thinner STO film, and thus would equilibrate faster. This explains why the KPFM contrast (S(t)) diminishes (grows) more rapidly in the 14-uc thick STO than in the 120-uc thick STO film (Fig. 1d-f).
Following the rationale above, we fit the time evolution of S(t) with Fick's 2nd law of diffusion (solid lines in Fig. 1f). From this fitting, we obtained the diffusion coefficient D = 9.4(3) × 10 −19 and 3.8(2) × 10 −18 cm 2 s −1 for the 14-uc and 120-uc-thick STO film, respectively. These values are well in the range of the bulk value D bulk ≈ 10 −(17±3) cm 2 s −1 (300 K-extrapolated) 34 , which validates the conceptual schematic depicted in Fig. 1c. Subsequently, we utilised this correlation between the vacancy concentration and KPFM signal to image the vacancy redistribution under an applied bias and force.
Oxygen vacancy redistribution under applied bias and force. In this section using the 120-uc-thick STO film, we compare the response of V ÁÁ o to the electrical bias and force from the SPM tip. For this study, we employed a sharp tip with a radius of curvature = 25 nm. As depicted in the KPFM image (Fig. 2a), we scanned the left side of a formerly V ÁÁ o -enriched region with an increasingly positive tip bias (0-5 V) at a contact force = 0.6 µN. Meanwhile, the right side was scanned with an increasing contact force (0.6-8.5 µN) by holding the tip at the ground potential. In either case, the tip crossed the border between the formerly V ÁÁ o -enriched and pristine regions. The central part of the KPFM image is uniformly dark, which indicates a uniform vacancy distribution. However, a change in the image contrast, which stems from the V ÁÁ o -redistribution is visible within both electrically and mechanically scanned areas. Notably, these scans together with acquiring the KPFM image took about 20 min. Because of their ultra-slow decay as evident from Fig. 1e, during this time V ÁÁ o can be assumed to be kinetically frozen. This implies that the V ÁÁ o -redistribution is due to the applied bias and force.
To quantify this V ÁÁ o -redistribution, we constructed the vacancy concentration map from Fig. 2a by defining the normalised vacancy concentration, Here, V b and V min represent the baseline and minimum KPFM signals, which were extracted from areas indicated by the black (top-right corner) and white (centre) squares, respectively. The reconstructed NVC map is shown in Fig. 2b. The positive (negative) NVC refers to a higher (lower) vacacny concentraion compared to the intrinsic surface concentration of V ÁÁ o (NVC = 0, at the top-right corner). Profiling the NVC map along lines E and M indicates that both the applied bias and force depleted the formerly V ÁÁ o -enriched region (Fig. 2c). We find that this depletion can be modelled with a phenomenological Boltzmann sigmoid function in the following 14 uc 14 uc where A 1 , A 2 , X o , and M refer to the initial value, final value, center, and decay rate, respectively. The solid lines in Fig. 2c are fits to this function. Table 1 lists the best-fit parameters. This functional analysis suggests that the applied positive bias and force have a same qualitative effect along the lines E and M, respectively. By comparing the X o ( = 0.5A 1 ) values from the fits, we obtain a force-voltage equivalence factor of 0.4 V µN −1 . The Boltzmann sigmoid fits in Fig. 2c also exhibit an onset. Defining this onset as 0.9A 1 (indicated by stars), on an ad hoc basis, yields a threshold voltage (V th ) and force (F th ) of 1.3 V and 2.7 µN or equivalently 1.1 V, respectively. These threshold values are consistent with the activation barrier potential ( = 1−1.4 V) for the electromigration of V ÁÁ o in STO 35,36 , which further validates the force-voltage equivalence relationship.
Surprisingly, electrical and mechanical scans yield the opposite image contrast across borders in the NVC map. As shown in Fig. 2d, along line E L , NVC monotonically decreases with the increasing bias. However, along line M L, NVC varies nonmonotonously: first, it increases as the force increases and then gradually drops; the horizontal line in Fig. 2d marks the background (NVC = 0.1 at 0 V and 0.6 µN). We, therefore, conclude that the electrical scan depletes the pristine surface outside the left border, while the mechanical scan enriches the pristine surface outside the right border.
The force-induced vacancy-enrichment of the pristine region implies that some of the depleted V ÁÁ o moved laterally with the tip across the border during the mechanical scan. Additional experiments elaborate that the mechanical scan does not alter the background within the pristine region, and thus rule out the formation of V ÁÁ o and triboelectric charging of the STO surface during the mechanical scan (see Supplementary Note 2). To quantify the lateral motion of V ÁÁ o , in Fig. 2e we, therefore, show the NVC with the background (NVC at 0.6 µN) subtracted (i.e., ΔNVC) along lines M and M L . Because KPFM is a surfacesensitive technique, the remnant V ÁÁ o on the surface contribute to Table 1 Best-fit parameters from the Boltzmann sigmoid fit to the normalised vacancy concentration (NVC) data in Fig. 2c Overall, our main observations are as follows. First, the contact force is bifunctional: it depletes the formerly V ÁÁ o -enriched region and simultaneously enriches the pristine surface. Second, the mechanical V ÁÁ o -redistribution involves a predominant surfacebulk migration and a relatively weaker lateral motion of V ÁÁ o with the tip.

Modelling the mechanical redistribution of oxygen vacancies.
To understand the mechanical redistribution of V ÁÁ o we first considered two mechanisms-the converse Vegard effect and the flexoelectric effect 15,16 . Recently, the magnitude of these two effects under an applied force from SPM tip has been compared in a PbTiO 3 (PTO) thin film 38 . This study suggests that the converse Vegard effect is much weaker than the flexoelectric effect. Notably, both the PTO and STO have comparable flexoelectric and Vegard coefficients 25,38,39 , which determine the relative contributions of these two effects for a given force. Based on these considerations, we thus conclude that the flexoelectric effect predominantly causes the mechanical redistribution of V ÁÁ o , and the contribution from the converse Vegard effect is marginal.
For gaining a mechanistic understanding of this flexoelectric effect-driven V ÁÁ o -redistribution, we performed phase-field simulations; whereby we incorporated the flexoelectric effect 40 and coupled the time-dependent Ginzburg-Landau and the Nernst-Planck equations 41,42 . The simulation was performed assuming that the STO is paraelectric (see Supplementary Note 3). Unlike in the experiment, the SPM tip was assumed to be static, and following the Hertzian model the contact radius was calculated to be 8 nm for a contact force of 4 µN. Initially, V ÁÁ o were assumed to be homogeneously distributed over the entire STO thickness, instead of being localised on the surface. Despite these oversimplified assumptions, our simulation still provides a qualitative insight into how V ÁÁ o respond to mechanical stimuli. A detailed explanation of our model and a discussion on the redistribution of V ÁÁ o under a positive tip bias are included in Supplementary Notes 4, 5.
The stress-gradient from the SPM tip locally polarises STO through the flexoelectric effect 24 , as evident from the out-of-plane polarisation vector map in Fig. 3a. Since in the simulation the STO film is assumed to be in the paraelectric phase, the polarisation in Fig. 3a purely stems from the flexoelectric effect. This flexoelectric polarisation reaches a maximum (~0.04 C m −2 ) underneath the tip and spatially varies both in magnitude and direction. The resulting polarisation bound charge and the associated depolarisation field accordingly redistribute vacancies around the tip-STO contact region. This redistribution can be visualised from the simulated NVC map in Fig. 3b, which plots the in-plane distribution of V ÁÁ o at the surface. Clearly, the vacancy concentration is decreased (enhanced) underneath the tip (around the contact edge).
An intuitive understanding of the simulated V ÁÁ o -redistribution can be gained from the component-resolved distribution of the depolarisation field. As shown in Fig. 3c, the z-component, E dep z , points downward (upward) below the tip (around the contact edge), whereas the x-component, E dep x , exhibits a parentheses-like structure: a node at the contact point and antinodes around the contact edge (Fig. 3d). Effectively, E dep x points inwards, as indicated by the white arrows in Fig. 3d. The y-component (not shown), E dep y , forms an analogous structure to E dep  Supplementary Figs. 9-11). These experiments also yielded a weaker lateral motion but a stronger surface-bulk migration of V ÁÁ o -highlighting the dominating influence of depolarisation field underneath the tip. In the following, we illustrate that the depolarisation field around the tip-STO contact junction can be tailored in favour of the lateral motion of V ÁÁ o , which enables controllably manipulating the vacancy distribution.
Controlled spatial modulation of oxygen vacancies. The basic concept of tailoring the SPM tip-induced depolarisation field can be understood with the aid of Fig. 4a, which compares the simulated surface deformation profiles under a static load of 4 µN using two different tip geometries. The upper panel corresponds to the spherical tip (contact radius = 8 nm) that is used in Fig. 3, and the lower panel corresponds to a flat-ended tip (contact radius = 15 nm). Compared to the spherical one, the flat-ended tip usually imparts a weaker stress on the STO surface underneath. Since the downward E dep z scales with the stress-gradient underneath the tip, the flat-ended tip induces a very small downward E dep z (Fig. 4b, upper panel). In contrast, the lateral deformation (indicated by curved arrows in Fig. 4a), which controls the depolarisation field around the contact edge, is alike for both geometries. Consequently, the flat-ended tip induces a depolarisation field distribution around the contact edge (Fig. 4b) similar to that of the spherical tip in Fig. 3c, d. Additionally, an enhanced contact radius enlarges its spatial extent. A selective suppression of E dep z underneath a flat-ended like tip should significantly reduce the surface-bulk V ÁÁ o migration. Meanwhile, the extended depolarisation field around the contact edge should improve the lateral transport of V ÁÁ o . To validate our proposition, we performed experiments with a sharp and blunt tip. The estimated radius of curvature of this blunt tip is larger than 200 nm (see Supplementary Fig. 7). Thus, it effectively yields a flat contact junction underneath the tip. We used the 120-uc thick STO film in these experiments. Figure 4c, d show the NVC maps after mechanical scans were performed with these tips at a contact force of 9.5 µN. Notably, we scanned both the left and right boundaries between the V ÁÁ o -enriched and pristine regions.
To compare the sharp and blunt tip-induced V ÁÁ o -redistribution we profiled the NVC maps, as indicated by vertical lines in Fig. 4c, d. Figure 4e, f show the corresponding NVC profiles. The overlapping NVC profiles along lines M1/M4 and M2/M3 demonstrate that the V ÁÁ o -redistribution is reproducible for both tip geometries. However, the response of V ÁÁ o to the applied force from these two tips are clearly different. The NVC profiles in Fig. 4e exhibit a maximum drop in NVC of Δ dec max ¼ À0:75 along lines M2 and M3 but no appreciable increase in NVC (Δ inc max ) along lines M1 and M4. This implies that the sharp tip strongly depletes the V ÁÁ o -enriched regions but barely enriches the pristine regions. This result is in qualitative agreement with that in Fig. 2e, which shows that the fraction of the depleted V ÁÁ o that laterally move with the tip progressively decreases for applied forces larger than 6 µN. The NVC profiles in Fig. 4f

Discussion
Through a combined experimental and theoretical approach, we demonstrated the flexoelectricity-mediated controlled manipulation of oxygen vacancies by the mechanical force from an SPM tip. A deterministic reconfiguration of spatial vacancy profile provides control over the electron density and related electronic correlation effects. This could enable, using an SPM-based all-inone platform, the investigation of mesoscale quantum phenomena in oxides 43 . Ultimately, this creates the opportunity for developing mechanically sketched oxide devices, and ambipolar mechanical control of device functionalities such as electroresistance states. The voltage-free operation of the SPM tip would thereby eliminate the possibility of surface charging. At this point, we want to emphasise that flexoelectricity is a universal phenomenon, which can occur in any dielectric 18,44 . However, the flexoelectric coefficients of few oxides, such as SrTiO 3 and BaTiO 3 , are currently known 45 . Thus, our work should motivate the study of flexoelectricity in other oxides.
Broadly speaking, our KPFM-based imaging approach offers a time-efficient way of characterising the activation barrier potential for oxygen vacancy migration at room temperature to complement the conventional Arrhenius analysis 46 . Combined with the feasibility of determining the diffusion coefficient, this technique could thus become an essential metrology tool for oxidebased energy and memory research. Furthermore, our theoretical model that couples the phase-field simulations to the Nernst-Planck equation, can be employed to elucidate how depolarisation fields cause oxygen vacancies., electrons, and holes to redistribute. Therefore, the model can be extended to the study of emergent problems such as the domain wall conductivity, high electrical conductivity of morphotropic phase boundaries, and leakage current in ferroelectric oxides [47][48][49][50] .

Methods
Thin film growth. SrTiO 3 thin films were homoepitaxially grown on TiO 2 -terminated Nb:SrTiO 3 (0.5% wt. doped) substrates using pulsed laser deposition technique. The growth dynamics and thickness were monitored by in-situ reflection high energy electron diffraction technique. The depositions were performed at 1000°C and using an oxygen partial pressure of 5 × 10 −7 torr. After deposition, films were annealed at 800°C for an hour in a 1 torr oxygen atmosphere and subsequently cooled down to room temperature at a cooling rate of 20°C per minute.
Kelvin probe force microscopy. KPFM measurements were carried out using the Asylum Research Cypher SPM at room temperature and under ambient conditions. Pt/Ir-coated metallic tips (NANOSENSORS™ PPP-NCHPt) with a nominal spring constant ≈ 40 N/m were used for electrical/mechanical scans and KPFM imaging. The KPFM measurements were obtained in the non-contact mode using a lift height of 30 nm and the typical scan parameters used are as follows: V ac = 1 V (peak-to-peak), f resonance = 250 kHz, and scan rate = 1 Hz. Before each experiment, the spring constant of the cantilever was accurately determined from force-distance measurements and thermal tuning methods. The contact force during mechanical scans was varied accordingly by controlling the set-point voltage.
Theoretical modelling. To model the oxygen vacancy redistribution by mechanical force, we performed phase-field simulations by coupling the time-dependent Ginzburg-Landau and Nernst-Planck equations.
In Eq. (2), P is the polarisation vector, L is the kinetic coefficient and F is the total free energy of the system, which includes Landau, electric, gradient and flexoelectric energy contributions. In Eq. To solve Eqs. (2) and (3), the system was discretised into 100Δx × 50Δy × 500Δz grid points to implement the semi-implicit Fourier method 45 . The 120-ucthick STO film is simulated to be 480Δz in thickness, while the substrate and air are each 10Δz in thickness. The parameters for the total free energy of STO were adopted from work of Y.L. Li et al. 46 . Following the work of R. Moos et al. 47 , the initial concentration of oxygen vacancies in STO is assumed to be 3.66 × 10 14 cm −3 based on our thin film growth conditions. In addition, the diffusion coefficient is assumed to be constant with regard to pressure and calculated as 1.23 × 10 −15 cm 2 s −1 at room temperature 47 . A sufficiently long simulation time is used to ensure that the induced polarisation and oxygen vacancy concentration reach a quasi-steady state. At each time step, the electrostatic and elastostatic equilibrium equations are solved under the electric short-circuit 48 and mechanical mixed boundary conditions 49 , respectively.
Data availability. The data that support the findings of this study are available from the corresponding authors upon reasonable request.