Creating polar antivortex in PbTiO3/SrTiO3 superlattice.

Nontrivial topological structures offer a rich playground in condensed matters and promise alternative device configurations for post-Moore electronics. While recently a number of polar topologies have been discovered in confined ferroelectric PbTiO3 within artificially engineered PbTiO3/SrTiO3 superlattices, little attention was paid to possible topological polar structures in SrTiO3. Here we successfully create previously unrealized polar antivortices within the SrTiO3 of PbTiO3/SrTiO3 superlattices, accomplished by carefully engineering their thicknesses guided by phase-field simulation. Field- and thermal-induced Kosterlitz-Thouless-like topological phase transitions have also been demonstrated, and it was discovered that the driving force for antivortex formation is electrostatic instead of elastic. This work completes an important missing link in polar topologies, expands the reaches of topological structures, and offers insight into searching and manipulating polar textures.

Introduction

Both spins and dipoles prefer alignment and often form uniform patterns that are topologically trivial. Nontrivial topologies such as vortices may arise[1,2], as schematically shown in the Supplementary Fig. 1a, often resulted from delicate energetic balance in confined structures and leading to exotic properties[3-6]. Such topological structures play important roles in condensed matter physics including fluid dynamics[7,8], superconductivity[9], and ferromagnetism[10,11], and they promise alternative device configurations for post-Moore spintronics and electronics[12,13]. Indeed, magnetic skyrmions are actively pursued for high-density data storage[14], while polar vortices with exotic negative capacitance[4] may enable ultralow power consumption in microelectronics. Following extensive investigations on a variety of magnetic textures including vortices[5,11], domain walls[15] and skyrmions[16] in the past decades, studies on polar topologies have taken off in recent years, resulting in discoveries of closure domains[17], vortices[18], skyrmions[19], and meron[20] in ferroelectric materials.

However, creating topologies in polar systems is usually more difficult as the dielectric anisotropy in polar materials is much stronger than magnetic ones[1], and there is tremendous energy penalty when polarization rotates to form polar topologies. As a result, the atomic-scale polar textures such as closure domains[17,21], vortices[2,18,22], and skyrmions[19] have only been observed in confined PbTiO3 (PTO) layers within appropriately designed (PTO)/(STO) superlattice (n unit cell (u.c.) thick PTO and m u.c. thick SrTiO3 (STO)). These works mainly focused on PTO within PTO/STO superlattices, which has large polarization, while little attention was paid to the STO. In fact, the cubic structure and paraelectric phase of STO are very delicate to external disturbance[23]. Therefore, it is insightful to examine what happens to the nominally paraelectric STO sandwiched between polar PTO layers. Equally important is if other topological structures often observed in magnetism, such as vortex–antivortex pairs as schematically shown in Supplementary Fig. 1, exist in such dielectric superlattice system or not. Following the groundbreaking work of Kosterlitz and Thouless, it is now well known that vortex–antivortex pair as schematically shown in Fig. 1a may form during Kosterlitz–Thouless transition[24], which substantially reduce the energy penalty arising from individual vortex and antivortex. Such a vortex–antivortex pair has indeed been observed in superconducting[9] as well as ferromagnetic systems[25,26], and it requires a pair of vortices with identical orientation. In (PTO)/(STO) superlattice, however, the neighboring vortices in each PTO layer are observed to possess opposite orientations[18,27], making the topology in between trivial. We thus turn our search for antivortex to STO sandwiched between two layers of PTO instead. This may appear counterintuitive at the first sight, though polar order has indeed been observed in STO before at the reduced thickness[23] or in a confined heterostrucure[28,29], and the weaker polarization induced in nominally cubic STO may exhibit weaker anisotropy, facilitating the formation of vortex–antivortex pair.

Fig. 1

Designing vortex–antivortex pair in (PTO)/(STO) superlattice.

a Schematic illustration of a topological antivortex sandwiched between two adjacent vortices. b–e Four typical polar structures exist in (PTO)/(STO) superlattices for different combinations of m and n, as predicted by phase-field simulation. b For 4-u.c. thick STO (m = 4) sandwiched between two 4-u.c. thick PTO (n = 4), antiparallel a-domain is observed in PTO, while polarization in STO is negligibly small, exhibiting no nontrivial topological structure. c For m = 20 and n = 10, vortex array emerges in PTO, while polarization in STO remains negligibly small. d For m = 10 and n = 10, sign of topological structure appears in STO, with modestly increased polarization, while antivortex appears irregular. e For m = 4 and n = 10, perfect antivortex array with relatively large polarization is observed in STO, sandwiched between two vortices in adjacent PTO. f Phase-field computed phase diagram of four typical polar structures in (PTO)/(STO) superlattices, as represented by b–e.

Motivated by such considerations and guided by detailed phase-field simulations, we design a series of (PTO)/(STO) heterostructures, and successfully create vortex–antivortex pairs in (PTO)10/(STO)4 system, where the vortices exist in the PTO layers and the antivortices exist in the STO layers. To our best knowledge, this is not only an observation of atomic-scale vortex–antivortex pair in a dielectric system, but also a realization of nontrivial polar topology in STO when embedded in a ferroelectric superlattice system. Furthermore, topological phase transition can be induced by either temperature change or electric field, and by examining the energetics of the superlattice, we conclude that the driving force for such antivortex formation is electrostatic, while misfit strain plays a negligible role. The successful creation of previously unrealized atomic-scale vortex–antivortex pair in PTO/STO superlattices expands the reaches of topological structures and completes an important missing link in polar topologies. The work thus sheds considerable insight into the formation of topological polar structures and offers guidance in searching for polar textures.

Result

Phase-field simulation

We first seek to create an antivortex as schematically shown in Fig. 1a, sandwiched between two vortices. Such topological structure has been predicted by Mermin from energetic point of view[30], though its realization has yet to be demonstrated in a dielectric system. We thus consider superlattices with configuration of (PTO)/(STO), wherein m-u.c. thick STO layer is sandwiched between two n-u.c. thick PTO layers. Array of polar vortices has recently been observed in PTO layer of such superlattices[18], giving us hope that under appropriate design antivortex may emerge in STO sandwiched between two vortices in two adjacent PTO layers. Based on systematic phase-field simulations, we have identified four typical polar configurations (Fig. 1b–e) for (PTO)/(STO), enabling us to construct a phase diagram to guide the design (Fig. 1f). When both PTO and STO layers are ultrathin, for example, for m, n = 4, a-domain is observed in two adjacent PTO layers, while STO layer is also slightly polarized (Supplementary Fig. 2), exhibiting no nontrivial topology (Fig. 1b). When they are both relatively thick, for example, n = 10 and m = 20, nontrivial vortex array emerges in PTO, while polarization in STO remains negligibly small (Fig. 1c). We thus keep n = 10 to maintain the desired vortex array in PTO, and reduce the thickness of STO. At m = 10, sign of antivortex pattern appears in STO (Fig. 1d), with its polarization magnitude modestly increases, though the topological structure is not regular, and two vortices in adjacent PTO is not well aligned. When thickness of STO is further reduced to 4 u.c., regular antivortex emerges in STO (Fig. 1e), sandwiched between two nicely aligned vortices in PTO, and the magnitude of its polarization increases further as well. This is precisely what we are looking for, fully consistent with theoretical expectation illustrated in Fig. 1a. The window for the vortex–antivortex pair is quite narrow in the phase diagram (Fig. 1f), with thickness of PTO ranging between 8 and 12 u.c. and thickness of STO smaller than 8 u.c. Note that similar superlattices have been studied by Hong et al. [29], though their focus was polar configuration of PTO.

(PTO)10/(STO) superlattice design

Encouraged by phase-field simulation, we design a gradient superlattice heterostructure of (PTO)10/(STO), with thickness of PTO fixed at 10 u.c., while that of STO varying at 4, 7, 10, and 15, as shown in Supplementary Fig. 3a. The superlattice heterostructures were then grown on DyScO3 (110) substrate by pulsed-laser deposition (PLD). The low-magnification high angle annular dark field (HAADF) scanning transmission electron microscopy (STEM) image of (PTO)10/(STO) heterostructure in Fig. 2a illustrates the stacked ferroelectric and dielectric layers of different contrast, whereas Supplementary Fig. 3b–d shows distributions of strain components estimated from geometric phase analysis based on the STEM image. The white-colored sinusoidal wave-like out-of-plane strain pattern is observed within PTO layers along the [100] direction, suggesting the existence of long-range vortex ordering consistent with previous reports[21,22,31]. The dark field transmission electron microscopy (TEM) image shown in Fig. 2b depicts periodic array of bright and dark intensity modulation, corresponding to the clockwise–anticlockwise vortex pairs previously reported in PTO layers[2,18,31]. We can also see such vortex ordering from the spatial distribution of polarization calculated from phase-field simulation (Fig. 2c) that closely resembles Fig. 2b, wherein zoomed-in examination at the interface between dark and bright contrasts clearly reveals a polar vortex.

Fig. 2

Polar vortex–antivortex pairs in designed (PTO)/(STO) superlattice.

a A low-magnification HAADF image depicts STO layers with varying thicknesses (4, 7, 10, and 15 u.c.) sandwiched between 10-u.c. PTO layers. b Dark field TEM image under two-beam conditions by selecting (002)pc g vector (subscript pc denotes pseudocubic). The periodic array of bright and dark intensity modulation corresponds to vortex arrays within PTO layers. c The spatial distribution of the out-of-plane polarization (unit: μC/cm2) was calculated from phase-field simulation. Inset: enlarged view of the polar vector configuration (black arrows). d An atomically resolved HAADF image for a 4-u.c. thick STO sandwiched between adjacent 10-u.c. PTO layers, colored for clarity. e Map of polar vectors between cations extracted from the HAADF image depicts vortex–antivortex texture. The cores of antivortices are highlighted by the dotted diamond boxes. Enlarged views of polar vectors overlaid with polar angle variation taken from dashed highlighted rectangle boxes in e for antivortex (red color) in STO (f) and vortex (black color) in PTO (g), respectively. Variation of polar displacement components within antivortex structure along A–B (h) and C–D (i) directions, as marked in e. D represents in-plane displacement and D represents out-of-plane displacement.

Atomic-scale antivortex structure

In order to confirm the polar structure in the superlattice at the atomic scale, we acquired high-magnification HAADF image for 4-u.c. STO sandwiched between 10-u.c. PTO, as shown in Fig. 2d. The Z-contrast sensitivity of HAADF imaging shows sharp and coherent interfaces between PTO and STO (Z is the atomic number), which is also confirmed by the atomically resolved energy dispersive X-ray spectra mapping incorporated in Supplementary Fig. 4. The polar map (Fig. 2e) of displacement vectors between A site (Pb, Sr) and B site (Ti) derived from HAADF image[32] illustrates a pair of antivortices within the STO layer, as highlighted by the dotted diamond boxes at their cores, and each antivortex is sandwiched between a pair of vortices in adjacent PTO layers, fully consistent with the theoretical expectation in Fig. 1a. To better appreciate the topology of polar structures, enlarged views of polarization vectors overlaid with polar angle variation for the marked rectangular boxes in STO (Fig. 2f) and PTO (Fig. 2g) are examined, revealing clearly antivortex structure in STO and vortex in PTO. Moreover, the variations of polar displacement within the antivortex along A–B (Fig. 2h) and C–D (Fig. 2i) directions show that out-of-plane (D) and in-plane (D) polar vectors reverse their directions when passing through the antivortex core, approaching and departing the core from two sets of opposite directions (head-to-head and tail-to-tail)[26,28,30,33]. Using the experimental data, we also obtain the distribution of winding numbers[34] (Supplementary Fig. 5a), confirming the topological nature of vortices and antivortices about their respective cores, where two antivortices with winding number −1 between four vortices with winding number 1 are revealed. Additional details on the polar topologies can be found in Supplementary Fig. 5b–d along with phase-field simulations (Supplementary Fig. 5e–g), which show good agreement between experiment and simulation.

The accurate quantitative measurement of polarization in STO, particularly at larger thicknesses, remains a challenge for HAADF image because the polarity in STO mainly arises from the displacement of the oxygen[23], while HAADF tends to underestimate the STO polarization relative to that of PTO[35]. Thus, we also acquired integrated differential phase contrast (iDPC) image[22], which presents the information of oxygen configurations with picometer precision and thus gives better accuracy for polarization measurements based on the atomic displacements between cations and oxygen[22] (see “Methods” for details). From the iDPC image (colored for clarity) in Fig. 3a, the atomic shift between Sr and O for 4-u.c. thick STO with respect to their respective centrosymmetric positions is up to ~20 pm (Supplementary Figs. 6 and 7), visible even with the naked eye. From the enlarged views of atomic structure shown in the inset of Fig. 3a, the octahedron shift in STO (red) are similar to that of PTO (shallow yellow) except less pronounced. The corresponding polar map in Fig. 3b illustrates three antivortices in the STO layer, with their cores highlighted by the dotted boxes. When passing through one of the antivortex cores along either A–B or C–D direction as marked, the polar vectors reverse their directions (Fig. 3c, d) in a similar manner as already revealed by HAADF image (Fig. 2h, i), demonstrating high fidelity of our analysis based on two independent techniques and data sets. The atomic structure of the polarized STO exhibits larger displacement between cations and O and smaller one between Sr and Ti, similar to that of PTO[36], though its magnitudes of displacements and thus polarization are much smaller (see Supplementary Fig. 7). Indeed, the polarization in the polarized STO is mainly contributed by the oxygen displacement, in good agreement with the previous study[23]. Here, the average magnitude of the polarization is estimated to be ~30 μC/cm2 (Supplementary Fig. 7g), consistent with phase-field simulation as well as previous first principle density functional theory calculations[28]. With increased STO thickness, the polarization decreases (Fig. 3e and Supplementary Fig. 8), with the antivortex-like polar topology remaining for 7-u.c. STO, though the structure is less ideal (Supplementary Fig. 9).

Fig. 3

Detailed polarization distribution of polar vortex–antivortex pairs in (PTO)10/(STO)4.

a An atomically resolved iDPC image for 4-u.c. thick STO sandwiched between two 10-u.c. PTO layers, colored for clarity. From the enlarged views taken from the marked regions within PTO and STO, the atomic shift between cations and oxygen is visible with the naked eye. b The corresponding unit-cell scale map of polarization vectors, calculated from the atomic displacements between cations and oxygen. Arrows denote the polarization orientation and the color represents the magnitude. The yellow dotted boxes highlight the locations of antivortex cores. Variation of polarization along A–B (c) and C–D (d) directions, as marked in b. P represents in-plane polarization and P represents out-of-plane polarization. e The comparison of measured (orange) and phase-field simulated (blue) average polarization versus m for (PTO)10/(STO). The error bar represents the standard deviation.

Formation mechanism of polar antivortex

These two sets of independent STEM data acquired using HAADF and iDPC techniques unambiguously established the existence of antivortex topology in STO, and it is compelling for us to examine its energetics, as shown in Fig. 4a, so that we can understand its formation mechanism. Interestingly, the electrostatic energy density in STO is found from phase-field simulation to be negative and decrease with reduced STO thickness, while the corresponding elastic energy is positive and does not change much with STO thickness. It suggests that the formation of antivortex in STO is largely driven by electric field, while misfit strain in superlattice plays a negligible role. This is in sharp contrast to corresponding analysis for PTO (Supplementary Fig. 10) showing that elastic energy is negative while electric energy is positive, so that the driving force for vortex formation in PTO is elastic, as commonly understood. An immediate implication of this finding is that we may be able to tune the antivortex electrically[37], as exhibited by the hysteresis loop of winding number versus electric field in Fig. 4b. Reversible field induced topological phase transition is observed, where the vortex–antivortex pair is turned into a single-domain state upon a modest electric field around 392 kV/cm, and is recovered when the electric field is reduced to 209 kV/cm. Similar phase transition can also be induced by heating and cooling (Supplementary Fig. 11), analog to Kosterlitz–Thouless transition[24]. The electric tuning also enables us to examine the stability of antivortex in STO with different thickness, measured by the critical electric field for topological phase transition (Fig. 4a). As expected, this critical field is largest for 4-u.c. STO, demonstrating its highest stability, while that of 10-u.c. STO is substantially reduced and thus is much less stable. Furthermore, the antivortex exhibits positive and much enhanced capacitance at its core, as shown in Fig. 4c, while vortex possesses negative capacitance in excellent agreement with previous report[4] (Supplementary Fig. 12). Interestingly, the field induced phase transition renders dielectric hysteresis as shown in Fig. 4d, where dielectric tunability as large as 50.7% is observed.

Fig. 4

Formation mechanism of antivortex in STO and its electric tuning via phase-field simulation.

a The electrostatic and elastic energy densities (on the left) for (PTO)10/(STO) heterostructures and the critical electric fields (on the right) under which the antivortex disappears. b The hysteresis loop of winding number of STO within (PTO)10/(STO)6 versus external electric field. The left and right insets show the polarization distributions of vortex–antivortex pair and single-domain state, respectively. c The spatial distribution of local permittivity in the middle plane of STO across the antivortex cores. The two peaks indicate the significant increase of the permittivity at the antivortex cores. d The average permittivity in the middle plane of STO versus the external electric field. The two abrupt changes of permittivity in the hysteresis loop are induced by the topological phase transition.

Discussion

In summary, our work completes an important missing link in polar topology, where stable antivortex is finally confirmed to exist at atomic scale after the discoveries of flux closure, polar vortex, and skyrmions in the artificial PTO/STO superlattices. The small and highly nonuniform polarization may explain why all previous studies did not notice the existence of antivortex in STO. More importantly, it points toward a direction for designing polar topologies in artificial dielectric systems. The energetics of such vortex–antivortex pair is rather delicate, making it nontrivial to realize experimentally, but also make it easier to manipulate via external mechanisms that can easily tip the energetic balances, as we have demonstrated. Our study thus offers a realistic roadmap forward to ultimately engineer and control the polar topologies for devices applications.

Methods

Fabrication of designed gradient heterostructures

(PTO)/(STO) superlattice heterostructures were grown via PLD (PLD-5000) equipped with a KrF excimer laser (λ = 248 nm)[22]. Heterostructures were deposited on (110)-DyScO3 substrates by alternately ablating ceramic targets of SrTiO3 and Pb1.1TiO3 at a laser energy of 340 mJ pulse−1, and a laser repetition rate of 10 Hz. The substrate was heated to 600 °C in a dynamic oxygen pressure of 200 mtorr for the growth of the PTO and STO layers. Thicknesses of the PTO and STO layers were held at desired thickness through controlling the number of laser pulse. Following growth, the (PTO)/(STO) superlattice heterostructures were cooled to room temperature in 200 mtorr oxygen pressure at 10 °C min−1.

TEM cross-sectional sample preparation

For image acquisition, the cross-sectional TEM specimen was thinned to less than 30 μm first by using mechanical polishing. The subsequent argon ion milling was carried out using PIPSTM (Model 691, Gatan Inc.) with the accelerating voltage of 3.5 kV until a hole was made. Low voltage milling was performed with accelerating voltage of 0.3 kV to minimize damage and remove the surface amorphous layer.

Electron microscopy characterization and image analysis

The dark field TEM image shown in Fig. 2b was carried out under the two-beam condition with g vector: g = 002pc from an aberration-corrected FEI Titan Themis G2 at 300 kV. HAADF and iDPC images were also recorded at 300 kV using an aberration-corrected FEI Titan Themis G2. The convergence semiangle for imaging is 30 mrad, the collection semiangles snap is 4–21 mrad for the iDPC imaging, and 39–200 mrad for the HAADF. The atom positions were determined by simultaneously fitting with two-dimensional Gaussian peaks using a MATLAB code[32]. The polar vectors in Fig. 2e were plotted from the offset between A site (Pb and Sr) and B site (Ti) sublattices based on the HAADF-STEM image in Fig. 2d. To determine the atomic shift for each atom column in the HAADF, the displacements of A (B) with respect to the center of their surrounding four B (A) columns are measured and decomposed into in-plane and out-of-plane components, respectively. For the iDPC image along [010] direction (Fig. 3a), each cation column [Pb(Sr)] and TiO is surrounded by four oxygen columns. The displacements of cations with respect to the center of their surrounding four oxygen columns can be measured along in-plane and out-of-plane directions, receptively (Supplementary Figs. 6 and 7). Based on the displacements, the unit-scale polarization (Fig. 3b, e and Supplementary Fig. 7g) can be calculated according to P = [38], where V is the volume of u.c. which for our case is a2c, δ is displacement/shift of atom (i) from their centrosymmetric position, and Z is the Born effective charge of atom (i) calculated by ab initio theory, 6.71 for Ti, 3.92 for Pb, −2.56 for O in PTO, and 7.12 for Ti, 2.54 for Sr, −2.00 for O in STO[38]. Taking the oxygen sublattice as the reference (standard position of each u.c.), the polarization is simplified to be P = . The vector and magnitude maps of displacement (Fig. 2e) and polarization (Fig. 3b) are plotted by Origin.

Winding number calculation

We carried out a local winding number analysis for polar angle distribution in order to authenticate the topological nature of vortices and antivortices about their respective cores. The two-dimensional winding number n along a closed loop C was calculated by the following line integral , where ∇θ is the angle gradient of polarization vectors along the integral loop C[30,34,39]. By performing a local winding number calculation on each closed loop, we are able to identify the existence of a single vortex and antivortex that gives a winding number equal to +1 and −1, respectively (for illustration simultaneously observe Figs. 2e, f, g and 4b and Supplementary Figs. 5 and 11). We used a condition of |Δθ| <180° to determine the angle rotation direction. We have accomplished the winding number quantification for each loop by taking 2 × 2 square loops tilted on the vector map sharing the boundaries.

In the phase-field modeling, the spatially dependent polarization vector P is selected as the order parameter to describe the polar states, and the total free-energy density of a PTO/STO superlattice thin film takes the following form:

where α, α, and α are the Landau expansion coefficients (sixth- and fourth-order forms for PTO and STO, respectively), c is the elastic constant, q is the electrostrictive coefficient, g is the gradient energy coefficient, ε0 is the dielectric constant of vacuum, and εr denotes the relative dielectric constant of the background material (cubic PTO and STO in this case). The summation convention for the repeated indices is employed, and the Latin letters i, j, k, and l take 1 and 2 in the present work. The detailed expression of each Landau energy forms can be found in the literature[40]. Based on the total free-energy density, the temporal evolution of the polarization field can be obtained by solving the time-dependent Ginzburg–Landau (TDGL) equation

where L represents the domain wall mobility, is the total free-energy, r is the spatial position vector, and t denotes time. Besides the TDGL equation, both the mechanical equilibrium equation

and the electrostatic equilibrium equation

must be satisfied simultaneously for a ferroelectric system without body force and space charge.

To solve the above equations, the nonlinear finite element method and backward Euler iteration method are employed for space discretization and time integration, respectively. For clearer illustration and computational simplicity, the carried-out simulations were restricted to the [100]-[001] crystallographic plane, which corresponds to the x-z plane in the Cartesian coordinate system. Discrete grids with Δx = Δz = 0.4 nm in real space were used for space discretization, and the step length for time integration was chosen as Δt/t0 = 0.2, where t0 = 1/(α0L) and α0 is the absolute value of α1 at room temperature. Periodic boundary conditions for the electric potential and polarization components were employed along the x direction. The material parameters for PTO and STO used in the simulations are given as follows, based on previous work[40,41]. For PTO, α1 = 3.8(T − 479) × 105 C−2 m2 N, where T is temperature in °C, α11 = −7.3 × 107 C−4 m6 N, α12 = 7.5 × 108 C−4 m6 N, α111 = 2.6 × 108 C−6 m10 N, α112 = 6.1 × 108 C−6 m10 N, α123 = −3.7 × 10 9 C−6 m10 N, Q11 = 0.089 C−2 m4, Q12 = −0.026 C−2 m4, Q44 = 0.03375 C−2 m4, C11 = 1.746 × 1011 N m−2, C12 = 7.937 × 1010 N m−2, C44 = 1.111 × 1011 N m−2, ε = 66, P0 = 0.757 C m−2 is the spontaneous polarization of PTO at room temperature. For STO, α1 = 7.06(T+238) × 105 C−2 m2 N, α11 = 1.70 × 109 C−4 m6 N, α12 = 1.37 × 109 C−4 m6 N, Q11 = 0.0457 C−2 m4, Q12 = −0.0135 C−2 m4, Q44 = 0.00975 C−2 m4, C11 = 3.156 × 1011 N m−2, C12 = 1.01 × 1011 N m−2, C44 = 1.19 × 1011 N m−2, ε = 100. Note that the electrostrictive coefficients Qij were transformed into qij using the formulas provided by ref. [41], and the Landau energy density coefficients were correspondingly modified. In addition, for both PTO and STO, g11/G110 = 0.4, g44/G110 = 0.4, where G110 = 1.73 × 10−10 C−2 m4 N. A normalization process for the material parameters was used to achieve better numerical stability in the simulations, which can be referred to elsewhere[42,43]. The domain structures of suppelattice system are evoluted from initial setup for polarization with small random fluctuation (<0.01P0). The cubic lattice constants for paraelectric PTO and STO were assumed as 3.955 and 3.905 Å, respectively, hence a misfit strain of −0.3% was applied on the PTO layers in order to take into account the resulted interlayer mechanical inhomogeneity.

The two-dimensional phase-field simulations were carried out for (PTO)/(STO) superlattices at room temperature. First, corresponding to the configuration of the gradient (PTO)10/(STO) superlattice in the experiment, the phase-field simulations were used to calculate its polarization and strain distributions. In addition, the phase-field calculations were further conducted to complete a phase diagram for the evolution of antivortex states in (PTO)/(STO) superlattices. The average polarization (blue color data in Fig. 3e) variation within the STO layer for (PTO)10/(STO) superlattice is calculated by , where P is the polarization magnitude of each node inside STO layer and N is the total number of the nodes. For energy density calculation we utilized the formula , where f is the total energy density of each element inside (PTO)10/(STO)/(PTO)10 superlattice system wherein the subscripts m = 4, 6, 8, and 10 u.c as can be seen in Fig. 4a. Furthermore, fi is calculated by Gaussian integration based on the Eq. (1) and M is total number of elements in the investigated (PTO)10/(STO)/(PTO)10 system. The local permittivity given in Fig. 4c and Supplementary Fig. 12 is calculated by ε = ΔE/ΔP after applying a small electric field to the initial stable domain structure without electric field.