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.
I. INTRODUCTION
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.
II. MODEL AND METHODS
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 nJm/C . The only difference with Refs. 15–17 is the four times smaller magnitude of the gradient coefficients used in the present work, so that pJm /C , pJm /C , and pJm /C . 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 nm 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 carried by the orthogonal set of unit vectors , , and , similarly as in Refs. 14–16.
III. CREATION OF ANTISKYRMION BY INHOMOGENEOUS FIELDS
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 -direction so that in terms of definitions adopted in Ref. 10, it is antiskyrmion.
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.
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.
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 antiskyrmion, with a diameter of about 9 nm. With the adopted model, the smallest stable antiskyrmion at zero field was that with , while in the previous atomistic shell model calculations, conducted at K, the smallest antiskyrmion was that with . Therefore, we have carried out also additional calculations with the Landau potential corresponding to K and verified that there the smallest stable antiskyrmion has , a closer value to the previous atomistic theory prediction.10
IV. ANNIHILATION OF THE ANTISKYRMIONS
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 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 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 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 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.
Annihilation of the 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 , and 9. The vertical line indicates the instant of the topological charge collapse ( ). (b)–(f) Energetically favored transient configurations corresponding to minima 0–3 and the inflection point 4 on panel (a) are projected on the 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.
Annihilation of the 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 , and 9. The vertical line indicates the instant of the topological charge collapse ( ). (b)–(f) Energetically favored transient configurations corresponding to minima 0–3 and the inflection point 4 on panel (a) are projected on the 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.
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.
V. TRANSPORT OF ANTISKYRMIONS BY THE ELECTRIC FIELD
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 ( MV/m) writing field, oriented opposite to the background polarization, and a weaker ( 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 plane.
Suitable electric tweezer field enabling the antiskyrmion transport in the positive sense of -axis (to the right). (a) Planar cross section of the field in the plane and (b) its profile along the dashed line indicated on (a).
Suitable electric tweezer field enabling the antiskyrmion transport in the positive sense of -axis (to the right). (a) Planar cross section of the field in the plane and (b) its profile along the dashed line indicated on (a).
Translation of the antiskyrmion by “tweezer” electric field of Fig. 3 in the direction . 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 component of the polarization in C/m , while the color code in panels (d)–(f) gives the component of the field in MV/m.
Translation of the antiskyrmion by “tweezer” electric field of Fig. 3 in the direction . 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 component of the polarization in C/m , while the color code in panels (d)–(f) gives the component of the field in MV/m.
VI. DISCUSSION AND CONCLUSION
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 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 . 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.
ACKNOWLEDGMENTS
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.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
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).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.