Atomistic model calculations recently revealed that the polarization pattern at the cross section through nanoscale ferroelectric nanodomains in rhombohedral barium titanate shows a clear antiskyrmion attributes such as the net invariant topological charge of minus two. The present work confirms that these topological defects also exist as local minima in the discretized Ginzburg–Landau–Devonshire model, which is widely used in phase-field modeling studies. We explore by phase-field simulations how an electric field can be used to create or move antiskyrmions within a monodomain ferroelectric matrix. The process of field-induced antiskyrmion annihilation reveals the existence of a discrete ladder of intermediate antiskyrmion states with distinct radii but common symmetry, topology, and geometrical properties.

The interest in nanoscale topological defects such as vortices and skyrmions is an important driving factor for the progress of current ferroelectric physics and materials science.1–9 Very recently, a peculiar type of topological defect, the so-called antiskyrmion, was found in the groundstate rhombohedral phase of the classical ferroelectric perovskite, barium titanate.10 More precisely, computer simulations using an ab initio based atomistic model, the shell model of Refs. 11 and 12, have shown that the antiparallel columnar nanodomains in this phase are spontaneously decorated by a particular non-collinear polarization texture, which gives rise to an unusual skyrmion number of 2 per each such nanodomain.10 The relatively thick non-collinear boundary region separating the monodomain ferroelectric background and the inverted nanodomain core could indicate an overall high mobility of this boundary and its shape instability. Instead, it has been found that once the regular antiskyrmion is formed, it is rather stable in its size, shape, and position at least up to temperatures of several tens of Kelvins.10 It has also been inferred that the formation of the antiskyrmion is related to the unusual anisotropy of the local and nonlocal interactions that are typical for certain rhombohedral ferroelectric perovskites, like barium titanate and potassium niobate.10 Considering the surprising geometry and topology of these antiskyrmionic polarization patterns, it is highly desirable to seek their confirmation within a radically different theoretical framework.

In the present work, we thus explored the properties of antiskyrmionic nanodomains using a phenomenological Ginzburg–Landau–Devonshire model based on the systematic symmetry-based free-energy expansion and a simulated annealing method. In this way, an independent verification of the previous predictions based on the molecular dynamics with the atomistic model was obtained. Moreover, the fact that the standard version on the discretized Ginzburg–Landau–Devonshire model is general enough to include the antiskyrmionic solutions also sheds a new light on the proposed stabilization mechanism.10 Last but not least, the availability of antiskyrmions within the phase-field approach gives an access to timescales and lengthscales larger than those accessible in the calculations with the atomistic models. In this work, we have exploited this advantage of the phase-field approach and explored in more detail the electric field driven process of formation and annihilation of antiskyrmions as well as the targeted transport of antiskyrmionic ferroelectric bubble domains.

The phase-field approach consists of a numerical simulation of the energy-optimizing relaxation process for the spatially varying field of the electric polarization, subject to the dissipative Landau–Khalatnikov dynamics and the direct interaction with the homogeneous and inhomogeneous strain fields, described by various interaction terms in a given Ginzburg–Landau–Devonshire functional.

The analytical definition of the model employed in the present work was fully described in Refs. 13 and 14. For the sake of comparison with our previous work, we have adopted Landau–Devonshire model parameters from Refs. 15 and 16, and so it is the model B of Ref. 17 with the temperature fixed to 118 K and q 1212 = 1.57 nJm/C 2. The only difference with Refs. 15–17 is the four times smaller magnitude of the gradient coefficients used in the present work, so that G 1111 = 127.5 pJm 3/C 2, G 1122 = 5 pJm 3/C 2, and G 1212 = 5 pJm 3/C 2. The reduction of the gradient terms was primarily motivated by the enhancement of the lattice pinning and the earlier comparisons of ab initio and phase-field domain wall profiles,18 which had already indicated that the gradient coefficients estimated from the experiments in the cubic phase of barium titanate are not appropriate for a quantitative description of thicknesses of domain walls in its low-temperature rhombohedral phase.

The simulations were carried out with a discrete time step of 1.25 ps under a constant volume and shape of the supercell, the latter corresponding to the mechanically free monodomain state with the spontaneous polarization along the pseudocubic [111] direction. The [111]-polarized monodomain state was used as the initial configuration. We have typically worked with a rectangular supercell with dimensions 32 × 32 × 32 nm 3 along the crystallographic axes of the parent cubic perovskite phase. The polarization field was discretized on a regular mesh with a spatial step of 0.5 nm and subject to periodic boundary conditions. The results are presented in a symmetry-adapted orthonormal coordinate system { r s t } carried by the orthogonal set of unit vectors r [ 111 ], s [ 1 ¯ 10 ], and t [ 1 ¯ 1 ¯ 2 ] ], similarly as in Refs. 14–16.

In contrast to the former work,10 where antiskyrmions were obtained by relaxing from an existing initial inverted antiparallel ferroelectric nanodomain of a hexagonal or triangular cross section, in this work, we have created the antiskyrmions by various local electric fields, comparable to the homogeneous thermodynamical coercive field value for 180 ° switching, which is about 60 MV/m for the adopted Landau–Devonshire model.

For example, we applied a writing field, of 80 MV/m, oriented opposite to the background monodomain polarization, within an infinitely long cylinder of a 6 nm diameter, as schematically shown in Fig. 1(a). This writing field reversed the polarization inside of the cylinder by 180 ° within a fraction of ps. In the next few ps, the non-collinear texture was spontaneously formed around the curved boundary of this nanodomain. This texture had the same qualitative geometrical and topological characteristics as the antiskyrmions reported in the shell-model molecular dynamics study of Ref. 10. The nanodomain and its antiskyrmionic character persist also when the field is removed. The final pattern with a core diameter of about 7 nm is shown in Fig. 1(b). This diameter corresponds to 17 spatial steps of the discrete grid in the s-direction so that in terms of definitions adopted in Ref. 10, it is N = 17 antiskyrmion.

FIG. 1.

Creation of the antiskyrmion by the electric field. (a) Schema showing the pseudocubic crystallographic axes and orientation of the locally homogeneous field used to create antiskyrmion. (b) Cross section through the stable antiskyrmion nanodomain formed after the application and removal of the 80 MV/m field, homogeneous within an infinite tube of 6 nm diameter. The arrows stand for the in-plane polarization, the segment in the corner indicates the in-plane magnitude scale, and colors show the out-of-plane polarization component.

FIG. 1.

Creation of the antiskyrmion by the electric field. (a) Schema showing the pseudocubic crystallographic axes and orientation of the locally homogeneous field used to create antiskyrmion. (b) Cross section through the stable antiskyrmion nanodomain formed after the application and removal of the 80 MV/m field, homogeneous within an infinite tube of 6 nm diameter. The arrows stand for the in-plane polarization, the segment in the corner indicates the in-plane magnitude scale, and colors show the out-of-plane polarization component.

Close modal

The same regular antiskyrmions could be created by the axially symmetric field with a Gaussian spatial profile. For example, the field with a full width at half maximum of 4 nm and the peak value of 175 MV/m applied during 0.1 ps also created the 7 nm-diameter antiskyrmion shown in Fig. 1(b). When the field was applied longer, the final diameter increased. For example, 1 ps exposure to the same Gaussian field resulted in the N = 23 antiskyrmion, with a diameter of about 9 nm. With the adopted model, the smallest stable antiskyrmion at zero field was that with N = 13, while in the previous atomistic shell model calculations, conducted at T = 1 K, the smallest antiskyrmion was that with N = 9. Therefore, we have carried out also additional calculations with the Landau potential corresponding to T = 1 K and verified that there the smallest stable antiskyrmion has N = 11, a closer value to the previous atomistic theory prediction.10 

The opposite process, that of the antiskyrmion contraction and its irreversible collapse, can be performed by an electric field applied parallel to the background polarization. The orientation of this field penalizes the core polarization and favors that of the background. Therefore, this field promotes a contraction of the antiskyrmion radius while keeping its macroscopic 3 m symmetry. Since the skyrmionic states are locally stable, the field needed to start this process needs to overcome a threshold value. This value is less than one order of magnitude lower than the thermodynamic coercive field.10 Therefore, the erasing field can be homogeneous throughout the entire simulation box. We have numerically determined that the threshold homogeneous field needed to initiate the process of the irreversible contraction of the N = 17 antiskyrmion in the present model is about 2.4 MV/m.

The antiskyrmion annihilation driven by a 10 MV/m homogeneous field is documented in Fig. 2. Although the simulated temporal evolution is a continuous process, it proceeds by a succession of faster and slower episodes, the latter implying passages through several distinct quasistationary states. These quasistationary states, enumerated in Fig. 2 by sequential numbers 0–4, correspond to symmetric antiskyrmionic configurations with odd size numbers N = 17 , 15 , 13 , 11 , and 9, respectively. Let us note that the [111]-projection of the three-dimensional discrete grid forms a two-dimensional lattice of points in which the local polarization is evaluated, and that the nearest neighbors form elementary equilateral triangular plaquettes (see Fig. 2). We can see that the odd values of N permit simultaneous location of the axis of the antiskyrmion on the lattice point and the six circumferential vortexes at the centers of the elementary triangular plaquettes. This clearly indicates that pinning of vortex cores to the discrete lattice is essential for the stability of antiskyrmions.

FIG. 2.

Annihilation of the N = 17 antiskyrmion by a 10 MV/m homogeneous field, parallel to the background polarization. (a) Time evolution of the overall average elastic Gibbs free energy density during the field-induced process of antiskyrmion collapse. Numbers 0–4 within the panel indicate the temporal location the quasi-stationary configurations with N = 17 , 15 , 13 , 11, and 9. The vertical line indicates the instant of the topological charge collapse ( 2 0). (b)–(f) Energetically favored transient configurations corresponding to minima 0–3 and the inflection point 4 on panel (a) are projected on the s t plane. Arrows show the in-plane components of polarization, while colors stand for the invariant skyrmion density calculated at each elementary triangular plaquette within the [111]-projection of the three-dimensional discrete simulation grid.

FIG. 2.

Annihilation of the N = 17 antiskyrmion by a 10 MV/m homogeneous field, parallel to the background polarization. (a) Time evolution of the overall average elastic Gibbs free energy density during the field-induced process of antiskyrmion collapse. Numbers 0–4 within the panel indicate the temporal location the quasi-stationary configurations with N = 17 , 15 , 13 , 11, and 9. The vertical line indicates the instant of the topological charge collapse ( 2 0). (b)–(f) Energetically favored transient configurations corresponding to minima 0–3 and the inflection point 4 on panel (a) are projected on the s t plane. Arrows show the in-plane components of polarization, while colors stand for the invariant skyrmion density calculated at each elementary triangular plaquette within the [111]-projection of the three-dimensional discrete simulation grid.

Close modal

In addition to the in-plane polarization components, the color circles depict the invariant skyrmion density.10 We can see that the vortex centers in the stationary states coincide with the sharp hot spots of the skyrmion density. As expected, in all (b)–(f) panels, the skyrmion density integrates the overall topological charge of minus two. During the passage between subsequent stationary states, the overall topological charge was preserved. The integral topological charge vanished abruptly only when the diameter was so small that the polarization at the center of the core itself vanished. This instant is indicated in Fig. 2(a) by a vertical dashed line.

Finally, we have explored the possibility of moving the antiskyrmion in a desired direction by suitably tailored electric fields. We have observed that a simple spatial translation of the writing field described above can lead to the formation of various elongated nanodomains with irregular or significantly deformed shape. An efficient method allowing the transport of the antiskyrmion without much shape distortion consisted of the usage of the bipolar electric field depicted in Fig. 3. This field was a superposition of two electric field distributions, a stronger ( E 52 MV/m) writing field, oriented opposite to the background polarization, and a weaker ( E + 18 MV/m) erasing field, pointing in the same direction as the background polarization. The writing and the erasing fields were both homogeneous within the circular diameter adapted to the antiskyrmion of interest (in the present case, 6 nm), but they were mutually shifted by about 1.5–2 nm. The resulting electric “tweezer” field enabled to shift the antiskyrmion core center toward the center of the writing field. When such a tweezer field is translated in space with a velocity slower than about 500 m/s, the antiskyrmion can travel together with the field over an arbitrary distance (see Fig. 4). By changing the direction of the mutual shift between the two fields, the antiskyrmion could be equally well dragged along any other direction in the s t plane.

FIG. 3.

Suitable electric tweezer field enabling the antiskyrmion transport in the positive sense of s-axis (to the right). (a) Planar cross section of the field in the s t plane and (b) its profile along the dashed line indicated on (a).

FIG. 3.

Suitable electric tweezer field enabling the antiskyrmion transport in the positive sense of s-axis (to the right). (a) Planar cross section of the field in the s t plane and (b) its profile along the dashed line indicated on (a).

Close modal
FIG. 4.

Translation of the antiskyrmion by “tweezer” electric field of Fig. 3 in the direction s. Panels (a)–(c) show polarization corresponding to the drifting electric tweezer field shown (d)–(f) at three subsequent instants of the simulation. The color code in panels (a)–(c) indicates P r component of the polarization in C/m 2, while the color code in panels (d)–(f) gives the E r component of the field in MV/m.

FIG. 4.

Translation of the antiskyrmion by “tweezer” electric field of Fig. 3 in the direction s. Panels (a)–(c) show polarization corresponding to the drifting electric tweezer field shown (d)–(f) at three subsequent instants of the simulation. The color code in panels (a)–(c) indicates P r component of the polarization in C/m 2, while the color code in panels (d)–(f) gives the E r component of the field in MV/m.

Close modal

In summary, the present phase-field calculations confirmed the existence of the string-like antiskyrmionic bubble domain solutions in the discretized Ginzburg–Landau–Devonshire model. The adopted Landau–Devonshire potential had standard parameters for rhombohedral barium titanate. The resulting antiskyrmion patterns were qualitatively similar to those encountered in the original atomistic simulations.10 In particular, the polarization pattern is composed from sextuplet of polar vortexes, it has an overall 3 m symmetry, and the overall invariant skyrmion topological charge is minus two. From this perspective, these low-temperature ferroelectric bubble domains are very different from lead-titanate ones.19,20

Present simulations focused on several field-induced processes. The spontaneous phase of the antiskyrmion formation process was activated by the application of a focused electric field impulse. Typically, the above-coercive field creates a narrow beam of the inverted polarization, which is subsequently decorated by the characteristic topologically nontrivial, noncollinear polarization pattern. The transport of antiskyrmionic nanodomains was achieved by a suitably designed bipolar tweezer field that combines off-centered writing and erasing fields of opposite sense. As pointed out in one of our review reports, a similar effect might be achieved by considering fields generated by scanning over the crystal surface by a mechanically oscillating sharp tip, carrying a synchronously alternating electric voltage.

The analysis of the process of irreversible annihilation of the antiskyrmion by a homogeneous sub-coercive field revealed a sequence of stationary states possessing odd values of the dimensionless diameter N. In overall, the reported annihilation process provides an instructive example of a topological collapse in an already energetically destabilized transient state.

The above findings confirm that the discrete lattice pinning of vortex centers plays an essential role in the stability of these antiskyrmionic nanodomains. Therefore, it can be anticipated that the presence of thermal fluctuations will further facilitate the transport and annihilation of the antiskyrmions at finite temperatures. At the same time, some of the numerical values obtained in this study, such as the critical anstiskyrmion diameters or the threshold field magnitudes are expected to be considerably sensible to the model parameters and in particular, to the magnitude of the gradient tensor components and to the size of the discrete spatial step. Similarly, the temporal rates should be considered with a caution since the value of the kinetic parameter was derived from a very simple assumption and from the room-temperature relaxation time.21 Nevertheless, we did not aim here to approach real barium titanate properties by a tight fine-tuning of all these parameters. We believe that in the near future, a more precise quantitative prediction for the nanoscale phenomena in the low-temperature of barium titanate will become available in a more straightforward and convincing manner from the recently developed ab initio based effective Hamiltonians,22–24 from the systematic second-principles models25,26 or also from other carefully designed and tested interatomic potentials.27 

The formal mathematical link between topological solitons in the theory of mesons by Skyrme28,29 and condensed matter topological defects of various dimensionality and nature involve abstract geometrical constructions such as the conformal mapping of the unit sphere onto a certain plane in the material. In the context of the vast literature on skyrmions and bubble domains in ferroic materials, it is worth to emphasize two aspects of the antiskyrmions addressed in this work, (i) they have the form of axisymmetric nanoscale right cylinders and (ii) they exist due to the intrinsic bulk interactions in the crystal. In that sense, the antiskyrmions reported here are closely analogous to the original magnetic vortices predicted in late 90s30 and, thus, to the first experimentally discovered magnetic skyrmions, evidenced in chiral magnets by neutron scattering studies.31 It would be interesting to explore also the finite-length antiskyrmions, terminated at the surface of the ferroelectric thin film. Resulting surface-termination effects can certainly bring additional interesting properties and possibilities that in the case of ferroelectric antiskyrmions were not explored yet.

In conclusion, the present phase-field calculations should give a realistic qualitative picture of antiskyrmion formation, sliding, and annihilation processes. Most likely, the present results also indicate experimentally relevant orders of magnitude of timescales and lengthscales of these processes. We believe that these findings will stimulate further research into these interesting topological ferroelectric bubble domains.

This work was supported by the Czech Science Foundation (Project No. 19-28594X). The authors acknowledge multiple inspiring scientific exchanges with Mauro A. P. Gonçalves, Pavel Márton, and Marek Paściak from FZU and with Florian Mayer from MCL Leoben.

The authors have no conflicts to disclose.

V. Stepkova: Conceptualization (lead); Data curation (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (lead); Methodology (equal); Project administration (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (supporting); Writing – review & editing (equal). J. Hlinka: Conceptualization (lead); Data curation (supporting); Formal analysis (equal); Funding acquisition (lead); Investigation (supporting); Methodology (equal); Project administration (lead); Software (lead); Validation (supporting); Visualization (supporting); Writing – original draft (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
J.
Junquera
et al., “
Topological phases in polar oxide nanostructures
,”
Rev. Mod. Phys.
95
,
025001
(
2023
).
2.
S.
Das
et al., “
Observation of room-temperature polar skyrmions
,”
Nature (London)
568
,
368
(
2019
).
3.
Y. J.
Wang
et al., “
Entangled polarizations in ferroelectrics: A focused review of polar topologies
,”
Acta Mater.
243
,
118485
(
2023
).
4.
G. F.
Nataf
et al., “
Domain-wall engineering and topological defects in ferroelectric and ferroelastic materials
,”
Nat. Rev. Phys.
2
,
634
(
2020
).
5.
X.
Guo
et al., “
Theoretical understanding of polar topological phase transitions in functional oxide heterostructures: A review
,”
Small Methods
6
,
2200486
(
2022
).
6.
C.
Halcrow
and
E.
Babaev
, “
Ferroelectric domain wall clusters in perfectly screened barium titanate type systems
,”
Phys. Rev. B
108
,
174101
(
2023
).
7.
C.
Halcrow
and
E.
Babaev
, “
Fractional skyrme lines in ferroelectric barium titanate
,”
Phys. Rev. Res.
6
,
L032011
(
2024
).
8.
R.
Zhu
et al., “
Dynamics of polar skyrmion bubbles under electric fields
,”
Phys. Rev. Lett.
129
,
107601
(
2022
).
9.
S.
Prokhorenko
et al., “
Motion and teleportation of polar bubbles in low-dimensional ferroelectrics
,”
Nat. Commun.
15
,
412
(
2024
).
10.
M. A. P.
Gonçalves
,
M.
Paściak
, and
J.
Hlinka
, “
Antiskyrmions in ferroelectric barium titanate
,”
Phys. Rev. Lett.
133
,
066802
(
2024
).
11.
S.
Tinte
et al., “
Ferroelectric properties of Ba xSr 1 xTiO 3 solid solutions obtained by molecular dynamics simulation
,”
J. Phys.: Condens. Matter
16
,
3495
(
2004
).
12.
M.
Sepliarsky
et al., “
Atomic-level simulation of ferroelectricity in oxide materials
,”
Curr. Opin. Solid State Mater. Sci.
9
,
107
(
2005
).
13.
J.
Hlinka
and
P.
Márton
, “
Phenomenological model of a 90 ° domain wall in BaTiO 3-type ferroelectrics
,”
Phys. Rev. B
74
,
104104
(
2006
).
14.
P.
Marton
,
I.
Rychetsky
, and
J.
Hlinka
, “
Domain walls of ferroelectric BaTiO 3 within the Ginzburg-Landau-Devonshire phenomenological model
,”
Phys. Rev. B
81
,
144125
(
2010
).
15.
V.
Stepkova
,
P.
Marton
, and
J.
Hlinka
, “
Stress-induced phase transition in ferroelectric domain walls of BaTiO 3
,”
J. Phys. Condens. Mat.
24
,
212201
(
2012
).
16.
V.
Stepkova
,
P.
Marton
, and
J.
Hlinka
, “
Ising lines: Natural topological defects within ferroelectric Bloch walls
,”
Phys. Rev. B
92
,
094106
(
2015
).
17.
J.
Hlinka
et al., “Ferroelectric domain walls and their intersections in phase-field simulations,” in Topological Structures in Ferroic Materials: Domain Walls, Vortices and Skyrmions, edited by J. Seidel (Springer International Publishing, Cham, 2016), Vol. 228, p. 161.
18.
M.
Taherinejad
et al., “
Bloch-type domain walls in rhombohedral BaTiO 3
,”
Phys. Rev. B
86
,
155138
(
2012
).
19.
F.
Gómez-Ortiz
et al., “
Liquid-crystal-like dynamic transition in ferroelectric-dielectric superlattices
,”
Phys. Rev. Lett.
133
,
066801
(
2024
).
20.
H.
Aramberri
and
J.
Íñiguez-González
, “
Brownian electric bubble quasiparticles
,”
Phys. Rev. Lett.
132
,
136801
(
2024
).
21.
J.
Hlinka
, “
Mobility of ferroelastic domain walls in barium titanate
,”
Ferroelectrics
349
,
49
(
2007
).
22.
F.
Mayer
et al., “
Improved description of the potential energy surface in BaTiO 3 by anharmonic phonon coupling
,”
Phys. Rev. B
106
,
064108
(
2022
).
23.
T.
Nishimatsu
et al., “
Fast molecular-dynamics simulation for ferroelectric thin-film capacitors using a first-principles effective Hamiltonian
,”
Phys. Rev. B
78
,
104104
(
2008
).
24.
M.
Marathe
et al., “
Electrocaloric effect in BaTiO 3 at all three ferroelectric transitions: Anisotropy and inverse caloric effects
,”
Phys. Rev. B
96
,
014102
(
2017
).
25.
J.
Zhang
et al., “
Structural phase transitions and dielectric properties of BaTiO 3 from a second-principles method
,”
Phys. Rev. B
108
,
134117
(
2023
).
26.
P.
García-Fernández
et al., “
Second-principles method for materials simulations including electron and lattice degrees of freedom
,”
Phys. Rev. B
93
,
195137
(
2016
).
27.
Y.
Qi
et al., “
Atomistic description for temperature-driven phase transitions in BaTiO 3
,”
Phys. Rev. B
94
,
134308
(
2016
).
28.
T. H. R.
Skyrme
, “
A non-linear field theory
,”
Proc. R. Soc. Lond. A
260
,
127
(
1961
).
29.
T. H. R.
Skyrme
, “
Particle states of a quantized meson field
,”
Proc. R. Soc. Lond. A
262
,
237
(
1961
).
30.
A. N.
Bogdanov
and
D. A.
Yabloskii
, “
Thermodynamically stable ‘vortices’ in magnetically ordered crystals, the mixed state of magnets
,”
Zh. Eksp. Teor. Fiz.
95
,
178
(
1989
).
31.
S.
Mühlbauer
et al., “
Skyrmion lattice in a chiral magnet
,”
Science
323
,
915
(
2009
).