Direct observation of polarization-induced two-dimensional electron/hole gases at ferroelectric-insulator interface

Two-dimensional electron gas or hole gas (2DEG or 2DHG) and their functionalities at artificial heterostructure interfaces have attracted extensive attention in recent years. Many theoretical calculations and recent experimental studies have shown the formation of alternating 2DEG and 2DHG at ferroelectric/insulator interfaces, such as BiFeO3/TbScO3, depending on the different polarization states. However, a direct observation based on the local charge distribution at the BiFeO3/TbScO3 interface has yet to be explored. Herein we demonstrate the direct observation of 2DHG and 2DEG at BiFeO3/TbScO3 interface using four-dimensional scanning transmission electron microscopy and Bader charge analysis. The results show that the measured charge state of each Fe/O columns at the interface undergoes a significant increase/reduction for the polarization state pointing away/toward the interface, indicating the existence of 2DHG/2DEG. This method opens up a path of directly observing charge at atomic scale and provides new insights into the design of future electronic nanodevices.


INTRODUCTION
Understanding the atomic-scale interface properties in heterostructures is critical for the engineering of oxide-based electronic devices. Such studies have recently become feasible due to advances in scanning transmission electron microscopy (STEM). The combination of high-quality thin film epitaxy and atomicresolution STEM has led to a deep understanding of a variety of novel phenomena emergent from interfaces, such as the twodimensional electron gas (2DEG) at the LaAlO 3 /SrTiO 3 (LAO/STO) interface. 1,2 Based on the polar catastrophe model, where alternatively charged sub-layers in LAO lead to free charge accumulation at the interface with neutral STO, the atomically sharp interface is crucial to the formation of the 2DEG. 3 Besides the LAO/STO interface, theoretical and experimental studies have shown that ferroelectric oxide interfaces can also host a 2DEG. [4][5][6][7][8] In cases where the ferroelectric has alternating charged sub-layers like the LAO/STO system, the same polar catastrophe mechanism is responsible for the appearance of the free-electron gas, but the spontaneous polarization further modifies the interfacial properties. When the polarization points toward the interface, additional electrons will move toward the interface to screen the positive polarization bound charge, resulting in an increased charge density of electrons. But when the polarization points away from the interface, electrons are pulled away from the interface, reducing the electron charge density. This polarization-induced difference can be big enough to generate a metal-insulator transition. First-principles calculations have predicted this phenomenon in KNbO 3 /ATiO 3 (A = Sr, Ba, Pb) 4 and BiFeO 3 /SrTiO 3 (BFO/STO). 5 If all sub-layers are charge-neutral, the spontaneous polarization alone will generate a free carrier gas at the interface where free holes or electrons aggregate to screen the negative or positive polarization bound charge, forming two-dimensional hole gas (2DHG) or 2DEG, respectively. In this case, the carrier gas at the interface can be switched between 2DHG and 2DEG by switching the polarization in the ferroelectric film. This has been predicted in BaTiO 3 /SrTiO 3 6 and PbTiO 3 /SrTiO 3 . 7 Experimental studies have shown indirect evidence indicating the existence of two-dimensional (2D) free carrier gases at BFO/TbScO 3 (TSO) 8 interfaces using conductive atomic force microscopy (C-AFM) and electron energy loss spectroscopy (EELS), supported by theoretical calculations. Energy-loss near-edge structure (ELNES) determined by the bonding characteristics, including bond length, symmetry, and coordination number, have been correlated with changes in formal oxidation state, but it does not measure the charge distribution which is a continuum. 9 Therefore, there has not been a direct measurement of the local charge distribution in the 2DHG/2DEG.
Recently, it has been shown that differential phase-contrast imaging and four-dimensional scanning transmission electron microscopy (4D STEM) can be used to directly image the electric field and charge distribution with sub-Å resolution without significant computational expense, 10-16 making it well suited for studying 2DHG and 2DEG. Unlike ELNES, the charge distribution measured from 4D STEM is derived from the probe's interaction with local electric field in the sample 10 which avoids correlations with spectroscopic features, thus simplifying interpretation of the results. By processing 4D STEM data with Bader charge analysis 12,17,18 , the boundaries of atomic columns in the continuous charge density map can be delineated and the charge states of individual atomic columns in ferroelectric heterostructures can be calculated. This technique offers the capability of directly showing how the electron distribution is coupled to the local polarization states at ferroelectric interfaces and how this leads to the formation of a 2DHG or 2DEG at ferroelectric/insulator interfaces, which are a critical part of understanding ferroelectric materials' properties, and provide guidance for developing nanoelectronic devices.
Here, we studied the charge distribution at a BFO/TSO interface using atomic-resolution STEM, 4D STEM, Bader charge analysis and EELS, and demonstrated the build-up of 2DHG/2DEG that is induced by upward (points away from the interface)/downward (points toward the interface) polarization. The specimen used is a 20

RESULTS AND DISCUSSION BiFeO 3 domain structures and chemical compositions
Atomic models of BFO unit cell (Fig. 1a) show the oxygen octahedra and the central Fe are displaced from their respective positions, resulting in a large spontaneous polarization along the <111 > p directions. 19 The rotation angle between different polarization directions can be 71 , 109 , or 180 , yielding three types of domain walls (DWs). [20][21][22][23] In 109 domain arrays, the alternating polarization states of these domains have their out-ofplane polarization component pointed either toward or away from the interface, inducing positive or negative bound charges, respectively, as shown in the schematic in Fig. 1b. Such domain structure is confirmed by dark-field TEM ( Supplementary Fig. 1) and STEM (Fig. 1b), in which the contrast clearly shows the crosssectional view of the high-quality BFO/TSO epitaxial thin films with ordered 109 domain arrays. An atomically resolved high-angle annular dark-field (HAADF) image, shown in Fig. 1c, and the corresponding atomic-resolution X-ray energy-dispersive spectroscopy (EDS) maps of the BFO/TSO interface reveal that the interface is atomically sharp and BFO's termination layer is Fe-O 2 . Oxygen columns cannot be seen in HAADF STEM image is because of their weak scattering of electrons.

Polarization and charge measurements
To determine the BFO polarization directions at the interface, high-resolution HAADF images of two local regions across the BFO/TSO interface (marked by white rectangles d and e in Fig. 1b) were collected as shown in Fig. 2a, d. We define a vector D FB , which is the displacement of the Fe atom from the center of the unit cell formed by the four adjacent Bi atoms. Because the D FB is always pointing toward the center of the oxygen octahedra, it is opposite to the polarization vector P that is in the same plane. Thus −D FB vector can be used to estimate P. 24 From Fig. 2a, d, Bi and Fe atomic positions were measured by using 2D Gaussian fittings on their atomic columns to locate the centers and record the offsets. 24 Fig. 2b, e shows the −D FB vector mapping that correspond to Fig. 2a, b, indicating an upward polarization state in Fig. 2b and a downward polarization state in Fig. 2e, respectively.
To study the charge distribution of the BFO/TSO interface with upward/downward polarization states, 4D STEM was performed on the same regions ( Fig. 2a, d) and the electric field and charge density were calculated from the 4D dataset (see "Methods" for details). Since the calculation of charge density from 4D STEM measurement is only valid when the sample is on zone axis and has thickness less than 6.6 nm 13 , the sample was carefully tilted to zone axis and thickness of the two regions was measured based on the position-averaged convergent beam electron diffraction (PACBED) patterns that are shown in Supplementary Fig. 5a. The thickness for both regions is measured as 6 ± 2 nm (see "Methods" for details), within the acceptable range. Real-space 2D chargedensity mappings derived from the 4D STEM data are shown in Fig. 2c, f. The positive nuclei, screened by the core electrons appear as positive charges in the charge density map, while the valence electrons appear as negative charge between the atomic columns. These features are also convolved with the probe size (0 .6 Å full-width at half-maximum). Due to the contribution of the nuclei, the maps show not only Bi and Fe columns but also oxygen columns, which is quite important for investigating oxygen octahedral rotation in ferroelectric materials. 25,26 However, due to a very large oxygen octahedral rotation in BFO (13.8 ) 19 , the oxygen columns are vague and elongated, even appearing merged with Bi columns. The atomic models are overlaid on the maps to confirm the positions of Bi, Fe, and O atoms.

Bader charge analysis on Fe/O columns
To investigate the possible existence of 2DHG or 2DEG at the BFO/ TSO interface, we measured the total charge of each Fe/O column by performing Bader charge analysis on the BFO side based on the charge density maps in Fig. 2c, f. In conventional Bader charge analysis used in first principles calculations, the saddle point in the 3D charge density map is defined as the boundary around each atomic column's site to separate the electrons that belonging to different atomic columns; within the boundary, the charge density is integrated to give a partial charge/valence charge for each atomic column. 17,18 Recent work has adapted technique for the projected 2D charge density maps that can be acquired from 4D STEM data. 12 As explained in previous publication and in the "Methods" section, the local minima of the charge density surrounding each atomic column in the 2D charge map corresponds to the saddle point in a 3D charge map and the charge state of the atomic column can be determined by integrating within this boundary.
As shown in Fig. 3a, b, each Fe/O column is defined by a black circle, including both the nuclei and the surrounding electrons.
The charge of each column that is 3 unit cells or further away from the interface is integrated and plotted in the histograms in Supplementary Fig. 3. The fitted curves show that the charge of the Fe/O column centers on 0.31e and 0.30e for upward and downward polarization states, respectively. By averaging the total charge of all the Fe/O columns in each unit cell row, the measured charge on Fe/O columns for the upward and downward polarization states as a function of a distance from the BFO/ TSO interface are plotted in Fig. 3c. The most remarkable feature is a change of the measured charge on Fe/O columns in the interfacial region of the BFO film that extends 2 unit cells from the interface. For the interfaces with upward polarization state, the measured charges on 0th and 1st unit cells are higher than that on the subsequent 8 layers which are far away from the interface. The charges in the first three unit cells decreases moving away from the interface; it is about 0.42e at 0th layer and then gradually decreases to 0.3e at 3rd layer. The measured charge stays around 0.3e for the remaining layers, which is in line with what has been previously reported from 4D STEM measurements of a BFO/STO interface. 12 These findings reveal the existence of a 2DHG in the first two unit cells on the BFO side. For the interfaces with downward polarization state, the measured charges on first two unit cells (0th: 0.17e, 1st: 0.23e) are lower than the rest (~0.3e) and the charge in the first three unit cells increases moving away from the interface, indicating there is a 2DEG. The same measurement was also applied to the Sc/O columns located on the TSO side of the interface, as shown in Supplementary Fig. 4. The result reveals the charge measured on Sc/O columns has no significant difference between the two polarization states, meaning the 2DHG/2DEG are located at the interface but primarily on the BFO side. In principle, bound charges at the interfaces may also contribute to the charge  Fig. 1b (a) and corresponding maps of −D FB vectors overlaid on the same HAADF STEM image (b) and charge density with atomic model overlaid (c). d-f An atomically resolved HAADF STEM image of downward polarized domain at the interface from the region marked by the white rectangle d in Fig. 1b (d) and corresponding maps of −D FB vectors overlaid on the same HAADF image (e) and charge density with atomic model overlaid (f). The HAADF STEM images and charge density maps are acquired separately. Scale bar, 1 nm. measurements on Fe/O or Sc/O columns. However, the negative/ positive bound charge induced by upward/downward polarization state will reduce/increase the amount of charge measured at the interfaces. This phenomenon is opposite to the experimental results that are shown in Fig. 3c, in which the interfaces with upward/downward polarization state have more/less charges, meaning that measurements of the charge of on Fe/O or Sc/O column are dominated by the free charges at the interface rather than bound charges.
Electron energy loss spectroscopy measurement of iron valence state EELS line scans across similar BFO regions to that shown in Fig. 2a, d were performed to confirm the existence of the 2D gases, as plotted in Fig. 4. The EELS fine structures in Fig. 4a show the Fe L 2,3 edge energy shift on the 0th and 1st unit cells as indicated by green arrows. The energy onset difference between O K and Fe L 3 edges is used to estimate the Fe oxidation state because of their linear relationship. 27 Fig. 4b shows that the energy onset difference is higher in the first two unit cells than the subsequent 8 layers for the interfaces with upward polarization state, indicating a higher Fe valence state at the interface than in the bulk. Such high Fe valence state is attributed to an accumulation of the free holes and therefore confirms the existence of the 2DHG. In contrast, for the interface with downward polarization state, the energy onset difference at the interface is lower than that in the bulk, indicating a lower Fe oxidation state and thus the existence of 2DEG. One may notice that the 2D gases are confined in the first two unit cells in our experiments but extend to 4 unit cells in previous studies by Zhang et al. 8 . This discrepancy is caused by the differences in BFO film thickness. Our BFO film is about 20 nm thick along the out-of-plane direction and its depolarization field is larger than that of the 400 nm thick BFO film used by Zhang et al. A larger depolarization field will result in smaller BFO polarizations and thus less build-up of bound charges at the interfaces. As a result, less free charges will be attracted, and the 2D gases will be contained in smaller regions in a 20 nm BFO film than in a 400 nm BFO film.
To summarize, in ferroelectric BFO thin films grown on TSO substrate, we have directly observed the build-up of 2D free carrier gases that are induced by spontaneous polarizations at BFO/TSO interface by using 4D STEM and Bader charge analysis to quantify the local charge distribution. Based on the analysis of the atomic-resolution charge density mapping and EELS analysis, we found that the upward polarization state will induce 2DHG while downward polarization state will induce 2DEG, both within the first two unit cells adjacent to the interface. This result further demonstrates the sensitivity of 4D STEM to both positive and negative free charge carriers in thin film heterostructures. This technique opens up a path of directly observing charges in atomic scale, contributing to establishing the design paradigm for nextgeneration electronic nanodevices.

METHODS Thin film growth
The 20 nm thick BiFeO3 films were grown on single crystal (110) O TbScO 3 surfaces by reactive-molecular beam epitaxy (MBE) in an EPI 930 MBE chamber equipped with reflection high-energy electron diffraction (RHEED). The Fe flux was about 2.0 × 10 13 atoms/(cm 2 ·s) and the Bi flux was 1.1 × 10 14 atoms/(cm 2 ·s). During the deposition, distilled ozone (~80% ozone) were used to create a background oxidant pressure of 1 × 10 -6 Torr. The substrate was maintained at a constant temperature of 610°C.

TEM sample preparation and characterization
TEM specimen was prepared by mechanical polishing followed by argon ion milling. HAADF STEM, 4D STEM, and EDS experiments were carried out on a JEOL Grand ARM300CF equipped with a cold field emission gun and double spherical aberration correctors with a spatial resolution of~0.6 Å operating at 300 keV in Irvine Materials Research Institute at the University of California, Irvine. The microscope was operated with a convergence angle of 32 mrad, a collection angle at 90-165 mrad, and a beam current of 36 pA. EDS maps were acquired using dual silicon-drift detectors (SDDs). 80 scans (each with a 0.1 ms dwell time and 0.05 nm pixel size) in the same area across the interface were summed. The 4D STEM data were collected with a Gatan OneView camera at a speed of 300 frames per second and a scanning step size of 0.4 Å.

Image processing
In order to improve the signal-to-noise ratio and more clearly show the polarization states, the STEM images in Fig. 2a, d were filtered in Fourier space by including all the diffraction frequencies and smoothed using a commercial software DigitalMicrograph. In the charge density maps, the drift along the horizontal scan direction was corrected by measuring the average displacement per line from a conventional fast-scanned STEM image and then shifting each line back. The vertical drift was corrected by rescaling the vertical axis to fit with the fast-scanned STEM image.

Thickness measurement
Sample thickness and orientation have significant influence on the intensity in CBED patterns. 28 Therefore, we used PACBED patterns to align the sample to the zone axis and estimated the sample thickness. PACBED patterns were collected by acquiring the average CBED patterns in diffraction plane while the beam scanned over the region on BFO side (Fig.  2a or d). During this process, a smaller convergence angle of 10.6 mrad was used because it made the sample tilt more apparent than in standard 33 mrad imaging condition, enabling precise alignment of the sample's zone axis 12,28 . Then the experimental PACBED patterns were compared with simulated PACBED patterns generated from multislice simulations using model structures with thickness ranging from 2 to 10 nm 29 . The intensity in experimental and simulated patterns was normalized to the same range and then compared quantitatively using least squares error to determine the sample thickness ( Supplementary Fig. 5). Both curves show lowest error at 6 nm, meaning the thicknesses for the regions in Fig. 2a, d are both 6 ± 2 nm.

Electric field and charge calculation
When transmits through a sample, the electron probe in the diffraction plane will be shifted negatively proportional to the local electric field, owning to the momentum change that is induced by the electric field. To quantify this, the center of mass (COM) of each diffraction pattern in the 4D STEM data was calculated. For each pattern, the deviation of the COM from the center of the diffraction pattern yields a vector field. Then we convert the unit of the pixel size in momentum space to mrad/pixel. The momentum change of electrons is based on where p is the momentum change, E is the electric field, and T is the time needed for one beam electron to pass through the specimen and can be calculated as where Δt is the sample thickness, v z is the speed of the electron along the beam direction. Based on the assumptions outlined in refs. 11,12 , we can rewrite the momentum change of electrons as where Δp xy is the momentum change in x-y plane, e is the charge of an electron and E xy is the electric field in the x-y plane. Because the diffraction patterns have the probe's momentum space information, Δp xy can be calculated based on the shift of the COM measured in the diffraction pattern where ΔCOM is the shift in the COM from the center of the diffraction pattern and p z is the momentum of the electron beam along the beam direction. Thus, the electric field can be calculated as Then we can calculate the divergence of the electric field (electric field maps are shown in Supplementary Fig. 2) and determine the charge density by using Gauss's law where ρ is the charge density and ε 0 is the vacuum permittivity. To fit the 2D map, charge density can be calculated by integrating both sides along the z axis.

Bader charge analysis
Based on the charge density map, we can calculate the total charge of each Fe/O and Sc/O atomic column by integrating the charge density within a region around each column. To define the size of the integration region, we applied Bader charge analysis, which is a method from DFT for calculating the charge state of atoms. 17,18 In conventional Bader charge analysis, the saddle point in the 3D charge density map is defined as the boundary around each atomic column's site to separate the electrons that belonging to different atomic columns; within the boundary, the charge density is integrated to give a partial charge/valence charge for each atomic column. In our case, 4D STEM measurements are generally projections of the sample along the beam direction, so it is not possible to calculate a saddle point as it would be in the full 3D charge distribution. However, as outlined by ref. 12 , the process can be adapted to work for 2D images: the minimum of the charge density surrounding each atomic column corresponds to the saddle point in the full 3D charge distribution. We first defined the charge density local minimum on left, right, top, and bottom sides for each atomic column. The minima on the left and right sides were determined based on a horizontal line profile of the charge density through the center of the atomic column. The minima on the top and bottom side were determined in the same manner with a vertical line profile. The width of the line profile is the averaged diameter of all the Fe/ O or Sc/O atomic columns measured from the corresponding HAADF STEM images. Then we defined a circular integration region as shown in Fig. 3a, b. Its radius is determined by averaging the distance to the four local minima. The slightly different sizes of the integration regions are to be expected in experimental data. Slight variations from one atomic column to another, scanning noise, and detector noise will all effect the exact shape of the region and thus the total integrated charge assigned to each atom. We quantified these uncertainties by calculating the standard deviation for each Fe-O and Sc-O layer and use them as error bars in Fig. 3c and Supplementary Fig. 4.

Electron energy loss spectroscopy
High spatial resolution EELS was carried out on an aberration-corrected monochromated Nion ultra-STEM that operated at 60 kV with a convergence semi-angle of 30 mrad and a beam current of 100 pA. The spectrum image was collected with a dispersion of 0.164 eV/channel and a dwell time of 1 s/pixel. The spectrum background was removed by fitting a power-law function to the pre-edge using the commercial software package DigitalMicrograph. The energy loss where the edge reaches 10% of its maximum intensity is taken as the energy onset value for each O K and Fe L 3 edge.