Oxygen vacancy-induced topological nanodomains in ultrathin ferroelectric films

Oxygen vacancy in oxide ferroelectrics can be strongly coupled to the polar order via local strain and electric fields, thus holding the capability of producing and stabilizing exotic polarization patterns. However, despite intense theoretical studies, an explicit microscopic picture to correlate the polarization pattern and the distribution of oxygen vacancies remains absent in experiments. Here we show that in a high-quality, uniaxial ferroelectric system, i.e., compressively strained BaTiO3 ultrathin films (below 10 nm), nanoscale polarization structures can be created by intentionally introducing oxygen vacancies in the film while maintaining structure integrity (namely no extended lattice defects). Using scanning transmission electron microscopy, we reveal that the nanodomain is composed of swirling electric dipoles in the vicinity of clustered oxygen vacancies. This finding opens a new path toward the creation and understanding of the long-sought topological polar objects such as vortices and skyrmions.


INTRODUCTION
In the past 10 years, an intriguing query in ferroelectricity studies has been the existence of topological arrangement of electric dipoles such as vortices and skyrmions in magnetic systems [1][2][3][4][5][6][7][8][9][10] . These topological objects hold considerable interests from both fundamental and practical aspects. However, unlike spins that are less coupled to the lattice strain and can rotate away from the high-symmetry direction to favor noncollinear, chiral arrangements, the arrangement of electric dipoles is strongly restricted by the crystal symmetry 9 . For example, ferroelectric systems usually exhibit Ising-type domain walls, of which the electric dipoles gradually vanish and then reverse along the polar axis, contrasting the chiral Bloch-or Neel-type domain walls in magnetic systems 11 . This difference renders the pursuit of ferroelectric topological objects particularly challenging. To date, successful cases have been only reported in a limited number of ferroelectric film systems, including PbTiO 3 and BiFeO 3 films, where polar vortices 2,4,5 , skyrmions 3,10 , and merons 12 have been discovered. These were mainly achieved via designing mechanical and electrical boundary conditions to flatten the energy landscape of ground-state polar structures and bring about the rotation of polar vectors. Another rarely investigated pathway lies in the exploitation of ionic point defects, which can break the lattice symmetry locally via electric and strain fields. For example, hedgehog-and vortex-like polarization states were recently discovered in the vicinity of charged point defects in BiFeO 3 thin films 13 . The charged defects were immobile cation impurities intentionally introduced during film growth. Despite some promise, this success is limited to specific materials with complex growth design, thus falling short of generality.
Oxygen vacancies, as a prevalent type of point defects in oxide materials, have been revealed to play an influential role in the properties of oxide ferroelectrics on various aspects 14,15 . For example, it has been well documented that oxygen vacancies can be attracted and accumulated at the domain walls, heterointerfaces, and other lattice defects as pinning sites to cause unwanted polarization fatigue 16,17 but also trigger intriguing emergent phenomena in some cases [18][19][20] . A fundamental aspect to comprehend the role of oxygen vacancies is the microscopic picture of the interplay between oxygen vacancies and the polar order [21][22][23][24][25] . Theoretical studies using first-principles calculations and force fieldbased methods have managed to depict atomic pictures, where oxygen vacancies can cluster and significantly impact local polarization profiles via lattice distortion, e.g., inducing 180°charged domain walls and forcing polarization to tilt away from the highsymmetry axis [22][23][24] . These findings imply that oxygen vacancies can be coupled to exotic polarization patterns with lower symmetries. However, the majority of experimental studies have so far only managed to identify the passive role of oxygen vacancies in their interaction with electric dipoles, e.g., electrostatic passivation on charged domain walls 19,25,26 . The more active role of oxygen vacancies in the creation of polarization domains has not been identified.
Here we show that high-quality, oxygen-deficient BaTiO 3 films (with a thickness below 10 nm), which are compressively strained on SrTiO 3 (001) substrates, can exhibit an exotic, dynamic evolution from a uniform polarization state to exotic nanoscale polarization structures in an aging process under ambient conditions. By piezoelectric force microscopy (PFM) and spectroscopic techniques, we point out the critical role of oxygen vacancy clusters in such a dynamic process. Scanning transmission electron microscopy (STEM) further reveals that oxygen vacancies clusters give rise to polarization vortices without any extended lattice defects. This is in sharp contrast to the common observation in a uniaxial ferroelectric system, e.g., compressively strained BaTiO 3 and PbTiO 3 films with ideal chemistry, of which the tetragonal c axis is strongly favored along the out-of-plane direction, thus exhibiting a pure out-of-plane polarization state 27,28 . Therefore, our finding will inspire future efforts to investigate and exploit the role of oxygen vacancies in the creation of exotic domain structures in ferroelectric materials.

RESULTS
Oxygen pressure-dependent polarization structures by PFM BaTiO 3 (BTO) films with a bottom SrRuO 3 (SRO) electrode layer were grown on TiO 2 -terminated SrTiO 3 (001) substrates by pulsed laser deposition (PLD) 29 . The oxygen partial pressure for BTO growth was varied from 20 to 5 mTorr with the purpose of modulating the density of oxygen vacancies in the film. The high quality of BTO films grown under 5 mTorr (henceforth 5-BTO) was inferred from the diffraction fringes of X-ray diffraction (XRD) and the step-terrace morphology with unit-cell step height by atomic force microscopy (AFM, Fig. 1a, b). By XRD and AFM, we also confirmed the consistent quality of BTO films grown under higher oxygen pressures up to 20 mTorr (Supplementary Figs. 1 and 2).
We first examined the polarization state of the BTO films by PFM. As shown in Fig. 1c, the pristine BTO film grown under 20 mTorr (henceforth 20-BTO) exhibits the same phase contrast as the negatively biased region (where the electric field pointed from the bottom electrode to the top surface of the BTO film), suggesting a pure upward self-polarization (pointing from the bottom SRO electrode toward the top surface). This result is consistent with previous studies on oxygen-rich, ultrathin BTO films 28,30 . Such a preferred polarization orientation was believed to originate from a built-in field generated by surface anionic adsorbates such as hydroxyl (see Supplementary Discussion) 27,31 . By contrast, fresh BTO films grown in lower oxygen pressures, namely 5-BTO, exhibit a uniformly opposite self-polarization in average, namely pointing from the top surface down to the bottom electrode (Fig. 1d). More interestingly, the uniform downward polarization state is metastable and will gradually evolve to a stable "raisin-in-cake" polarization structure after aging in the air or vacuum for over one day (Fig. 1e). As shown in Fig. 2a, the aged polarization structure is composed of nanoregions with dark PFM phase contract (owning a lateral dimension of 100 to below 10 nm), embedded in an upward polarization matrix. Figure 2b shows the cross-section line profiles of the PFM phase and amplitude of two representative regions containing the nanostructures. As the nanostructures exhibit an opposite PFM phase compared to that of the upward polarization matrix, it seems to suggest that these nanostructures possess a downward polarization. However, we also note that the PFM amplitude of these nanostructures is significantly weaker than that of the upward polarization matrix. As PFM estimates an average polarization in the depth direction with a relatively large lateral dimension (tens of nanometers), this weaker amplitude may suggest an inhomogeneous polarization state in the nanostructures (which we will discuss in detail in the STEM results later) 29 . But for convenience, we will use "up domain" and "down domain" to refer to these two polarization structures by PFM in the following discussion, respectively.
The non-uniform polarization structure of the 5-BTO film should point to inhomogeneous distribution of local built-in fields (E in ), which can be revealed by local PFM hysteresis loops (Fig. 2c). The PFM hysteresis loops were measured at different regions, from which we extracted their corresponding coercive fields (E c ) and E in (Fig. 2c). As expected, the "up domain" typically exhibits an asymmetric hysteresis loop shifting to the positive field side, suggesting an E in pointing to the top surface. The corresponding E in can be then estimated from the difference of the positive coercive field (E c + ) and the negative coercive field (E c -), i.e., , and E in of the up domains at 30 different locations, showing an average E in of about 0.3 V. The uniform up-polarization state of the 20-BTO film exhibits similarly asymmetric hysteresis loops with comparable E in (not shown here). Therefore, we believe that the upward polarizations in the two films share the same origin, namely surface anionic adsorbates. In comparison, the hysteresis loop of the "down domain" in the 5-BTO film is oppositely shifted with an average E in of about −0.3 V (Fig. 2d). These results clearly suggest an additional factor in the "down domain," which produces an opposite E in to cancel the field from surface anionic adsorbates. A variety of factors besides surface-adsorbed ions have been accounted for the self-polarization of epitaxial oxide ferroelectrics in previous studies. These include strain gradients 32-34 , interfacial electrostatic potential caused by valence discontinuity 35 or Schottky junction at the electrode/ferroelectric interface 36 , and ionic point defects 31,37 . First of all, since our ultrathin BTO films are fully strained (Fig. 3b), strain gradients resulted from misfit strain relaxation can be ruled out as a possible factor in our study. Second, both BTO and SRO along the [001] direction have alternating stack of charge-neutral AO and BO 2 layers (A is referred to the A site cation and B for B site cation), thus proscribing the existence of interfacial electrostatic potential due to valence discontinuity. In addition, a Schottky junction usually generates an electrostatic potential uniformly in the film plane, thus favoring a homogeneous polarization state. Therefore, it is unlikely that interfacial electrostatic potentials play a vital role in the observed nanodomain structure of the 5-BTO film. Based on these considerations, we are inclined to ascribe the nanodomain structure to ionic point defects, or more specifically, oxygen vacancies. This supposition is consistent with the distinctive domain structures of the two films by just varying the oxygen pressure from 20 to 5 mTorr. A more detailed examination reveals that the nanodomains grow in size and density with the decrease of O 2 partial pressure from 20 to 5 mTorr (Supplementary Fig. 3).
To further testify the critical role of oxygen vacancies, we post-annealed the 5-BTO film at 500°C in O 2 to remove possible oxygen vacancies. As expected, the non-uniform polarization structure disappeared and a pure upward self-polarization was created. Similarly, by post-annealing the 20-BTO film in ultrahigh vacuum, an analogous nanoscale polarization structure can be also built ( Supplementary Fig. 4). These results strongly imply the critical role of oxygen vacancies in the current study.
Then an interesting question is how oxygen vacancies can bring about the dynamic evolution of the domain structure in the 5-BTO film. In the beginning, oxygen vacancies generated during film deposition tend to distribute uniformly in the film at the growth temperature (700°C) as a result of high entropy 38 . Meanwhile, depending on the formation energy of oxygen vacancies in BTO and the oxidizability of the environment in the growth chamber, which may lead to two different cases. (1) Oxygen may escape from the top surface of BTO if the formation energy of oxygen vacancy is too small (that is difficult to be reoxidized). (2) Or considering a slightly oxidizing environment and the diffusivity of oxygen vacancies at such a high temperature, oxygen vacancies may be further attracted to the top surface by their affinity to oxygen on the free surface. In both cases, a density gradient of oxygen vacancies should form along the film thickness direction. Then after a relatively short cooling procedure (~1 h with a cooling rate of 10°C min −1 ), the distribution of oxygen vacancies in the as-prepared 5-BTO film should remain in such a metastable state: namely uniform in the film plane but more populated at the top surface along the thickness direction. The population of oxygen vacancies with positive charges at the top surface of the film can then generate a built-in electric field, pointing downward to the substrate. This built-in field can well explain the as-observed uniform downward polarization in the fresh 5-BTO film 31 . Furthermore, oxygen vacancies have been well known for their aging and clustering behavior at room temperature as a result of large energy gain during clustering 24,38,39 . This aging and clustering process usually occurs with a relatively long relaxation time due to the low diffusion coefficient of oxygen vacancies at room temperature. Consequently, oxygen vacancies gradually evolve into inhomogeneously distributed clusters accompany with locally varying fields, which is consistent with the gradually emerged inhomogeneous polarization configurations observed by PFM.

Examination of oxygen vacancy population
To verify the existence of oxygen vacancy, we measured the photoemission spectra of O 1s core levels by X-ray photoemission spectroscopy (XPS). As displayed in Fig. 3a, three peaks can be observed: namely one main peak centering at~529.7 eV, one satellite peak at~531.8 eV, and another at~532.5 eV. According to the literature, these three peaks have been assigned to the lattice bonded O in BTO, oxygen vacancies, and surface-adsorbed OHand H 2 O, respectively 40,41 . As expected, besides the main peak of lattice bonded O, the 5-BTO sample shows a satellite peak mainly contributed by oxygen vacancies. In comparison, the 20-BTO sample exhibits a satellite shoulder with a major component of surface-adsorbed oxygen and a minor component of oxygen vacancies. This result is consistent with the scenario of surface anionic adsorbates-induced upward self-polarization in the 20-BTO sample, and confirms populated oxygen vacancies in the 5-BTO sample. In addition, we extracted a c-axis lattice constant of 0.424 and~0.416 nm for the 5-BTO and 20-BTO films, respectively, from the reciprocal space mapping of their (103) planes (Fig. 3b). Such a significant elongation of the out-of-plane axis in the 5-BTO film in comparison to the 20-BTO film can also be ascribed to the well-known chemical expansion effect due to oxygen vacancies 15 .

Atomic-scale observation of polarization structures by STEM
Rich local fields including electrostatic forces and distorted lattice strains are usually accompany with oxygen vacancies 15 , we were thus intrigued to investigate the microscopic configuration of the electric dipoles in the vicinity of clustered oxygen vacancies in the ultrathin 5-BTO film by STEM. A large-scale low-angle annular dark-field (LAADF)-STEM image is shown in Fig. 4a, which displays cloudy contrast in BTO and suggests the existence of rich point defects, such as oxygen vacancies 39 . By further adjusting the contrast (Fig. 4b), the processed image shows cloudy regions of 10 nm or above in the BTO film. The contrast variation is consistent with the picture of oxygen vacancy clusters. Figure 4c displays a high-angle annular dark-field (HAADF)-STEM image of the 5-BTO/SRO/STO (001) heterostructure. The heterostructure shows sharp heterointerfaces with no extended defects throughout the films, demonstrating its high quality. The displacement vectors of the Ti cations with respect to the Ba cage in BTO are mapped in Fig. 4d with overlaid arrows. We note here that, due to rich oxygen vacancies, it is challenging to resolve oxygen positions via angular bright-field STEM, thus unable to acquire the local polarization vector accurately. Nevertheless, as the polarization of BTO is contributed by a dominant Ti displacement and a much smaller, opposite O displacement (relative to Ba), the geometry of Ti displacement vectors should resemble the polarization patterns to a large extent (see Supplementary Discussion) 42 . Contrary to the well-established picture of 180°domains with Ising-type domain walls in ultrathin films of uniaxial ferroelectrics, the electric dipoles in the 5-BTO film form complex rotation patterns (marked with the blue squares). Figure 4e shows the magnified vortex (Domain III) or half vortex (Domain I and II) 5,13 . It is noteworthy that these polarization patterns are prevailing throughout the sample in our STEM observation (Supplementary Fig. 6). To validate our STEM observation, we plot the distribution of the Ti displacement of both BTO and STO ( Supplementary Fig. 5). As expected for a non-polar material, STO exhibits nearly vanishing Ti displacement. This result certifies the quality of our STEM image. By contrast, the Ti displacement of BTO has a rather dispersed distribution in both magnitudes (ranging from nearly 0 up to 30 pm) and directions, in consistence with the vortex-like polarization patterns.

Origin of polarization nanodomains
To reveal the origin of such exotic polarization pattern, we analyze the local strain state by lattice constant mapping and the oxidation state of Ti by electron energy loss spectroscopy (EELS). Here we note that, although O-K spectrum can directly provide information of oxygen vacancies, the much lower sensitivity of oxygen to the electron probe makes it difficult for detecting. Nevertheless, oxygen vacancies can be traced via their footprints in oxides including local lattice strain and electron doping (usually accompany with a valence reduction of Ti from 4+ to 3+) 19,39 . Due to a coherently strained state of the heterostructure, the in-plane lattice strain (calculated by d ½100 À d ½100 À Á =d ½100 , where d ½100 is the lattice constant of each BTO unit cell measured from the STEM image and d ½100 is the corresponding lattice constant of bulk BTO) varies slightly throughout the BTO film (Supplementary Fig. 7). By contrast, the out-of-plane strain of the same region in Fig. 4d exhibits  significant inhomogeneity (Fig. 4f). Considering the direct connection between oxygen vacancies and lattice expansion, the inhomogeneous out-of-plane lattice constant distribution thus reflects the clustering state of oxygen vacancies. Further comparing Fig. 4d, f, we notice that oxygen vacancies tend to cluster at the vicinity of vortices. This phenomenon should be contributed by complex local fields accompanying oxygen vacancies. First of all, density functional theory calculations suggested that the local polarization vector neighboring a single oxygen vacancy can be reformed to point outward from the oxygen vacancy and away from the inherent polar axis due to lattice relaxation 22,23,43 . Second, the inhomogeneous strain coupled with oxygen vacancy clusters naturally gives rise to strain gradients and accompany flexoelectric fields, which also impact the arrangement of electric dipoles 32,33 . Lastly, a radial electrostatic field can be generated depending on the geometry of local charge arrangement, e.g., clustered ionized oxygen vacancies, which can contribute to the non-uniaxial polarization pattern as well 13 .
It is particularly worth mentioning that charged domain walls in the center of our BTO film are observed as well, which can also be induced by oxygen vacancies as reported previously 44 . The different patterns of electric dipoles induced by oxygen vacancies may depend on the specific geometries of oxygen vacancy clusters. For example, oxygen vacancy plates may give rise to charged domain walls 44 while spherical oxygen vacancy clusters can induce swirling electric dipoles (as Domain I in Fig. 4d). Another interesting finding is that the region where oxygen vacancies populate show diminishing polarization. This contrasts the general belief of strong polarization-tetragonality coupling, but in consistence with previous calculations that predict the polarization suppressing role of oxygen vacancies 24 . Therefore, this finding calls for particular attention to the relation between lattice strain and local polarization values under the existence of oxygen vacancies. Figure 4g plots the Ti-L 2,3 spectra of two representative BTO regions with different polarization states, namely region-ofinterest (ROI)-1, -2, and -3 in Fig. 4d. While all spectra show the characteristic t 2g -e g splitting peaks of L 2 and L 3 edges for Ti 4+ , an additional peak emerges between t 2g and e g in the spectrum of ROI-1 with large polarizations and has been ascribed to the Ti 3+ state 39 . This result indicates admixture of Ti 3+ and Ti 4+ in ROI-1. It is worth noting that this feature is prevailing throughout the film, namely regions with large polarizations, e.g., the periphery of polar vortices, exhibit a reduced Ti valence state; regions with diminishing polarizations, e.g., the core of polar vortices (ROI-3) and charged domain walls (ROI-2), exhibit pure Ti 4+ . The dissociation of the oxygen vacancy clusters (from the local strain analysis) and the Ti 3+ -accumulated regions suggests certain delocalization of electrons doped by oxygen vacancies, which could be driven by polarization screening 18,24 . This finding again verifies the existence of oxygen vacancies and further suggests considerable inhomogeneity due to oxygen vacancy clusters in 5-BTO films.

DISCUSSION
In the end, it is worth discussing the connection between the microscopic polarization patterns observed by STEM and those by PFM. The length scales of the polarization configurations by two methods may seem inconsistent, namely a typical scale of several nanometers in STEM and a scale of several to tens of nanometers in PFM. However, we note that PFM usually probes collective piezoelectric response over 100 nm deep into the target sample and evaluates the average local polarization under the tip 45 . Therefore, in the case of complex polarization structures, e.g., 180°domains with tail-to-tail domain walls and vortex-like polarization structures, PFM can only extract the sum of polarization from a large probing volume. Despite this drawback, we may still infer the existence of a complex polarization configuration in the 5-BTO film from PFM. For example, in Fig. 1d, both positively and negatively poled regions exhibited comparable PFM magnitudes whereas the pristine region with the same PFM phase as that of the positively poled region shows a much weaker magnitude. The weaker PFM magnitude can be interpreted as canceling of polarizations with opposite directions in an ensemble of complex polarization configurations such as rotating, head-tohead, or tail-to-tail polarization patterns, as observed by STEM 46 . Nevertheless, by STEM, we indeed observe some regions with a mainly upward polarization ( Supplementary  Fig. 6a). Moreover, by analyzing the average polarizations of each unit-cell column in Fig. 4d, we can also observe a region with an average downward polarization and a lateral dimension of~10 nm (lower panel in Fig. 4d and Supplementary Fig. 8). Therefore, our STEM and PFM results are mutually consistent. In this regard, STEM should provide a better knowledge of the polar structure, regardless of its limited observation area.
In summary, we report a close examination on the role of oxygen vacancies in the creation of swirling electric dipoles in a prototypical ferroelectric, namely BaTiO 3 , ultrathin film, which have been usually reported with a uniaxial polarization pattern due to the epitaxial strain. There are two important findings that have been overlooked in previous experimental studies. First, the dynamic clustering process of oxygen vacancies results in an evolution of a seemingly single domain to complex nanoscale polarization structures, as observed by PFM. This result calls for particular attention in understanding the self-polarization behavior of ferroelectric thin films with the presence of oxygen vacancies. Second, clustering oxygen vacancies can lead to nonuniaxial polarization structures such as vortices, which can be taken as embryo of topological polar objects. Notably, it was also reported that oxygen vacancies cluster at the core of vortex in the well-studied PbTiO 3 /SrTiO 3 superlattice 26 . However, the formation of vortices in the PbTiO 3 /SrTiO 3 superlattice has been regarded as a result of competition between the strain and the electrostatic boundary condition, and the positive role of oxygen vacancies has not been identified. Our study reveals a more active role of oxygen vacancies in producing rotating polarization structures. Moreover, the mobility of oxygen vacancies under external mechanical and electrical stimuli implies promising controllability over the coupled polarization structures 47 . From the aspect of practical applications, the electron-donor role of oxygen vacancies combined with the coupled topological polarization patterns should exert profound influences on the transport of ferroelectric oxide-based electronic devices, especially thin ferroelectric tunnel junctions and domain wall memresistors 19,40 . Therefore, our study should inspire future efforts on the exploration of topological polar objects by exploiting point defects.

Materials preparation and characterization
Details of the PLD growth have been reported elsewhere 29 . AFM and PFM were performed using a Cypher scanning probe microscope (MFP-3D, Asylum Research) with Ir-Pt-coated tips (PPP-EFM, Nanosensors). A Dual Alternating Current Resonance Tracing mode was used for PFM measurements. XRD was performed with a Bruker D8 Discover. XPS was carried out in a VG ESCALAB 220i-XL system using a monochromatic Al Kα source with the pass energy of 10 eV.

STEM characterization
Cross-sectional STEM specimens were prepared by focused ion-beam milling (Helios 650, FEI) with low-energy ion beams (<2 kV) at the final milling process and further thinned by focused low-energy Ar-ion milling (NanoMill 1040; Fischione Instruments). HAADF-STEM images for measuring atomic positions in this work were obtained by spherical aberration corrected STEM with a probe corrector (JEM-ARM 200F; JEOL Ltd) equipped with a cold field emission gun at an acceleration voltage of 200 kV installed at the National Center for Inter-university Research Facilities at Seoul National University (South Korea). The semiconvergence angle for STEM analyses was set to 24 mrad. The inner semi-collection angle range was set to 80 and 30 mrad for HAADF-and LAADF-STEM, respectively. For suppressing spatial incoherency and improving signal-to-noise ratio of HAADF-STEM images, a series of images was acquired, with a fast scan speed of 1 μs pixel −1 for 30 frames, and registered to the atomic-resolved STEM images. Both rigid and non-rigid registration methods were applied and showed consistent results. The atomic positions were determined by fitting and subtracting them in a regular sequence of atomic number as 2D-Gaussian function using customized MATLAB code. Then polarizations of each BTO unit-cell were determined by using the column offsets of each titanium from the center of Ba cage. For all images, the horizontal and vertical lines represent the [100] and [001] orientations, respectively. EELS analyses were carried out at the corresponding locations to investigate the relationship between polarization domain structure and electronic structure using an EEL spectrometer (Model 965 GIF Quantum ER, Gatan, USA) with an energy resolution of 0.8 eV.