Theoretical guidelines to create and tune electric skyrmion bubbles

See allHide authors and affiliations

Science Advances  15 Feb 2019:
Vol. 5, no. 2, eaau7023
DOI: 10.1126/sciadv.aau7023


Researchers have long wondered whether ferroelectrics may present topological textures akin to magnetic skyrmions and chiral bubbles, the results being modest thus far. An electric equivalent of a typical magnetic skyrmion would rely on a counterpart of the Dzyaloshinskii-Moriya interaction and seems all but impossible; further, the exotic ferroelectric orders reported to date rely on specific composites and superlattices, limiting their generality and properties. Here, we propose an original approach to write topological textures in simple ferroelectrics in a customary manner. Our second-principles simulations of columnar nanodomains, in prototype material PbTiO3, show we can harness the Bloch-type structure of the domain wall to create objects with the usual skyrmion-defining features as well as unusual ones—including isotopological and topological transitions driven by external fields and temperature—and potentially very small sizes. Our results suggest countless possibilities for creating and manipulating such electric textures, effectively inaugurating the field of topological ferroelectrics.


Magnetic skyrmions or chiral bubbles (MSBs) are spin structures with unusual topological, dynamical, and response properties (15). Skyrmion bubbles are characterized by a nonzero integer topological chargeEmbedded Image(1)where the Pontryagin density q(x, y) is given byEmbedded Image(2)

Here, u = u(x, y) is a vector field that describes the spin order in the xy plane in an idealized continuum limit. The MSB sketched in Fig. 1A has Q = 1; in contrast, the most usual spin arrangements (e.g., ferromagnetic and spin spirals) all present Q = 0. Beyond their fundamental interest, skyrmions hold definite technological promise, e.g., for racetrack memories (6, 7), and constitute a very exciting field in today’s condensed matter physics and materials science. (As mentioned in Supplementary Discussion 1, the distinction between skyrmions and chiral bubbles is a subtle one and usually involves aspects—e.g., long-range chiral order—not considered in this work. Here, we use the terms “skyrmion” and “skyrmion bubble” indistinctly.)

Fig. 1 Concept to create ESBs.

(A) Sketches of a typical Bloch-like MSB. (B) Unit cell of PTO, where arrows mark the displacements yielding a local polarization Pz or, equivalently, a vector field u||(0, 0, 1). Structure of the 180° FE DW of PTO at high (C) and low (D) T, as predicted in (23). (E) ND within a matrix of opposite polarization investigated in this work.

MSBs are typically found in ferromagnets featuring competing interactions whose combined action, often in the presence of thermal activation and external fields, results in nontrivial spin arrangements. Ferroelectrics (FEs) form another well-known family of ferroics where competing couplings abound (8, 9). Hence, one would expect to find in FEs an electric analog of MSBs, with electric dipoles in place of spins. However, electric skyrmion bubbles (ESBs) remain elusive.

The apparent lack of ESBs may be partly due to existing differences between spins and electric dipoles. For example, it is proving all but impossible to find an electric analog of the Dzyaloshinskii-Moriya interaction, which favors noncollinear spin arrangements and is a common ingredient to obtain small MSBs. Electric dipoles are the result of local symmetry-wise polar lattice distortions, whose amplitude can vary continuously (Fig. 1B). Hence, instead of accommodating competing interactions by forming skyrmion bubbles, electric dipoles can just vanish, while spins (typically) cannot.

Nevertheless, recent studies suggest where to find interesting dipole textures, e.g., mediating nucleation (10) or switching (11) of FE domains, or induced by defects (12) or magnetic order (13, 14). Further, vortex-like dipole structures and bubble domains have been observed in heterostructures that combine FE [PbTiO3 (PTO) and PbZr1−xTixO3] and paraelectric [SrTiO3 (STO)] layers (1519), suggesting that exotic orders may occur under appropriate electrostatic boundary conditions. Electrostatics is also responsible for the ESBs predicted to occur in FE (BaTiO3) pillars inside a paraelectric (STO) matrix (2022). While encouraging, these examples rely on complex artificial nanostructures, which limit the generality and properties of the prospective ESBs.

Our approach to ESBs is different. Wojdeł and Íñiguez (23) showed that the common 180° domain walls (DWs) of PTO have a Bloch-like character at low temperatures, with a spontaneous electric polarization confined within the DW plane (Fig. 1, C and D). Simulations further predict that such Bloch DWs occur in the PTO layers of PTO/STO superlattices and cause them to be chiral, which explains recent experimental observations (24). Here, we show how these Bloch DWs also allow us to create ESBs with tunable and unique properties.


Strategy to stabilize ESBs

In view of Fig. 1D, consider the situation in Fig. 1E, where a columnar nanodomain (ND) is embedded in a big matrix of opposite polarization, and the corresponding nanodomain wall (NDW) forms a closed surface. (Let us stress that, here, we assume that these NDs can be written, and do not discuss the conditions for their spontaneous formation.) We run second-principles simulations of NDs like that in Fig. 1E, starting from a configuration where all electric dipoles align strictly along −z (matrix) or +z (ND) and relaxing the structure using the model potential previously applied to PTO (23, 25) and PTO/STO (16, 24, 26). Figure 2A shows the lowest-energy solution thus obtained, which features a Bloch-like NDW. Figure 2D shows the corresponding q(x, y), which presents local maxima at the NDW corners and a total topological charge Q = 1. The simulated ND is thus an ESB. In Supplementary Discussion I and fig. S1, we discuss in some detail how this ESB fits the classification schemes in the MSB literature.

Fig. 2 Prediction of an ESB.

Calculated polarization (A to C) and Pontryagin density (D to F) maps for our ND within a matrix in its ESB ground state (A and D), the same ND-ESB subject to an in-plane electric field along (1,1) (B and E, the field is indicated by a shadowed arrow) and the NDW-polar state stabilized for large enough field values (C and F). In (A) to (C), the color scale gives the out-of-plane Pz component, while the arrows correspond to the in-plane Px and Py. (G) Probability distribution for Q as a function of T. (H) Polarization as a function of in-plane electric field; black-filled squares give |Pz| as obtained at the middle of either matrix or ND (right vertical axis; the results for matrix and ND are essentially identical and very close to those for a monodomain state); blue-filled circles give the Px = Py components (left axis), obtained from a supercell average and normalized to the supercell volume; green open circles give the monodomain result for Px = Py. (I) Energy difference ΔE = ENDW-polarEND-ESB between the NDW-polar and ND-ESB states as obtained in a 16 × 16 × 1 supercell. Note that the ND-ESB and NDW-polar solutions are (meta)stable in the whole field range considered here. In (B) and (E), we show the ND-ESB solution at 500 kV/cm to better visualize the shift of the ESB center; in (C) and (F), we show the NDW-polar solution at zero field.

Figure 2 shows results for a square ND with a small section of about 6 × 6 perovskite cells (~2.32 nm2); yet, as shown in fig. S2, the predicted ESB is robust upon variations of the shape and size of the ND and could be made much bigger. Likewise, the ESB can be as small as the smallest stable ND (see Supplementary Discussion II and fig. S3 for more details). This suggests the possibility of reaching ESB radii of a few nanometers, well beyond what is typical for MSBs. [Very small MSBs have been reported, though (2, 27).]

We check the stability at low temperatures of our predicted ESBs by computing their associated dynamical matrices and verifying that all phonon frequencies are real. Monte Carlo simulations (Fig. 2G) indicate the electric dipoles at the NDW disorder upon heating, yielding a topological transition between skyrmionic (Q = 1) and normal (Q = 0) ND states at TQ ~ 235 K.

Response to electric fields

We now consider the response of our ESB to applied in-plane electric fields (Fig. 2H). We can distinguish two regimes. For moderate fields, we find a seemingly trivial linear dielectric behavior. Yet, the reaction of the ESB (Fig. 2, B and E) is peculiar: For a field applied along (1, 1), the ESB center moves along the perpendicular Embedded Image, which is reminiscent of the response to fields and currents observed in some MSBs (2).

For larger fields, the ESB undergoes a discontinuous topological transition (Fig. 2, H and I) to a state with Q = 0, in which the NDW becomes polarized in-plane (Fig. 2, C and F). Further simulations show that this NDW-polar state is metastable at zero field, with a dipole moment of about 1.4 × 1028 C m and an energy of about 6.7 meV per DW cell above the ESB ground state. Note that our ESB behaves as an antiferroelectric (28).

Control by epitaxial strain

Let us now consider the effect of epitaxial strain, which can be imposed by growing films on suitable substrates and is known to tune the properties of FE perovskites (29). For example, if we start with a PTO film in a z-polarized monodomain state and apply an isotropic tensile strain in the xy plane, then the polarization will eventually rotate by developing Px = Py components. To test whether such a perturbation affects our ESB, we run simulations imposing a square substrate with lattice constants asub = bsub. Our results are summarized in Fig. 3.

Fig. 3 Strain-driven transitions.

(A) Polarization (top) and energy difference between the NDW-polar and ND-ESB states (bottom) as a function of the epitaxial constraint asub; details as in Fig. 2. (B) Position of the ESB center (top) and related susceptibility (bottom) as a function of asub (see the main text for definitions). In the top figure, the dashed line marks the regime where the polar NDW becomes the ground state and the polar ESB is a metastable solution. (C) Polarization map for the polar ESB state that appears for asub ≳ 3.95 Å. (D) For the same state, color map for the orientation of the in-plane polarization (Px, Py), evidencing the formation of 90° DWs.

For asub close to the theoretical bulk value (abulk = 3.93 Å), the ESB remains essentially identical to the solution presented above in bulk-like conditions (i.e., with no elastic constraint). The surprises commence for larger asub values. As shown in Fig. 3A, asub ≈ 3.95 Å marks the onset of the said polarization rotation, in both the monodomain and ND-ESB cases; yet, the Px = Py components are smaller when the ESB is present. Figure 3 (C and D) shows representative results in this regime: Locally, the electric dipoles develop an in-plane component, but they form 90° domains that are compatible with the ESB topology, yielding a Q = 1 multidomain structure with nearly null in-plane polarization. The small Px = Py values obtained in the ESB case for asub ≳ 3.95 Å are mainly due to a break in symmetry that affects the skyrmion bubble itself: The ESB center moves away from the midpoint of the ND, thus developing an in-plane polarization. This polar ESB is very similar to the one obtained above (Fig. 2B) by applying a field to the (high-symmetry, nonpolar) skyrmion bubble. Epitaxial strain allows us to stabilize the polar ESB state in the absence of electric field.

To better characterize the transition at asub ≈ 3.95 Å, we define the position of the ESB center as (30)Embedded Image(3)where α = x, y; rx = x and ry = y. We also introduce the susceptibilityEmbedded Image(4)where Eβ is the β component of an applied electric field. As shown in Fig. 3B, this susceptibility nearly diverges at the skyrmion-skyrmion transition, which reflects its second-order character and the very soft (low-energy) vibrations of the ESB center. Accordingly, as shown in fig. S4, for asub ≳ 3.95 Å, we have a region in which moderate fields can be used to switch the polar ESB among its four symmetry-equivalent states, suggesting a novel possibility for storing information.

Last, for asub ≈ 3.98 Å, we observe a first-order topological transformation to a normal (Q = 0) state that is strongly polar. This solution is all but identical to the NDW-polar state discussed above, obtained under relatively large in-plane electric fields (Fig. 2, C and F). By favoring the occurrence of in-plane dipoles, the tensile strain reverses the relative stability of the ESB and NDW-polar configurations, yielding the latter as the ground state even with no electric field applied.


In sum, our results show that, by writing columnar FE domains within a matrix of opposite polarization, one can stabilize ESBs owing to the Bloch-like structure of the DWs. These skyrmions resemble the soft bubbles of some magnetic materials (1) and have an analogous origin, albeit important differences (electric dipoles can vanish, so it is not obvious that a bubble domain will yield a skyrmion). The predicted skyrmion bubbles show the expected topological properties and some novel ones, including various isotopological and topological transitions. For our skyrmion bubbles to move, one should seek conditions favoring DW mobility, as, e.g., in the “domain liquid” of (16). An experimental characterization of our skyrmion bubbles will be a great challenge; yet, recent studies of related cases using advanced microscopy and diffraction techniques (15, 24, 26) suggest promising strategies to pursue. With plenty of challenges and opportunities ahead, we hope that the present work will propel the field of electric skyrmions.


For the second-principles simulations, we used the same methodology and potentials described in previous works (23, 25), as implemented in the SCALE-UP package (25, 31). The used PTO potential was fitted to density functional theory results within the local density approximation (LDA) (32, 33) and inherits LDA’s well-known overbinding problem; thus, as customarily done in these cases (25, 34), we included an expansive hydrostatic pressure (−13.9 GPa) to compensate for it. The resulting potential yields a qualitatively correct description of the lattice dynamical properties and structural phase transitions of PTO; further, it has been explicitly checked against direct first-principles simulations with regard the Bloch-type structure of PTO’s 180° DWs (23), a result critically important for the present work.

To solve the models, we used standard Monte Carlo and Langevin molecular dynamics methods. Typically, we ran simulated annealing calculations to perform structural relaxations and Metropolis Monte Carlo to compute thermal averages. For the former, we typically worked with periodically repeated supercells made of 16 × 16 × 1 elemental five-atom cells (i.e., 1280 atoms), our ND being formed by approximately 6 × 6 × 1 cells. For the latter, we worked with a 16 × 16 × 10 supercell (i.e., 12,800 atoms) with an embedded 6 × 6 × 10 ND; at a given temperature, we typically ran 20,000 sweeps for thermalization, followed by 20,000 (100,000 around TQ) additional sweeps to compute averages.

Our simulations indicate that the lowest-energy states here discussed have a trivial periodicity of one cell along the z direction. In the case of our simulations at 0 K, the periodicity along z is strict. In the case of our Monte Carlo simulations at finite T, the periodicity pertains to the thermal average, instantaneous departures from it being, of course, possible as the DW dipoles fluctuate.

We computed local polarizations within a linear approximation, using the atomic displacements from the (cubic) reference perovskite structure and the corresponding Born charge tensors. We computed the topological charge Q by (i) processing our local dipoles to obtain a normalized polarization field (| u(x, y) |= 1) and (ii) applying the scheme in (35); we found that this yields well converged results even for our rapidly changing polarization fields. To compute the probability distribution for Q from our Monte Carlo runs, we worked with polarization maps u(x, y) and associated topological densities q(x, y) corresponding to 16 × 16 × 1 slices of our 16 × 16 × 10 supercell.

Note that, typically, a given u(x, y) map can be viewed as a linear superposition of modes with different characteristic topology. For example, around TQ, the system is characterized by fluctuations that resemble the ESB (Q = 1) and polar (Q = 0) solutions discussed in our article, which constitute its lowest-energy excitations. Yet, it is important to realize that Q is not linear in u(x, y); for any given state, even if the structure is the result of a superposition of vibrational modes, our calculation procedure (35) yields an integer Q corresponding to the dominant mode, instead of a linear combination of Q’s. Hence, the corresponding probability distribution pertains only to integer Q values.


Supplementary material for this article is available at

Supplementary Discussion I: Terminology considerations

Supplementary Discussion II: Minimal size of ESBs

Fig. S1. Sketches of simulation supercell and skyrmion structures.

Fig. S2. ESBs of different sizes and shapes.

Fig. S3. Smallest ESB.

Fig. S4. ESB as a four-state memory.

References (3638)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This work was supported by the Luxembourg National Research Fund, through the CORE (C15/MS/10458889 NEWALLS to M.A.P.G. and J.Í.) and AFR (grant 9934186 to C.E.-S.) programs. We were also funded by the Spanish Ministry of Economy and Competitiveness through grants FIS2015-64886-C5-2-P (to P.G.-F. and J.J.) and RyC-2013-12515 (to P.G.-F.). Author contributions: The work was mainly conceived and executed by M.A.P.G. and J.Í. C.E.-S. and P.G.-F. contributed to resolve technical issues, and together with J.J., discuss the results, and prepare the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article