Revealing the three-dimensional arrangement of polar topology in nanoparticles

In the early 2000s, low dimensional ferroelectric systems were predicted to have topologically nontrivial polar structures, such as vortices or skyrmions, depending on mechanical or electrical boundary conditions. A few variants of these structures have been experimentally observed in thin film model systems, where they are engineered by balancing electrostatic charge and elastic distortion energies. However, the measurement and classification of topological textures for general ferroelectric nanostructures have remained elusive, as it requires mapping the local polarization at the atomic scale in three dimensions. Here we unveil topological polar structures in ferroelectric BaTiO3 nanoparticles via atomic electron tomography, which enables us to reconstruct the full three-dimensional arrangement of cation atoms at an individual atom level. Our three-dimensional polarization maps reveal clear topological orderings, along with evidence of size-dependent topological transitions from a single vortex structure to multiple vortices, consistent with theoretical predictions. The discovery of the predicted topological polar ordering in nanoscale ferroelectrics, independent of epitaxial strain, widens the research perspective and offers potential for practical applications utilizing contact-free switchable toroidal moments.


Introduction
Ferroelectric crystals have the ability to exhibit spontaneous polarizations that can be switched by external electric fields.A range of physical properties is related to ferroelectricity, including piezoelectricity, pyroelectricity, flexoelectricity, and dielectricity, which make ferroelectrics useful for various electronic devices, such as memories, sensors, and capacitors [1][2][3][4][5][6][7] .Typically, bulk ferroelectric crystals exist in a paraelectric state at high temperatures, but as the temperature decreases below the Curie temperature, they make a transition into ferroelectric states.The ferroic states are often accompanied by multiple domains with varying spontaneous polarizations, and the properties of ferroelectrics are significantly affected by the internal domain structures.
Ferroelectric nanostructures differ from their bulk counterparts; uncompensated charges at their surfaces create a depolarizing field, which can be strong enough to suppress spontaneous polarization and prevent ferroelectric nanostructures from exhibiting a ferroelectric state even at low temperatures 8- 11 .To induce polarizations, external electric fields or conducting boundary conditions are required to screen the depolarizing field 12,13 .Therefore, the occurrence and behavior of ferroelectric ordering in low-dimensional systems have been a topic of interest.
With the development of effective Hamiltonian and phase-field methods, the existence of vortices 14,15 or other topological structures such as flux-closures 16 , bubbles 17 , skyrmions 18 , merons 19 , polar waves 20 , and topological defects 21 was predicted from simulations for different types of nanostructures including films and nanodots.Many experimental efforts were also devoted to fabricating various ferroelectric nanostructures and observing polar topologies, and the use of advanced transmission electron microscopy (TEM) and piezoresponse force microscopy techniques has enabled the observation of a periodic array of full flux-closure quadrants in PbTiO3/SrTiO3 multilayer/superlattice systems 22 , as well as other predicted structures 23,24 .
However, these techniques are often restricted to surfaces or lower dimensional projections, or resolutions at a non-atomic scale 25 , and are inadequate for probing three-dimensional (3D) polarization structures within the ferroelectric nanostructures.Due to this limitation, most of the polar topology structures have been determined from cross-sections of thin films for which the epitaxial strain plays a crucial role.Although the existence of polar vortices in ferroelectric nanostructures, including nanodots, nanowires, and nanodisks, has been predicted without the need for epitaxial strain 14,15,26 , experimental observation of such ordering and their size-dependent structural transitions has not yet been achieved.Furthermore, a complete classification of 3D topological textures would require a true 3D mapping of local polarizations, because predicting higher dimensional structures from lower dimensional projections can often be misleading 27,28 .Therefore, it is necessary to directly image the 3D polarization structures to fully understand the polar topology.
Typically, ferroelectric polarizations are directly related to local atomic displacements within the unit cell (inversion symmetry breaking).If full 3D atomic coordinates of all the atoms within ferroelectric nanostructures can be determined at high precision, atomic scale 3D local polarization can be fully mapped.Here we report the experimentally observed topological ferroelectric ordering within BaTiO3 nanoparticles of approximately 8.8 nm and 10.1 nm diameters at an unprecedented 3D individual atom level, which was achieved via precise determination of the 3D coordinates of cation atoms using atomic electron tomography (AET) [27][28][29][30] .

Determination of 3D atomic structures of BaTiO3 nanoparticles
The BaTiO3 nanoparticles were fabricated using hydrothermal synthesis.Their structure and composition were confirmed using energy dispersive X-ray spectroscopy (EDS) mapping and powder X-ray diffraction (PXRD) (see Supplementary Fig. 1).After drop-casting the particles onto a SiN membrane, an amorphous carbon layer of approximately 2.0 nm thickness was deposited, followed by a vacuum cleaning of 10 hours at 150 °C (Methods).Two nanoparticles, with average diameters of ~8.8 nm (Particle 1) and ~10.1 nm (Particle 2), were selected for annular dark-field scanning transmission electron microscopy (ADF-STEM) tilt-series measurements to conduct AET (see Supplementary Fig. 2 and Methods).The GENFIRE algorithm 31 was used to reconstruct the atomic resolution 3D tomograms of both nanoparticles (Methods).Note that the total electron dose for ADF-STEM data acquisition was optimized to minimize the beam damage (see Supplementary Fig. 3) and the contrast from oxygen atoms was too weak to be detected under the given dose condition.
The 3D atomic coordinates and chemical species of Ba and Ti atoms in each nanoparticle were determined using atom tracing and classification techniques, with a precision of 38.3 pm for Particle 1 and 34.4 pm for Particle 2, as determined via multislice simulations (Methods).The reliability of the structures was supported by the reasonable level of R-factors, which were calculated by comparing the measured tilt-series data with the forward-projections from the determined atomic structures (Methods).
The resulting 3D atomic models of Particle 1 and Particle 2, depicted in Fig. 1a and Supplementary Fig. 4a, revealed the alternating Ba and Ti atomic layers along the <001> directions.Figure 1b-d and Supplementary Fig. 4b-d further show the 1.07 Å thick internal slices of the 3D tomogram along the [001] (Fig. 1b-d) or [010] (Supplementary Fig. 4b-d) directions for three consecutive atomic layers near the core of the particles.The intensity distribution of the tomogram slices and traced atom positions clearly represent the atomic planes of square lattices, consistent with expected ABO 3 perovskite cation ordering.Over 99.5% of the 3D cation atomic coordinates obtained could be successfully fitted into a body-centered cubic (bcc) lattice (Methods), which is consistent with the expected cubic phase of BaTiO3 with sizes smaller than 30 nm 32 .
Before the 3D displacement analysis, we first conducted the two-dimensional (2D) Ti displacement mapping of linearly projected 3D volumes (along [100], [010], and [001]) generated from the experimentally determined 3D atomic structures (Fig. 1e-g and Supplementary Fig. 4e-g).The grayscale background indicates the projected intensity, where Ba and Ti columns can be distinguished from the relative difference in contrasts (Ti columns are weaker).To quantify the local Ti displacement directions, we calculated the deviation of the identified Ti column positions from the geometric centers of the four neighboring Ba column positions 33 .Note that the projections along the [100] direction exhibit minor indications of vortex-shaped ordering (Fig. 1e and Supplementary Fig. 4e), but no such ordering can be found from the projections along the other directions (Fig. 1f-g and Supplementary Fig. 4f-g).We also conducted several 2D-ADF-STEM measurements for other BaTiO3 particles along <001> zone axes, and they also showed no vortex-like features (Supplementary Fig. 5), consistent with previous 2D-projection-based studies 34,35 .

3D local polarization maps of the nanoparticles
Figure 2 and Supplementary Fig. 6 present sliced maps showing the 3D distribution of Ti atomic displacements for the Particle 1 (8.8 nm BaTiO3).In Fig. 2a, we show the in-plane atomic displacement direction maps of representative Ti atomic layers sliced along the [100] direction.The displacements were calculated as the deviation of each 3D Ti position from the geometric center of the eight neighboring Ba atomic positions and subsequently averaged using a Gaussian kernel (Methods).Note that the application of kernel-based local averaging effectively improves precision, reducing the error in the locally-averaged displacement field to about 4.9 pm for Particle 1 and 4.3 pm for Particle 2. This level of precision is adequate to describe the global polar ordering behavior, considering the typical ferroelectric local displacement level (tens of picometers).However, as a trade-off, it results in a slightly lower 3D spatial resolution of approximately 0.9 nm 30 .The polarization maps before kernel averaging can also be found in Supplementary Fig. 6.
The in-plane displacement direction maps clearly reveal that unidirectional polar order dominates in the top half of the particle while a counterclockwise vortex exists in the bottom half.This corresponds to the theoretically predicted polar-toroidal multiorder (PTMO) state, in which the polar and toroidal orders coexist 36 .In real samples, due to the presence of various screening sources such as charge carriers arising from oxygen vacancies, the electrostatic boundary condition can exhibit a partially screened state, rather than being perfectly short-circuited or fully open-circuited.In this case, the vortex core can be off-centered rather than being at the center of the particle, resulting in toroidal order near the vortex center and polar ordering away from it.Note that no clear vortex structures can be found if the particle is sliced along [010] or [001] directions (Supplementary Fig. 6), indicating that the vortex line is predominantly oriented along the [100] direction.
In Fig. 2b, the magnitudes of the kernel-averaged 3D Ti atomic displacements and calculated local polarization values are plotted (Methods).The range of 3D Ti displacements (and the resulting polarization values) is consistent with the previously observed values in the cubic phase of nanosized BaTiO3 37 .The region near the core of the vortex (indicated by blue arrows) shows relatively smaller local polarizations compared to other regions where the vortex does not appear, which is also in agreement with the result from phase-field simulations 26 .
Based on the 3D displacement vectors of all Ti atoms, we calculated the average displacement vector, which was found to be 18.7 ± 0.50 pm in magnitude along the direction of (0.27, -0.96, 0.07), resulting in the estimated net polarization per unit cell of 0.40 ± 0.13 C m -2 (Methods).This value is slightly larger than that of bulk BaTiO3 38 (0.25 C m -2 ).Several experimental and simulation studies have indicated that the magnitude of Ba-Ti displacements increases as the size of the system decreases up to a certain point (beyond which the ferroelectric ordering can be suppressed) 35,37,39 .For instance, a mean Ba-Ti displacement of ~18.9 pm was observed from a 26 nm BaTiO3 nanoparticle 35 , and ~30 pm was found from a 4 × 4 nm 2 BaTiO3 nanocluster 37 .Considering that these experimental observations are for 2D projected displacements, our measured 3D values are consistent with these previous studies.Having the full 3D local polarization distribution further enables us to calculate the net toroidal moment vector per unit cell G (Methods).The magnitude and direction of the vector were 0.22 ± 0.07 e Å -1 and (0.89, -0.40, 0.20), respectively.Here, the direction of the toroidal moment is not fully aligned with the [100] direction.This is also responsible for the PTMO state mentioned above, resulting in the shift of the vortex center positions in different 2D slices (eventually hitting the surface of the particles), as shown in Fig. 2. Similar analyses were conducted for the Particle 2 (10.1 nm), as described in Fig. 3 and Supplementary Fig. 7. Notably, the in-plane displacement direction maps along the [100] direction of this particle show multiple vortices, as can be seen in multiple consecutive layers near the center of the nanoparticle (Fig. 3a).This behavior is precisely what is expected for ferroelectric nanoparticles, transitioning from a single vortex structure to multiple vortex structures as the particle size increases 26 .The magnitudes of 3D Ti atomic displacements and local polarization values also exhibit reasonable values, with expected suppression near the core of the vortices (Fig. 3b).
The average displacement vector for the Particle 2 was found to be 29.0 ± 0.34 pm in magnitude along the direction of (0.61, -0.78, 0.10), and the net polarization per unit cell was estimated to be 0.61 ± 0.20 C m -2 .These values are slightly larger than those of Particle 1, but still in agreement with the values in the literature 37 .For Particle 2, the magnitude and direction of the net toroidal moment vector per unit cell G were 0.51 ± 0.16 e Å -1 and (0.88, -0.46, -0.04), respectively.While the direction is more or less consistent with that of Particle 1, the magnitude is about twice as large as that of Particle 1.However, the distribution of the magnitude of local toroidal moments is similar for both particles (Supplementary Fig. 8a-b).Therefore, the difference in the net moments originates from the fact that the local toroidal moments of Particle 2 show a stronger tendency to align parallel to the direction of the net moment (Supplementary Fig. 8c-d).Moreover, the contribution of the surface atoms to the net moments is higher compared to that of the core atoms for both particles (Supplementary Fig. 8e-f), resulting in a larger net moment for the larger particle, which has a greater number of surface atoms.
Since the nanoparticles show substantial cation displacement and consequent local polarization which typically accompany local tetragonal distortion, we further performed local structural analysis regarding the local tetragonality (Methods).Although no long-range tetragonal distortion can be observed from the PXRD analysis (see Supplementary Fig. 1g and Methods), there are localized areas with high tetragonality that cannot arise from a cubic phase (Supplementary Fig. 9).Our findings, where the reduced size of BaTiO3 exhibits distortion in the short-range but simultaneously shows a cubic phase in the long-range, are consistent with the results in previous studies 40 .Additionally, a comparison of the tetragonality map (Supplementary Fig. 9) with the 3D displacement magnitude map (Figs. 2 and 3) reveals that areas with low tetragonality (i.e., c/a ratio being close to 1) can also exhibit large Ti atomic displacement and resulting polarization.In nanostructured ferroelectric systems, unlike in bulk materials, an unusually strong atomic displacement has been observed in previous STEM studies, and it does not always correlate with the tetragonality, consistent with our findings 41,42 .

Effective Hamiltonian simulations
To corroborate the experimental findings, we conducted effective Hamiltonian simulations for similarly sized nanoparticles (Particle 1 � and Particle 2 � ; see Methods).For the Particle 1 � whose size is comparable to that of Particle 1, the toroidal moment is predicted to be 0.16 e Å -1 in magnitude along (-0.81, 0.44, -0.39), which agrees reasonably well with the experiment in terms of both magnitude and direction (sign of each component is not important here because there are eight degenerate solutions of equal free energy arising from the + or -sign).For Particle 2 � (similar in size to Particle 2), the calculations give the toroidal moment of magnitude 0.18 e Å -1 (about a third of the experimental value) along the direction of (0.86, 0.44, -0.25).
Note that the toroidal moment, in both simulations, depends on temperature and grows in amplitude, when cooling the systems down (Supplementary Fig. 10), eventually reaching (-0.20, 0.20, -0.19) e Å - 1 and (0.21, 0.22, -0.22) e Å -1 for Particle 1 � and Particle 2 � , respectively, at 5 K.This confirms that the larger particle exhibits a greater magnitude of toroidal moment than the smaller one, as observed in the experiments, although the difference is somewhat less pronounced.The discrepancies between the experiment and simulation can be attributed to the fact that the calculations were performed for ideal nanoparticles without considering surface reconstruction, defects, and surface adatoms (carbon atoms can easily contaminate the surface during the experiment, which cannot be detected in our ADF-STEM images).Moreover, the depolarizing field results in the absence of polarization in our simulations, contrasting with this and other experiments 34,35,43 which show the presence of polarization.This is not surprising because there can be an effective screening of the depolarizing field in the experimental system, which may result from charge carriers induced by the presence of oxygen vacancies [44][45][46] , partially screened electrostatic boundary condition, or/and lateral strain, leading to the aforementioned PTMO states.Such a comparison tells us that the performed experiment brings a new challenge for the future theoretical investigations.
We also found a single vortex structure in the Particle 1 � and multiple vortices in the Particle 2 � from the 2D slices perpendicular to the [100] direction, as illustrated in Fig. 4, showing an excellent agreement with the experimental findings (Figs. 2 and 3).These vortices arise from a connection of the domains having different directions of the local polarization.For example, the vortex in the Particle 1 � connects 4 domains having in-plane dipoles along  � +  �,  � −  �, − � −  �, and − � +  �, where  � and  � are unit vectors along the y-and z-axis, respectively.Moreover, the presence of more vortices in the Particle 2 � tells us that the domain structures observed in these experiments and simulations are more complex than the texture resulting from a simple elongation of dots 47,48 .

Discussion and outlook
In this work, we have determined the full 3D atomic coordinates of cation atoms for ferroelectric BaTiO3 nanoparticles of two different sizes via AET.This has unveiled the theoretically predicted toroidal ordering of local polarizations in the absence of epitaxial strain, together with a hint of the sizedependent phase transition in terms of the number of vortices.Our approach provides a full 3D polarization distribution for ferroelectric nanostructures, allowing for many more analyses to be conducted on various nanosized ferroelectric systems.For instance, topological polar structures and their transitions can be experimentally measured and classified by tuning the size, shape, and surface boundary conditions of nanostructures, leading to a better understanding and control over topological orderings.Additionally, a curl of in-plane polarization (the curl of the normalized in-plane Ti displacement fields for each atomic plane; see Methods) can be calculated to qualitatively visualize the tendency of the local rotations of polarization vectors in the 2D slices through the particles.A few representative maps are shown in Supplementary Fig. 11, where local curl vectors exhibit large magnitudes near the vortex cores of both particles.Interestingly, the curl vectors in Particle 2 alternately flip the direction from +x to -x multiple times around a small region, and the polarization ordering pattern near this region certainly resembles that of an antivortex (Supplementary Fig. 11).In addition to the vortices and antivortices, convergent hedgehog-type structures can also be found in some regions of Particle 2 (Supplementary Fig. 11).The difference in the number of vortices (one vortex for Particle 1 and two vortices for Particle 2) originates from the size-dependent topological transitions from single vortex-like structures to those of multiple vortices, as predicted from effective Hamiltonian and phasefield simulations 15,26,49 .In particular, in Particle 2, an antivortex is formed between two neighboring vortices of identical orientation (both are counterclockwise in this case) to create an energetically stable configuration, as observed in other ferroelectric 50,51 , ferromagnetic 52 , and superconducting 53 systems.
Our 3D polarization measurements also provide the out-of-plane components of the local polarizations at the vortex cores, which show clear out-of-plane components along the [100] direction at the vortex cores of both particles (Supplementary Figs.12-14), indicating that the vortices are chiral.The normalized helicity density (Hn) maps also clearly show non-zero helicity density values near the vortex centers, confirming their chirality (Supplementary Fig. 15).The helicity density map further allows us to determine the handedness of each chiral vortex by examining its sign 54 .Positive helicity density can be clearly observed for all the vortices appearing in both Particles, indicating that all the vortices are right-handed, which is also evident from the fact that the direction of the curl is parallel to the axial component of the polarization (Supplementary Figs.12-15).The antivortex and hedgehog-type structures are expected to show zero helicities in ideal configurations 54 , consistent with our experimental findings (Supplementary Fig. 15).Although the existence of vortices and antivortices, along with their chirality/handedness, can be qualitatively analyzed based on the curl of in-plane polarization and helicity density maps, these methods cannot quantitatively provide polar topological information of measured 3D polarization distributions.Therefore, we further calculated topological invariants (e.g., winding number, skyrmion number, and Hopf invariant) with the circle-equivalent and sphere-equivalent order-parameter spaces 54 (Methods).First, we calculated the topological invariants defined for 2D polarization vector fields (2D vectors in a 2D plane).From the 2D polarization vector fields sliced from our 3D polarization distributions along [100] direction, local winding numbers were calculated by taking a line integral of the change in the arguments of the in-plane polarization vectors over a closed path 54 .As can be seen in Fig. 5a, b, we successfully identified the positions of vortices (winding number of +1), convergent hedgehog-type structures (winding number of +1), and antivortices (winding number of -1), along with their 3D distributions, which are consistent with the qualitative findings described above.Second, we continued to calculate the topological invariants defined for 3D vectors in a 2D plane, starting from the Pontryagin charge density maps, again for the sliced polarization distributions along the [100] direction (Methods).Figure 5c, d shows the presence of clearly finite Pontryagin charge density near the vortex and antivortex areas, indicating that each vortex and antivortex possesses a non-trivial polar topology.To achieve the complete topological classification, we further determined the skyrmion numbers for each vortex and antivortex by taking a local 2D surface integral of the Pontryagin charge density 55,56 , as shown in Supplementary Fig. 16.In the case of Particle 1, the skyrmion number was found to be approximately +0.3 near the vortex core region, which can be classified as a Bloch-type fractional skyrmion.For Particle 2, as observed in the first representative slice of Supplementary Fig. 16b, the skyrmion numbers of the two vortices were +0.7 and +0.1, respectively, again characterizing them as Bloch-type fractional skyrmions.The skyrmion number of the antivortex, only found in Particle 2, was also identified as -0.4.Additionally, in the case of Particle 2, we confirmed the presence of a polar topology classified as a convergent hedgehog-type fractional skyrmion with a skyrmion number of +0.3 (Supplementary Fig. 16b).The topological analyses described above require 2D slicing of the 3D polarization distributions, and the results can vary depending on the choice of slicing direction.Moreover, due to the finite size and asymmetric geometry of the system, the integration ranges for skyrmion number calculations are limited, and occasionally, two or more neighboring topological objects coexist within the integration range.The electrostatic boundary conditions of our experimental system can also be non-ideal, exhibiting partial screening.The observation of topological charge values deviating from integer or half-integer values, which are generally not allowed in a 2D continuum model, can be attributed to these issues.Note that the emergence of such fractional topological charges under non-ideal conditions was found in previous experiment 57 .Furthermore, we observed that the skyrmion number undergoes an abrupt sign change (from -0.5 to +0.6) through the region marked with black squares in Supplementary Fig. 16b.This typically indicates the existence of a Bloch point, where the polarization vanishes 58 .The 3D polarization distribution near the Bloch point is shown in the Supplementary Fig. 16c, and the topological charge for the surrounding area is +1, indicating that the Bloch point structure has a diverging spiral configuration (Methods).
Finally, since we have obtained the full 3D polarization distributions, the Hopf invariant 54 , indicative of the integer topological charge of a hopfion (3D counterparts of skyrmions), can be calculated.Particles 1 and 2 showed the Hopf invariant of -0.002 and -0.011, respectively, indicating that neither of the particles is likely to possess a hopfion.Through these analyses, for the first time, we were able to observe topologically non-trivial 3D polar structures at the atomic scale and quantitatively calculate their topological invariants to classify them into different categories including fractional skyrmions and spiral Bloch point structures.
The current implementation of AET can only resolve cation atoms with relatively large electron scattering cross-sections.However, to fully understand the surface boundary condition, as well as the resulting 3D polarization distribution and underlying physical mechanism, precise information regarding the 3D locations of low-Z surface adatoms (such as carbon) and internal oxygen atoms is required.With the development of the four-dimensional STEM (4D-STEM) scheme based on highspeed pixelated electron detector technology, it is now possible to fully utilize the information contained in scattered electron beams during STEM measurement 59,60 .Several experimental and computational studies suggest that the 3D atomic coordinates of both low-Z and high-Z elements can be accurately measured together via electron tomography based on 4D-STEM [61][62][63] .We foresee that the full 3D atomic structures (including oxygens and other low-Z elements on the surface) and underlying mechanisms of ferroelectric orderings in nanosized systems will be deciphered in the near future.Additionally, a 4D-STEM-based multislice electron tomography algorithm 62 allows measurements to be taken with the sample in the out-of-focus state, overcoming the limitation imposed by the depth of focus on the sample size.Since the multislice-based method can also compensate for the nonlinear effects caused by multiple scatterings, we anticipate that it will allow us to reveal the toroidal orderings in much larger nanoparticles and the transition mechanisms from the vortex states towards the bulk polar states.
Lastly, it has been predicted that storing information as the handedness of the toroidal ordering can substantially enhance the density of stored information by up to five orders of magnitude (due to the greatly suppressed cross-talk between the information carriers) 15 .Moreover, due to the unique topological nature, polar skyrmions are expected to show enhanced stability over time, and their distinct topological protection renders them highly metastable in environments with field polarization with minimal energy dissipation 64 .If combined with the PTMO state we observed, our observation of nontrivial polar structures (fractional skyrmions and a spiral Bloch point structure) in 3D ferroelectric nanostructures without epitaxial strain can open a door towards even higher capacity storage devices, since it allows the doubling of the memory capacity by separately controlling the linear polarity and chirality, which can carry a bit of information 36 .We hope that our approaches will provide a more comprehensive understanding towards the nature of 3D polar orderings in strain-free nanostructures, and open a pathway leading to functional devices along with further development of switching and reading schemes.

Sample preparation
The BaTiO3 nanoparticles were fabricated using the hydrothermal synthesis.For this work, BaTiO3 nanoparticles were donated by Samsung Electro-Mechanics Co. Ltd.The resulting BaTiO3 nanoparticles were dispersed in ethanol with 3 hours of bath sonication and drop-casted onto SiN membranes of 5 nm thickness.To prevent electron beam damage and charge accumulation, an amorphous carbon film of approximately 2 nm thick was coated on the surface of the specimen.The grids were subsequently heat-treated in a vacuum at 150 °C for 10 hours to minimize hydrocarbon contamination during TEM imaging.

Initial characterization
The atomic structure and chemical composition of the nanoparticles were initially characterized using two techniques: EDS mapping with an FEI Titan cubed G2 60-300 double Cs-corrected and monochromated TEM at the KAIST Analysis Center for Research Advancement (KARA), and PXRD with a Rigaku SmartLab at the KARA.The EDS signals were collected by a Super-X EDS detector in ADF-STEM mode under the following conditions: 300 kV acceleration voltage, 52 pA screen current, and 512 × 512 scan points with 36.65 pm step size.For elemental mapping, Ba Lα, Ti Kα, and O Kα lines were used.The ADF image and EDS maps revealed alternating columns of Ba and Ti atoms, consistent with the <001> zone axis of cubic BaTiO3 (Supplementary Fig. 1).A PXRD measurement was conducted for the nanoparticles using the Cu Kα1 source (8.04 keV) in the two-theta angular range of 10° to 80° with a 0.001° angle step utilizing a D/teX Ultra 250 detector.By comparing the peaks of the PXRD pattern with the standard reference card (PDF #01-084-9618), the powder exhibited a longrange cubic phase (i.e., no (002) peak splitting and c/a of unity), but the precision of our (c/a) value was estimated to be 0.006, considering the PXRD peak widths (Supplementary Fig. 1g).We cannot rule out the possibility that our particles have tetragonal distortion within this c/a ratio range.The lattice constant of a five-atom unit cell (pseudocubic cell) was determined to be 4.026 ± 0.001 Å based on the angular positions of (211), ( 220), (300), (310), and (311) diffraction peaks.

STEM data acquisition
Two tomographic tilt-series of the BaTiO3 nanoparticles were acquired using an FEI Titan Themis3 double Cs corrected and monochromated TEM at the Korea Basic Science Institute in Seoul (for Particle 1), and an FEI Titan cubed G2 60-300 double Cs corrected and monochromated TEM at the KARA (for Particle 2), respectively.The images for each tilt-series were obtained using the ADF-STEM mode at 32 different tilt angles, ranging from −73.0° to +71.6° for Particle 1 and −73.8° to +71.6° for Particle 2, respectively (Supplementary Fig. 2).All the images were acquired under the following microscope parameters: acceleration voltage of 300 kV, detector inner and outer semi-angles of 38 mrad and 200 mrad, respectively, convergence beam semi-angle of 18.0 mrad, and beam current of 7 pA.At each tilt angle, three consecutive images of 1024 × 1024 pixels were collected using a pixel dwell time of 3 μs with a pixel size of 0.354 Å for Particle 1 and 0.364 Å for Particle 2. Note that the pixel sizes were calibrated to match the lattice constant determined from the PXRD with the 3D structures obtained from electron tomography analysis.The total electron dose for the entire tilt series was 1.01 × 10 5 e Å -2 for Particle 1 and 9.52 × 10 4 e Å -2 for Particle 2, respectively.Existing studies show that the electron beam exposure of up to 8.78 × 10 5 e Å -2 does not alter the cation displacement of a BaTiO3 nanocluster of approximately 4 × 4 nm 2 size even at an elevated temperature of 473 K 37 .Therefore, no significant structural damage is expected under the electron dose used in this study.Furthermore, the domain formation or switching induced by electric fields resulting from charge accumulation due to the electron beam is also not a significant concern in this study, as it requires at least ten times larger electron dose 65 .
To ensure minimal beam-induced structural changes during the tilt series acquisition, the zero-degree projections were measured three times throughout the experiment (at the beginning, in the middle, and at the end of the tilt series acquisition).Supplementary Fig. 3a-f illustrates that the internal structures observed in the experiments exhibit consistency.While mild beam-induced atomic diffusion can be observed at the surfaces, particularly for the smaller nanoparticle (Particle 1), this phenomenon was expected and considered during our analysis.The signal from the diffused surface atoms is likely to be averaged out during the reconstruction process 30 , and thus its impact on the reconstructed internal structure becomes insignificant.We artificially introduced surface amorphization to a known structure and simulated a tomography experiment on the changing surface structure to check the effect of surface diffusion.This indeed verified that its effect on the reconstructed internal structure is negligible compared to the precision of AET (More details regarding the evaluation of the effect of surface diffusion are described in later sections of Methods below).

Image post-processing and 3D reconstruction
We performed image post-processing for the tilt series using a series of procedures, including drift correction, scan distortion correction, noise reduction 66,67 , background subtraction, and tilt-series alignment, following the methods described in previous works 27,28,30,[68][69][70] .
(I) Drift correction.We applied drift correction to the images collected at each tilt angle.The stage drift during acquisition was estimated from the three consecutively measured ADF-STEM images and subsequently compensated for using an affine transform that inverts the estimated linear drift effect.
(II) Scan distortion correction.We corrected the ADF-STEM scan distortion resulting from the miscalibrations of fast and slow scan directions by measuring the distortion matrix using a single crystal silicon (110) standard sample.The scan distortions were corrected by applying the proper shear and scaling operations corresponding to the distortion matrix in the Fourier domain 71 .
(III) Image denoising.To address the Gaussian-Poisson mixed noise in ADF-STEM images, we first estimated the Gaussian-Poisson noise parameters from a statistical analysis based on the three consecutively acquired images.The images were then denoised using the block-matching and 3D filtering (BM3D) algorithm 66 .To convert Gaussian-Poisson noise into pure Gaussian noise, as required by the BM3D algorithm, we applied the Anscombe transformation 67 and its inverse transformation to the images before and after the BM3D de-noising process with the estimated noise parameters, respectively.
(IV) Background subtraction.After denoising, we defined a mask slightly larger than the boundary of the nanoparticles for each projection image.We determined the background inside the mask by solving the Dirichlet boundary value problem of the discrete Laplace's equation, and subtracted it from the image.
(V) Tilt-series alignment.Each tilt series was then aligned based on the center-of-mass 72 and commonline alignment 27 methods with sub-pixel accuracy.The images were shifted to a consistent center-ofmass position along the direction perpendicular to the tilt axis.Then, the common-line was computed from each image by projecting it onto the tilt axis.The tilt series was subsequently aligned along the direction parallel to the tilt axis by shifting the projections to ensure complete overlap of the common lines.
(VI) 3D reconstruction.After the image post-processing, we reconstructed 3D tomograms from the post-processed tilt series using the GENFIRE algorithm 31 with the following parameters: discrete Fourier transform (DFT) interpolation method, number of iterations of 500, oversampling ratio of 4, and interpolation radius of 0.1 pixels.

Atom tracing and species classification
The 3D atomic coordinates of internal cation atoms in each nanoparticle were determined using the following atom tracing procedures 27,68,73 .
(i) The 3D local maxima positions in the tomograms were identified and sorted based on the peak intensity in descending order.Starting from the highest intensity, a volume of 1.8 × 1.8 × 1.8 Å 3 (5 × 5 × 5 voxels) centered on each local maximum was cropped, and a 3D Gaussian function was fitted for the cropped volume to find the peak position with sub-pixel accuracy.The fitted peak position was added to the traced atom list only if it did not violate the minimum distance constraint of 2 Å compared to any of the fitted positions already in the list.By repeating this process for all the 3D local maxima, a list of potential atom positions was obtained.
(ii) The potential atom positions were classified into three types of atoms (non-atom, Ti, and Ba) using a k-means clustering algorithm described in previous works 27,68,73 .With this classification procedure, we obtained an initially classified atomic model.
(iii) Due to the imperfectness of the experimental images (non-linear effect, noise, misalignment, etc.), the initial atomic structure contained some atoms that were inconsistent with the ABO3 perovskite structure.As the cation atoms (Ti and Ba) of BaTiO3 nanoparticles should exhibit a bcc configuration, bcc lattice fittings were applied to the initially classified atomic models to filter out the atoms not compatible with the perovskite structure.More details regarding the bcc lattice fitting procedures are described in later sections of Methods below.
(iv) The Ba and Ti layers were defined for each atomic plane perpendicular to the <001> direction, based on the initial tracing and classification results.First, histograms were drawn based on the atom positions along the reference directions ([010] direction for Particle 1 and [001] direction for Particle 2), and the atomic layers along the reference directions were defined based on the peak positions in the histograms.The atomic positions were then assigned to corresponding atomic layers of the smallest distance.The layers dominated by Ba atoms were identified as the Ba layers, while the layers with dominant Ti atoms were identified as the Ti layers.We then added the new best-fit position, obtained from the volume size showing the smallest mean squared residuals, to the traced atom list, provided that the minimum distance constraint of 2 Å between neighboring atoms was met.We further applied bcc lattice fittings again to eliminate some of the newly added fit positions that did not conform to the perovskite structure.
(vi) To finalize the 3D atomic structures, we manually removed physically unreasonable atom candidates and added candidates that were missed during the automated process.During this process, we used a minimum distance constraint of 2 Å between the nearest atoms.In total, we manually added and removed 406 and 39 atoms (for Particle 1) and 265 and 28 atoms (for Particle 2), respectively, resulting in the final 3D cationic atomic models of 11,126 atoms (Particle 1) and 17,260 atoms (Particle 2).Note that the manual adjustment of this level has been commonly employed in the field of AET 27,73 .
(vii) After manually adjusting the atomic positions, we classified the final atomic structures into Ba and Ti atoms based on the layer labellings determined in step (iv), producing 3D models of 11,126 atoms (viii) To assess the consistency of the determined atomic structures with the measured data, we calculated the R-factors 27,28,69 by comparing the experimentally acquired tilt series with the simulated tilt series forward-projected from the final 3D atomic models.The averaged R-factors were 0.22 for Particle 1 and 0.17 for Particle 2. Given that R-factor minimization was not carried out during the structural determination and only cation atoms were used for the forward-projection calculation (excluding oxygens), the R-factors we obtained can be seen as indicative of the reliability of our atomic structures.In fact, R-factors in the range of 0.20 to 0.25 are generally considered acceptable in the crystallography community [74][75][76][77] .

Assignment of experimental 3D atomic positions to ideal bcc lattices
Since the cation atoms of BaTiO3 nanoparticle display a bcc configuration, the 3D atomic positions obtained from the tomograms were assigned to ideal bcc lattice sites using the following procedures.(e) For each calculated site, if an atom is found within 25% of the nearest neighbor distance of the bcc lattice (initially 0.87 Å), the component of the lattice vector (i.e., the vector from the origin to the calculated site) along the reference direction was compared with the atomic layer assigned to the found atom.If the position of the atomic layer (relative to the origin) is consistent with the lattice vector component, the candidate atom was assigned to that bcc lattice site.
(f) The nearest and the next-nearest neighbor search was repeated for all the newly assigned bcc lattice sites.This process was repeated until no further atoms could be assigned to the lattice.
(g) New bcc lattice vectors were fitted to the atoms assigned to the lattice using fitting parameters of translation, 3D rotation, and lattice constant parameters, to minimize the error between the measured atomic positions and the corresponding lattice sites of the fitted bcc lattice.
(h) The iterative processes from step (c) to (g) were repeated until there were no further changes in the fitted lattice vectors.After this procedure, 99.63% (Particle 1) and 99.97% (Particle 2) of the target atoms were successfully assigned to the bcc lattices.The root-mean-square deviations (RMSDs) between the assigned atom positions and the fitted bcc lattices were 92.5 pm (Particle 1) and 75.1 pm (Particle 2), respectively.

Local tetragonal lattice fitting
For the analysis of local tetragonality (c/a ratio), we performed a local tetragonal lattice fitting using the following procedures.First, for each Ti atom, we identified eight nearest-neighbor Ba atoms based on the globally fitted bcc lattice configuration (see 'Assignment of experimental 3D atomic positions to ideal bcc lattices' section of Methods).Second, we performed three tetragonal lattice vector fittings to the atoms assigned to the lattice, for [100], [010], and [001] as potential c-axis directions, respectively.This fitting involved parameters for translation, 3D rotation, and scaling, and was done to minimize the error between the measured atomic positions (which includes only the 8 nearest neighbor atoms for each Ti atom) and the corresponding lattice sites of the fitted tetragonal crystal lattice.Among these, the c-vector direction that minimized the errors between the measured atomic positions and the corresponding lattice sites of the fitted tetragonal crystal lattice was determined as the c-axis for that particular local structure.This procedure was applied to all Ti atoms which have eight nearest neighbor atoms.As a result, the average c/a ratio for each nanoparticle was determined to be 1.002 ± 0.002 for Particle 1 and 1.002 ± 0.001 for Particle 2. Note that the standard error propagation, based on the precision of the atomic coordinates, was used to determine the error bars in this process.

Precision analysis using multislice simulation
To evaluate the reliability of our final 3D atomic structures of the two BaTiO 3 nanoparticles, we conducted a precision analysis.From the determined 3D atomic structures for Particle 1 and Particle 2, 32 projection images were calculated using multislice simulation 78,79 at the experimental tilt angles under the following parameters: incident electron energy of 300 keV, convergence semi-angle of 18 mrad, detector inner and outer semi-angles of 38 mrad and 200 mrad, respectively, slice thickness of 2 Å, −775 nm C3 and +378 μm C5 aberration for Particle 1, and 130 nm C3 and 0 μm C5 aberration for Particle 2. Eight frozen phonon configurations were used to consider the thermal vibration effect at room temperature, which corresponds to the random spatial displacement with a standard deviation of 7.4 pm for Ti and 5.6 pm for Ba 35 .To consider the impact of the electron probe size and other incoherent effects, we applied a Gaussian kernel with a standard deviation of 0.46 Å to each multislice simulated image.
We then employed the GENFIRE algorithm 31 to reconstruct the 3D tomograms from the multislicesimulated tilt series by using the following parameters: DFT interpolation method, number of iterations of 500, oversampling ratio of 4, and interpolation radius of 0.3 pixels.
The atomic structures were determined from the 3D tomograms using the same methods described above (see above for more details regarding atom tracing and species classification).To compare the experimentally determined atomic structures with those from multislice simulations, we calculated the distances between atoms in the experimental structure and the multislice-simulated structure.Any pairs of atoms whose distance was less than a specified threshold distance (half of the nearest-neighbor distance of the ideal fitted bcc lattice: 1.7 Å) were considered as common atom pairs.The result showed that 90.2 % of the atoms in Particle 1 and 94.0 % of the atoms in Particle 2 were properly retrieved in the simulated experiments.The RMSDs of all common atom pairs (i.e., the precision of the atomic coordinates; see Refs. 27,30,68) were found to be 38.3 pm (Particle 1) and 34.4 pm (Particle 2), respectively.

Nanoparticle size estimation
To estimate the size of the BaTiO3 nanoparticle, the volume of the nanoparticle was calculated as the product of the number of unit cells (estimated as the number of Ti atoms) and the volume of a pseudocubic unit cell (65.5 Å 3 ).If we assume a spherical shape, the diameters of the nanoparticles can be estimated to be 8.8 nm for Particle 1 and 10.1 nm for Particle 2, respectively.

Ti atomic displacement field and polarization field calculation
The Ti atomic displacements from the centrosymmetric position of each unit cell were calculated as the deviation of each 3D Ti position from the geometric center of the eight neighboring Ba atomic positions.
The neighboring Ba atoms for each Ti atom were determined as the Ba atoms assigned to the nearest neighbor sites around the Ti atom based on the fitted bcc lattices.Furthermore, the calculation of Ti atomic displacement was performed only for Ti atoms that had all eight nearest neighbor Ba atoms.
In order to calculate the local polarization distribution of BaTiO3 38,41 , the Ti atomic displacement of each unit cell and the displacement of the oxygen octahedron must be identified, as shown in equation (1) : Here,  is an empirical constant with a value of 1.89 ± 0.15 (C m -2 ) Å -1 , obtained using the measured spontaneous polarization and the displacement of the Ti atom from the center of the oxygen octahedron in bulk BaTiO3 38 .  and   are the Ti atomic displacement and the displacement of the oxygen octahedron with respect to the centrosymmetric position of each unit cell, respectively.Since current AET methods can only identify cation atoms, an alternative expression for spontaneous polarization 80 using structural information was additionally employed.This expression is given by where V is the volume of a unit cell,  ⊥ * ,  ∥ * , and   * are the Born effective charges of O perpendicular to the Ti-O direction, O parallel to the Ti-O direction, and Ti atom, respectively.The Born effective charges of each case ( ⊥ * = -2.13, ∥ * = -5.75, and   * = 7.29) we used in this study were obtained from Ref. 81 , which were calculated through ab initio theory with a cubic phase BaTiO3 having a lattice constant of 4.00 Å.
By integrating equations ( 1) and ( 2) to express   as   , and using the value of  ,  * , and V, spontaneous polarization can be calculated by simply multiplying the experimentally observed   with the scale factor (): To obtain the Ti atomic displacement field, we first interpolated the Ti atomic displacements (  ) onto a 3D Cartesian grid with 0.25 Å pixel size by applying a 3D Gaussian kernel 30,82 with a standard deviation of 4.03 Å, equivalent to the nearest neighbor distance between Ti atoms.Subsequently, the polarization field was obtained by multiplying the scaling factor () between Ti atomic displacement (  ) and the spontaneous polarization (  ).The standard error propagation, based on the precision of the atomic coordinates and the  value, was used to determine the error bars in this process, as well as in other calculations conducted in this study.

Analysis of swirling characteristics of the polarization vector fields
A vortex is defined by the rotation of polarization around a core, and its rotational behavior can be described by a toroidal moment (), curl of in-plane polarization, and chirality.The toroidal moment 15,83 , which is determined by the cross product of local polarization with their atomic positions, is given by where N is the number of unit cells,   is the atomic position of the i th Ti atom,   is the local polarization of the i th Ti atom, and  is the averaged polarization.
Curl of in-plane polarization (also called as axial current 84 ) is a concept in fluid mechanics that describes the rotational motion of fluid around a common centerline 85 .In the context of polar structures, it is also used to characterize the clockwise or counterclockwise toroidal ordering 23,86,87 .The curl of the in-plane polarization map was obtained by slicing the 3D polarization vector field along the [100] direction and then calculating the x-component of the curl for the normalized in-plane polarization fields.Note that, in this study, the in-plane polarization fields were normalized before the curl of in-plane polarization calculation since the unnormalized polarization approaches zero in the vicinity of the vortex center, resulting in the emergence of a singularity.
Chirality can be determined through the normalized helicity density (Hn), a mathematical concept adapted from fluid dynamics, which is defined as 88 where  is local polarization.

Topological structures analysis
The characteristics of topological structures can be described by topological invariants (e.g., winding number, skyrmion number, and Hopf invariant).Winding number is defined as a topological invariant in the circle-equivalent order-parameter space, and is used to distinguish between vortices (winding number of +1) and antivortices (winding number of -1) in a 2D system.The winding map was calculated by taking a line integral of the change in the arguments of the in-plane polarization vectors over a closed path (square loop with a side length of 3.75 Å) near different positions in the 2D polarization vector fields (2D vectors in a 2D plane) sliced from our 3D polarization distributions along [100] direction 89 .
Skyrmion number (Nsk) is defined as a topological invariant in the sphere-equivalent order-parameter space and is used to characterize the swirling structure of a skyrmion.To obtain skyrmion number, we first calculated the Pontryagin charge density (qx) maps for the 2D polarization vector fields (3D vectors in a 2D plane) sliced from our 3D polarization distributions along the [100] direction (Fig. 5c, d).This expression is given by 55,56 where   is the Levi-Civita symbol, and  � is the unit vector representing the direction of local polarizations.For the topological characterization of 3D polarization vector fields, we further calculated the skyrmion numbers for each vortex and antivortex, which are given by The skyrmion map was determined by performing a local 2D surface integral of the Pontryagin charge density over a circular area with a radius of 10.00 Å for Particle 1 and 13.75 Å for Particle 2, centered around each position in the (100) plane (Supplementary Fig. 16).Note that skyrmion numbers at each vortex and antivortex can vary depending on the chosen integration areas.Therefore, we calculated multiple skyrmion numbers at the locations of vortices, antivortices, and hedgehog-type structures for several different choices of the integration area and checked their trends.These trends clearly indicate that the skyrmion numbers saturate beyond certain integration areas, which were chosen as the area for our analysis given above.To calculate the topological charge (Q) at the Bloch point, the left-hand side of equation ( 7) can be generalized into 3D based on the divergence theorem 90 .
Hopf invariant (NH), representing the integer topological charge for a hopfion (3D counterparts of skyrmions), was calculated using the following equation 54,91 .
Here,   = 1 8    � • (   � ×    � ) and  is defined as a gauge field satisfying  = ∇ × .The vector potential A was determined in momentum space using the Coulomb gauge condition (∇ •  = 0), and based on this, the Hopf invariant was calculated in momentum space 91 .

Evaluation of the effect of possible surface diffusion during tilt series measurement via simulation
By comparing the experimental projections at zero-degree before and after tomography measurement (Supplementary Fig. 3a-c), we could observe mild structural changes near the surface, especially for Particle 1.To evaluate the effect of surface diffusion for the determination of internal polarization structures, we simulated 3D reconstruction considering the surface diffusion effects and calculated the precision (RMSD) of the retrieved internal structure.
From the determined 3D cation atomic structure of the Particle 1 (total 11,126 atoms), 5,452 atoms were selected as surface atoms by repeatedly applying the Alpha-shape algorithm 92 with shrink factor 1.
To simulate the effect of surface diffusion due to the electron beam damage, the position of each surface atom was randomly perturbed, resulting in the RMSD of 300 pm when we compare the structure before and after the random perturbation (Supplementary Fig. 3g-i).A tilt series of 32 projections was generated by linearly projecting a 3D potential volume from the determined cation atomic structure of Particle 1 (we named this Tilt-Series S1) along the experimental tilt angles.Another tilt series was similarly generated, but for this tilt series, half of the projections were randomly chosen and replaced with the projections obtained from the structure perturbed by the random surface displacement (we named this Tilt-Series S2).
Two 3D reconstructions were obtained from the tilt series (Tilt-Series S1 and S2) using the GENFIRE algorithm 31 with the following parameters: DFT interpolation method, number of iterations of 500, oversampling ratio of 4, and interpolation radius of 0.3 pixels.Using the same atom tracing method, two traced 3D atomic models (model S1 and S2) were obtained from the two 3D reconstruction volumes, respectively.Although the structures of the surface atoms were very different between the two models (due to the perturbation of 300 pm RMSD we applied), the structures of the internal atoms (i.e., nonsurface atoms) were well-preserved and the RMSDs of the atom pairs between the internal atoms in the experimental 3D atomic structure and the traced atomic models (model S1 and S2) were calculated to be 27.3 pm and 39.0 pm, respectively, comparable to the precision our technique can provide.
We subsequently performed the bcc lattice fitting and Ti atomic displacement field analysis (using the same methods described in the above sections) for the internal atoms (i.e., non-surface atoms) of the traced atomic models (model S1 and S2).As can be seen in the Supplementary Fig. 3j-l, consistent internal vortex structures could be obtained even from the tomographic reconstruction which suffered from the surface diffusion effect.

Effective Hamiltonian calculations
The effective Hamiltonian method described in Ref. 83 was employed to calculate the local polarization distributions of ellipsoidal BaTiO3 nanoparticles whose sizes are similar to the experimentally measured ones.In this simulation, the coefficient in front of the square of the local mode within the self-mode energy was slightly changed to place the computational critical temperature of the tetragonal-toorthorhombic phase transition of bulk BaTiO3 to the experimental value of 278 K 93 .This ensures that the room-temperature equilibrium phase of the BaTiO3 bulk is tetragonal.This Hamiltonian has three degrees of freedom, namely the local mode at each 5-atom site (which is directly proportional to the local electric dipoles at these sites), and the homogeneous and inhomogeneous strains.Two different sizes were considered in line with our experiments.The first particle (Particle 1 � ) has the dimensions of (23, 22, 21) in terms of primitive 5-atom unit cells along the x, y, and z-directions, respectively (chosen to be along the [100], [010], and [001] pseudocubic directions).The second particle (Particle 2 � ) is slightly larger, with dimensions of (23, 24, 26).The pseudocubic lattice parameter is 4.01 Å for both Particle 1 � and Particle 2 � at 300 K.The corresponding diameter of each particle used in the simulation was estimated by taking the geometric mean of the dimensions along each axis, resulting in values of 8.8 nm (Particle 1 � ) and 9.7 nm (Particle 2 � ), respectively.
During the simulation, the samples were gradually cooled from 600 K within a Monte-Carlo procedure under mechanically and electrically free boundary conditions, while the 3D local polarization distributions were obtained for each temperature.The toroidal moment 83 given by  =  *  2 ∑   × (  − 〈〉)  (10)   were also calculated at each temperature, where the dipole vortex has started to form around 300 K (Supplementary Fig. 10a, b).Here,  * is the dynamical charge,  is the electron charge, a is the lattice constant for the 5-atom pseudocubic unit cell,  is the volume of the nanodot,   is the position vector for the i th unit cell,   is the local mode of the i th unit cell, and 〈〉 is the averaged local mode over the entire particle.

Figure 1 |
Figure 1 | Determination of the 3D atomic structure of a 10.1 nm BaTiO3 nanoparticle (Particle 2).a, Cation atomic layers along the [001] crystallographic direction, showing the 3D atomic arrangements of alternating Ba and Ti layers.b-d, 1.07 Å thick internal slices of the 3D tomogram along the [001] direction for three consecutive atomic layers near the core of the particle.The intensity of the sliced tomogram is depicted in grayscale, while the atomic coordinates of Ba and Ti are represented with blue and red dots, respectively.Scale bar, 2 nm.e-g, Linear projections of a 3D intensity volume representing the determined 3D atomic structure, projected along [100] (e), [010] (f), and [001] (g) directions, respectively.The grayscale background indicates the projected intensity, where Ba and Ti columns can be distinguished by their relative intensities (Ti columns show weaker contrasts).The arrows represent the Ti displacements, calculated as the deviation of the identified Ti column positions from the geometric centers of the four neighboring Ba column positions.The arrows are colored based on the displacement direction, as given in the color wheel in (e).

Figure 2 |
Figure 2 | Sliced maps showing the 3D distribution of Ti atomic displacements for the 8.8 nm BaTiO3 (Particle 1).a, In-plane atomic displacement direction maps of the representative Ti atomic layers sliced along the [100] direction.Note that the spacing between the plotted layers is 2 unit cells.The displacements were calculated as the deviation of each 3D Ti position from the geometric center of the eight neighboring Ba atomic positions and subsequently averaged using a Gaussian kernel (Methods), and their in-plane directions are plotted (the color of arrows also indicates the direction as given in the color wheel).The existence of a single counterclockwise vortex can be clearly seen (marked with blue arrows).Scale bar, 2 nm.b, Maps showing the magnitude of the kernel-averaged 3D Ti atomic displacements and local polarization for the corresponding Ti atomic layers in (a).The local polarizations were obtained through a linear relation with 3D Ti atomic displacements.The in-plane directions of the displacements are overlaid as black arrows, and the positions of counterclockwise vortices given in (a) are consistently marked with blue arrows.

Figure 3 |
Figure 3 | Sliced maps showing the 3D distribution of Ti atomic displacements for the 10.1 nm BaTiO3 (Particle 2).a, In-plane atomic displacement direction maps of the representative Ti atomic layers sliced along the [100] direction.Note that the spacing between the plotted layers is 2 unit cells.The displacements were calculated as the deviation of each 3D Ti position from the geometric center of the eight neighboring Ba atomic positions and subsequently averaged using a Gaussian kernel (Methods), and their in-plane directions are plotted (the color of arrows also indicates the direction as given in the color wheel).The existence of multiple counterclockwise vortices can be clearly seen (marked with blue arrows).Scale bar, 2 nm.b, Maps showing the magnitude of the kernel-averaged 3D Ti atomic displacements and local polarization for the corresponding Ti atomic layers in (a).The local polarizations were obtained through a linear relation with 3D Ti atomic displacements.The in-plane directions of the displacements are overlaid as black arrows, and the positions of counterclockwise vortices given in (a) are consistently marked with blue arrows.

Figure 4 |
Figure 4 | Simulated local dipole patterns in BaTiO3 nanoparticles at 300 K. a, b, Simulated local dipolar structure maps at 300 K, representing the slices perpendicular to the [100] direction for the smaller (Particle 1 � , 8.8 nm) (a), and larger (Particle 2 � , 9.7 nm) (b) BaTiO3 nanoparticles.The existence of a single vortex and multiple vortices (marked with red arrows) can be clearly observed from the smaller and larger nanoparticles, respectively, consistent with the experiments.The lattice constant of 5-atom unit cell was 4.01 Å in the simulation.

Figure 5 |
Figure 5 | Representative 2D slices through the nanoparticles showing topological invariants.a, b, Winding number maps of representative Ti atomic layers along the [100] direction of Particle 1 (8.8 nm) (a), and Particle 2 (10.1 nm) (b).The in-plane directions of the Ti displacements are indicated by black arrows.The red, yellow, and blue colors represent the vortex (winding number of +1), hedgehog-type structure (winding number of +1), and antivortex (winding number of -1), respectively.Note that the winding number itself cannot distinguish between stripe domains, flux-closure domains, hedgehog-type structures, and vortex (or antivortex) structures; the identification of vortex structure requires more information regarding the local polarization distributions.Accordingly, we marked in green the areas where the winding number is +1 or -1, but which are neither vortices, hedgehog-type structures, nor antivortices.The distance between the arrows is 3.75 Å. c, d, Pontryagin charge density maps of representative Ti atomic layers along the [100] direction of Particle 1 (8.8 nm) (c), and Particle 2 (10.1 nm) (d).The in-plane directions of the Ti atomic displacements are indicated by overlaid black arrows.Vortices, antivortices, and hedgehog-type structures are marked with red, blue, and yellow circles, respectively.The distance between the arrows is 3.75 Å.
(v) Due to the intensity elongation effect of the tomogram along the missing-wedge direction and slight imperfections caused by the data imperfectness described above, several connected intensity blobs appear in the tomogram.Since the 3D local maxima of the connected intensity blobs are not welldefined, an additional atom tracing procedure was necessary to fully determine the atomic structures.The 3D tomogram was sliced [slice thickness of 1.1 Å (3 pixels)] for every atomic layer in either the [010] direction (for Particle 1) or the [001] direction (for Particle 2), and each slice was summed along the slice directions to obtain 2D images.For all the 2D slices, the positions and intensities of 2D local maxima were extracted.Starting from the highest intensity, the 3D Gaussian fitting procedure was repeated using the 2D local maxima positions as initial positions.The size of the cropped volume for the 3D Gaussian fitting was adjusted by changing the side length, ranging from three to seven pixels.

with 5 ,
497 Ti and 5,629 Ba atoms for Particle 1 and 17,260 atoms with 8,268 Ti and 8,992 Ba atoms for Particle 2, respectively.
(a) First, histograms were drawn based on the atom positions along the reference directions ([010] direction for Particle 1 and [001] direction for Particle 2), and the atomic layers along the reference directions were defined based on the peak positions in the histograms.(b) Each atom was classified into the atomic layer with the smallest distance to it.(c) A Ba atom closest to the mean position of the 3D atomic coordinates was chosen as the origin of an ideal bcc lattice.(d) Next, the nearest and next-nearest bcc sites of this atom were calculated based on the initial pseudocubic lattice vectors with a lattice constant of 4.036 Å (the bulk lattice constant of cubic BaTiO3).

Supplementary Figure 1 |
Characterization of a BaTiO3 nanoparticle.a, An ADF-STEM image of a representative BaTiO3 nanoparticle viewed along the <001> direction.Scale bar, 5 nm.b-d, The EDS intensity maps after Gaussian smoothing (standard deviations of 2 pixels) for Ba (b), Ti (c), and O (d). e, The ADF-STEM line scan profile within the region marked with a green box in (a); the alternating Ba and Ti columns can be clearly identified (weaker peaks represent the Ti columns).f, The atomic fraction line profile of Ba (blue) and Ti (red) within the region marked with the same green box in (a), confirming that the peaks in (e) (i.e., the atom columns) are indeed the alternating Ba and Ti columns.g, PXRD pattern of the BaTiO3 nanoparticles.The five peaks used for the lattice constant determination are indexed based on the pseudocubic convention in the figure.The inset shows the PXRD pattern in the two-theta angular range of 44° to 46° corresponding to the (002) diffraction peak.The black dots represent the measured PXRD data and the blue solid lines indicate the measured data after Gaussian smoothing.