# Controllable skyrmion chirality in ferroelectrics

Yu Tikhonov, S. Kondovych, J. Mangeri, M. Pavlenko, L. Baudry, A. Sene, A.

Galda, S. Nakhmanson, O. Heinonen, A. Razumnaya, et al.

► To cite this version:

Yu Tikhonov, S. Kondovych, J. Mangeri, M. Pavlenko, L. Baudry, et al. Controllable skyrmion chirality in ferroelectrics. Scientific Reports, 2020, 10 (1), pp.8657. 10.1038/s41598-020-65291-8. hal-03321562

HAL Id: hal-03321562

https://hal.science/hal-03321562

Submitted on 20 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Controllable skyrmion chirality in ferroelectrics

Yu. Tikhonov<sup>1,2</sup>, S. Kondovych<sup>2,3</sup>, J. Mangeri<sup>4,5</sup>, M. Pavlenko<sup>1</sup>, L. Baudry<sup>6</sup>, A. Sené<sup>2</sup>, A. Galda<sup>7</sup>, S. Nakhmanson<sup>5,8</sup>, O. Heinonen<sup>9</sup>, A. Razumnaya<sup>1</sup>, I. Luk'yanchuk<sup>2,10</sup> & V. M. Vinokur<sup>9,11</sup> ❙

Chirality, an intrinsic handedness, is one of the most intriguing fundamental phenomena in nature. Materials composed of chiral molecules find broad applications in areas ranging from nonlinear optics and spintronics to biology and pharmaceuticals. However, chirality is usually an invariable inherent property of a given material that cannot be easily changed at will. Here, we demonstrate that ferroelectric nanodots support skyrmions the chirality of which can be controlled and switched. We devise protocols for realizing control and efficient manipulations of the different types of skyrmions. Our findings open the route for controlled chirality with potential applications in ferroelectric-based information technologies.

Chirality, a fundamental asymmetry property describing systems that are distinguishable from their mirror images, remains in the focus of modern science<sup>1-4</sup>, and chiral materials find diverse applications<sup>5-8</sup>. Chiral topological textures set the stage for a new generation of chiral materials, where the chirality is extended over nano- and micro-scales. Nonuniform chiral states, helical, blue, and twist grain boundary (TGB) phases have been observed in cholesteric liquid crystals<sup>9,10</sup>. Skyrmions, which are the chiral texture of a vector order parameter, such as magnetization or polarization density, have been attracting considerable attention in magnetic materials during the past decade<sup>11-13</sup> due to their potential applications in information technologies. However, a salient feature of these materials is the specific non-chiral symmetry, carried either by non-mirror-symmetric molecules in cholesterics or the antisymmetric spin exchange in magnetic systems leading to the Dzyaloshinskii-Moriya spin interaction. Recently, an extension of the class of magnets hosting skyrmions onto systems without Dzyaloshinskii-Moriya spin interaction has been reported<sup>14,15</sup>. However, the possibility of tuning the chirality of skyrmions in these systems remains an open question.

Although a pre-defined chiral symmetry is absent in ferroelectric materials, they were recently found to host a wealth of chiral topological excitations, including Bloch domain walls<sup>16-19</sup>, coreless vortices with a skyrmion structure<sup>20-22</sup>, single skyrmions<sup>23,24</sup>, skyrmion lattices<sup>25</sup>, and Hopfions<sup>26</sup>. A distinct feature of ferroelectrics is that the chirality appears as a result of the spontaneous symmetry breaking due to specific interplay of confinement and depolarization effects when the depolarization charges $ρ = \nabla \cdot P$ rearrange to reduce their interaction energy, leading to the chiral twisting of the polarization. Importantly, the different chiralities ("left" and "right" states) are energetically degenerate and hence inter-switchable. However, performing such chirality-switching poses a challenge because of the non-chiral nature of the fundamental fields that could serve as a control parameter. We find that ferroelectric nanodots can provide rich phase diagrams as depolarization effects lead to an abundance of topological excitations, and we demonstrate that ferroelectric nanodots harbor polarization skyrmions. In particular, we devise a system in which controlled switching between the opposite chiralities may be implemented by the applied electric field.

<footnote n="1"><sup>1</sup>Faculty of Physics, Southern Federal University, 5 Zorge str., 344090, Rostov-on-Don, Russia. <sup>2</sup>University of Picardie, Laboratory of Condensed Matter Physics, Amiens, 80039, France. <sup>3</sup>Life Chemicals Inc., Murmanska st. 5, Kyiv, 02660, Ukraine. <sup>4</sup>Institute of Physics, Academy of Sciences of the Czech Republic, Na Slovance 2, 18221, Praha 8, Czech Republic. <sup>5</sup>Department of Physics, University of Connecticut, Storrs, CT, USA. <sup>6</sup>Institute of Electronics, Microelectronics and Nanotechnology (IEMN)-DHS Département, UMR CNRS 8520, Université des Sciences et Technologies de Lille, 59652, Villeneuve d'Ascq Cedex, France. <sup>7</sup>James Franck Institute, University of Chicago, Chicago, Illinois, 60637, USA. <sup>8</sup>Department of Materials Science & Engineering and Institute of Material Science, University of Connecticut, Storrs, Connecticut, 06269, USA. <sup>9</sup>Materials Science Division, Argonne National Laboratory, 9700S. Cass Avenue, Argonne, Illinois, 60637, USA. <sup>10</sup>L. D. Landau Institute for Theoretical Physics, Akademika Semenova av., 1A9, Chernogolovka, 142432, Russia. <sup>11</sup>Consortium for Advanced Science and Engineering (CASE) University of Chicago, 5801S Ellis Ave, Chicago, IL, 60637, USA. ❙e-mail: vinokour@anl.gov</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Ferroelectric switch circuit and skyrmion states. (a) The circuit is controlled by the external switching voltage *U*. The top electrode is separated from the ferroelectric nanodot carrying the polarization topological states. (b) Four types of the skyrmions, differing by their chirality and polarity. The hand pictograms define the classification of the skyrmions.</figcaption></figure>

**The model.** Our target system is a ferroelectric nanodot in a shape of the disk deposited on a substrate (see Fig. 1a). We choose the lead titanate pseudo-cubic perovskite oxide, PbTiO$_3$, as the model ferroelectric material. The typical nanodots that we use for observation of the skyrmion have diameter about 40 nm and thickness about 20 nm. For calculations we employ the phase field approach, described in details in the Methods section. The ferroelectric nanodot is located in a capacitor, the upper plate of which is separated from the sample either by vacuum or by a dielectric material with low dielectric constant. The thin lower electrode forms an interface between the nanodot and the substrate, which induces a weak strain caused by the ferroelectric lattice mismatch with the substrate. The strain results in an out-of-plain polarization anisotropy (see Methods for relevant parameters). This system enables the creation and manipulation of topological nonuniform textures confined in the nanodot by applying voltage *U* to the electrodes.

The emergence of these topological excitations is related to the effect of depolarization charges that arise because of the abrupt termination of the spontaneous polarization at the nanodot surface. Induced depolarization fields destabilize the uniformly-polarized state, resulting in a polarization texture corresponding to a self-consistent local energy minimum. In thin-film geometry, the depolarization effect yields Kittel domains$^{27}$, structured as soft polarization stripes$^{28,29}$ (or parallel-anti-parallel vortices$^{30,31}$), or bubble domains in a skyrmion lattice structure$^{25}$. Confinement of these textures to within a nanodot gives rise to the competition between striped and cylindrical domains$^{32}$. We demonstrate that the latter configuration generates chiral skyrmions, provided that the substrate-induced uniaxial anisotropy is weak and that the polarization draws up a Bloch twisting of the polarization inside the domain wall. The left and right chiral skyrmions, shown in Fig. 1b, are monotonically visualized by left and right hands, where thumbs point in the direction of the polarization orientation in the core, which is referred to as skyrmion polarity, and the fingers curl along the orientations of the winding of the polarization around the skyrmion core corresponding to the skyrmion chirality. Accordingly, we introduce the notations $L_+$, $L_-$, $R_+$, and $R_-$, where the “L” and “R” correspond to the left and right-hand skyrmions and “+” and “-” subscript denote the “up” and “down” polarities, respectively.

**The structure of a skyrmion.** The polarization distribution inside the skyrmion confined within the nanodot is shown in Fig. 2a. The polarization texture preserves the structure of the polarization rotation inside the chiral Bloch-like circular domain wall over almost entire height of the sample. It is within these chiral domain walls where the chirality of the nanodot is concentrated. At the top near-surface layer the polarization configuration assumes the sinc-like shape to form a Néel-type non-chiral skyrmion to maintain the polarization tangential to the surface, see the top view in Fig. 2b, in order to prevent the formation of surface depolarization charges. The entire texture resembles the structure of the bubble domain in a double-periodic domain structure observed in ferroelectric superlattices$^{25}$, and has a topology of the Hopf fibration similar to that in the ferroelectric nanoparticles$^{26}$. In what follows we will construct protocols that enable formation of a skyrmion as well as electric field-tuned transitions between different polarization configurations. The switching between the different states will be described as switching between different mean chiralities, which we define as

$$ \bar{x} = \frac{1}{V} \int_V \chi(r) dr, \quad \chi = \mathbf{P} \cdot [\nabla \times \mathbf{P}], \qquad (1) $$

where $\chi(r)$ is a chirality density, and integration is performed over the nanodot volume $V$, so that the $\bar{x} > 0$ for right- and $\bar{x} < 0$ for left skyrmions. It is important to clearly distinguish between chirality defined by Eq. (1) and
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Field-tuned topological states in the nanodot. (a) Cross section of the nanodot displaying the polarization distribution, with white arrows showing the direction of the polarization. The polarization rotation over 180° in the plane of the Bloch domain walls results in the chirality distribution, $\chi(r)$, shown by the colour map. The legend to the map is given below in low-right corner of the figure. The crossed circle, $\otimes$, at the domain wall denotes the polarization vector going into the cross-section plane, and the circle with the central dot, $\odot$, stands for the out-of-plane polarization. (b) The top view of the sinc-like distribution of the polarization at the near-surface layer of the nanodot. (c) Hysteresis behaviour of the polarization of the nanodot as a function of the applied field. The blue and red branches correspond to the up-down and down-up sweeps of the applied field. The numbers mark the different topological states of the polarization. (d) Hysteresis protocols of the chirality switching that allow to come to the $L_\pm$ and $R_\pm$ skyrmion states. The arrows show the direction of the sweep. The gray branch corresponds to the virgin curve of the poling of the nanodot reflecting both, positive and negative directions of the electric field variation. The blue and red branches again correspond to the up-down and down-up sweeps of the applied field. The violet and yellow branches correspond to the reversal of the field sweep from the blue and red branches, respectively. (e) The distribution of the polarization and chirality in the original polarization state (0) and in the sequence of the topological states arising during the re-polarization of the nanodot by the applied field (view from the bottom) that follows the blue branch of panels c and d. The yellow points mark the cross-section of the Bloch lines piercing the nanodot.</figcaption></figure>

identifying the objects that cannot be mapped to their mirror images by rotations and translations¹, and another swirling characteristics of the vector fields, vorticity, quantified as $\nabla \times P$ and toroidal moment $\int r \times (P - \bar{P}) dV$. The controlled manipulation by four distinct skyrmion states, $R_\pm$ and $L_\pm$, offers the opportunity for implementing a platform for ferroelectric-based multivalued logic units³²,³³.

**The behaviour of a skyrmion in the electric field.** In order to elaborate on the protocols for controlling and tailoring skyrmions, we investigate the response of the polarization in the nanodot to applied electric field $U$. The mean polarization ($\bar{P}(U)$) and the chirality ($\bar{\chi}(U)$) hysteresis curves are shown in Fig. 2c,d respectively. The corresponding stages of the process are displayed in Fig. 2e. We first set the zero-field ferroelectric state at room temperature by quenching it from the paraelectric state with randomly oriented small-amplitude polarization.
The resulting state, which we denote as 0-state, has in Fig. 2e a structure of four-band stable domain stripes³². The domains are separated by the Bloch domain walls in which the direction of the polarization rotations determines the domain wall chirality that can be either positive (shown in red) or negative (shown in blue). At the loci where the chiralities of the opposite sign meet, a linear topological defect, so-called Bloch line (denoted as yellow dot in the cross-section images) penetrating the nanodot forms. The Bloch lines were observed in ferromagnetic domain walls³⁴ and predicted to appear in ferroelectrics³⁵. It is the dynamics of the Bloch lines that eventually controls the chirality switch in the nanodot. The precise distribution of the chiralities in the Bloch domain walls in the initial state of course depends on which random paraelectric configuration was quenched, but does not affect the polarization evolution after the initial poling in a large electric field.

To create and manipulate the skyrmion, we ramp up the electric field by applying a positive voltage $U$ to the circuit, necessary to uniformly polarize the sample. The virgin curve (gray color line) passes through the polarization stripe states differing by the structure of domain walls which move and interswitch their chiralities, see Supplementary Video. The similar effect is achieved by varying the field in the negative direction. The evolution concludes with the jump into the skyrmion state at the latest stage before complete poling is achieved (state 1 in Fig. 2d.). We then start reducing the field strength to zero and, subsequently, having reversed the field direction, increase its magnitude. By analyzing the evolution of the polarization, we establish a protocol that allows us to switch the chirality of the system, as shown in Fig. 2d,e. We consider first the evolution of the system under decreasing field from state 1, depicted by the blue line in Fig. 2d. The details of the transformation of the polarization distribution are visualized in the video in the SI. Upon decreasing the electric field below the threshold value, the polarization in the central region of the nanodot switches its orientation to the opposite one. As a result, the bubble domain (state 2) forms, for which the cylindrical domain wall partitions into two half-cylindrical segments with opposite chiralities separated by Bloch lines, the total chirality remaining zero. As the field strength is further reduced, the skyrmion $L_-$ (state 3) with negative chirality $\bar{\chi}$ forms as a result of merging the Bloch lines and a concomitant collapse of the positive-chirality segment, and a skyrmion appears (state 3). Further reducing the field strength, the thickness of the skyrmion core grows, leading to repartitioning of the up- and down-oriented polarization regions (state 4) and, hence in the $\bar{P}(U)$-dependence shown in Fig. 2c.

The skyrmion texture remains as the field is completely removed (state 4), and even as the field direction is reversed, although the texture now becomes metastable. At some negative field (state 5) the skyrmion becomes unstable and decays into a multi-domain state (state 6) with zero chirality. The arc-shape domain walls form three pairs of segments with opposite chirality that connect to the nanodot sides. The two segments of each pair are separated by a Bloch line, and the total chirality of each pair is zero. As the field magnitude is further increased, one of the domain walls pairs is rearranged to favor the nucleation of a cylindrical domain at the surface of the nanodot (state 7) with the abrupt propagation into the interior of the sample. Two other pairs of domain walls disappear from the sample. The cylindrical bubble domain with two opposite-chirality domain walls settles at the center of the nanodot (state 8) and further transforms to the $R_+$ skyrmion with $\bar{\chi} > 0$ and negative mean polarization (state 9) through the collapse of the negative-chirality domain wall. As the field magnitude continues to increase, the core of the skyrmion shrinks and the system arrives at state 10. Finally, the system jumps to the uniformly polarized state 11 with the negative polarization orientation. If, however, the field at point 9 is reversed with a decreasing magnitude to reach zero (the violet branch in Fig. 2d), the system remains in the $R_+$ skyrmion state at $U = 0$. With the further field increase in the positive direction, the system repeats the sequence 3–11 but with the opposite polarization and chirality. In other words, the system finally returns into the uniform up-polarized state 1.

One can make the system evolve having started with the negative fields and the polarization poled in the negative direction. In this case, the system will follow the branches denoted by the red-yellow traces in Fig. 2d. The emerging hysteresis branches are symmetric to the blue-violet ones with respect to $U \to -U$ reversal. The red branch corresponds to the blue one and the yellow branch corresponds to the violet one. The corresponding polarization states for the potential $U$ are obtained from their counterparts from the blue-violet branches at the potential $-U$ by the reversing the sign of the $P_z$ component of the polarization **P**. By the proper sweep protocol of the electric field, one arrives therefore to the $R_-$ and $L_-$ skyrmions at $U = 0$. Therefore, one sees that by the appropriate set of the protocol, one can obtain and switch between all of the four skyrmion states, $L_\pm$ and $R_\pm$ with different chirality and polarity orientations.

**The spontaneous chirality breaking.** It is important to note that the applied electric field does not possess its own chirality. This implies that the direction of the switch is determined rather by underlying local fluctuations in the chirality of the material which then serve as the nucleation centers of the emergent skyrmions. We describe this effect by introducing the fluctuating chirality field $\Lambda(r)$, which enters the system energy functional as the additional term $-\Lambda(\mathbf{P} \cdot \nabla \times \mathbf{P})$.

On the verge of the spontaneous symmetry breaking, even the slightest fluctuations would push the system into the either of degenerate, “left” and/or “right” free energy minima. In nature, one or another type of fluctuations would arise spontaneously but randomly. It is natural, in the case of numerical simulations, to quantify and control this effect, which is achieved by introducing the conjugate chiral field $\Lambda$ which impersonates random fluctuations of chirality that would appear in a real experimental system. In our numerical experiments, these spatial fluctuations are implemented via generating random tetrahedral configurations, maintaining the approximately constant mesh size. Altering a particular mesh changes the sign of emerging 1–2 and/or 6–8 jumps at the blue branch. The triggering effect of mesh fluctuations is verified by the mirror reflection of the discretization mesh leading to changing the sign of the chirality jump to the opposite one. Remarkably, even small fluctuations in $\Lambda$ lead to switching chirality. This implies an opportunity for laser-activated manipulation of the polarization³⁶
employing the circular polarized irradiation of the optical tweezers for controlling the direction of the chirality switch.

## Methods

**Functional and coefficients.** The polarization of the strained nanodot is obtained from the minimization of the free energy functional, depending on the polarization, $\mathbf{P} = (P_1, P_2, P_3)$, and the electrostatic potential, $\varphi$,

$$F = \int \left( [a_i^*(u_m, T)P_i^2 + a_{ij}^*P_i^2P_j^2 + a_{ijk}^*P_i^2P_j^2P_k^2]_{i \le j \le k} + \frac{1}{2}G_{ijkl}(\partial_i P_j)(\partial_k P_l) + (\partial_i \varphi)P_i - \frac{1}{2}\varepsilon_0\varepsilon_b(\nabla \varphi)^2 - \Lambda(\mathbf{P} \cdot \nabla \times \mathbf{P}) \right) d^3r, \quad (2)$$

where the summation over the repeated indices $i, j, \dots = 1, 2, 3$ (or $x, y, z$) is performed. The numerical parameters are specified for the PbTiO$_3$ nanodot, strained by the substrate with the compressive misfit strain $u_m \simeq -0.002$. The first square brackets term of (2) stands for the Ginzburg-Landau energy of the strained ferroelectric film$^{37}$, written in the form given in$^{33}$. The 2nd-order coefficients depend on the misfit strain $u_m$ and temperature $T$ and are expressed as $a_1^* = a_2^* = 3.8 \times 10^5(T - 479^{\circ}\text{C}) - 11 \times 10^9 u_m \text{ C}^{-2} \text{ m}^2 \text{ N}^{-1}$ and $a_3^* = 3.8 \times 10^5(T - 479^{\circ}\text{C}) + 9.5 \times 10^9 u_m \text{ C}^{-2} \text{ m}^2 \text{ N}^{-1}$. The strained renormalized 4th-order coefficients partially account for the elastic interactions, obey the tetragonal symmetry conditions and equal to $a_{11}^* = a_{22}^* \simeq 0.42 \times 10^9 \text{ C}^{-4} \text{ m}^6 \text{ N}$, $a_{33}^* \simeq 0.05 \times 10^9 \text{ C}^{-4} \text{ m}^6 \text{ N}$, $a_{13}^* = a_{23}^* \simeq 0.45 \times 10^9 \text{ C}^{-4} \text{ m}^6 \text{ N}$ and $a_{12}^* \simeq 0.73 \times 10^9 \text{ C}^{-4} \text{ m}^6 \text{ N}$. The 6th-order coefficients conserve the cubic symmetry, $a_{111} = a_{222} = a_{333} \simeq 0.26 \times 10^9 \text{ C}^{-6} \text{ m}^{10} \text{ N}$, $a_{112} = a_{113} = a_{223} \simeq 0.61 \times 10^9 \text{ C}^{-6} \text{ m}^{10} \text{ N}$, and $a_{123} \simeq -3.7 \times 10^9 \text{ C}^{-6} \text{ m}^{10} \text{ N}$. The second term of (2) corresponds to the gradient energy. The gradient energy coefficients $G_{ijkl}$ are obtained by the cubic symmetry permutations of the non-equivalent representatives $G_{1111} = 2.77 \times 10^{-10} \text{ C}^{-2} \text{ m}^4 \text{ N}$, $G_{1122} = 0$, and $G_{1212} = 1.38 \times 10^{-10} \text{ C}^{-2} \text{ m}^4 \text{ N}^{38}$. The next two terms in (2) correspond to the electrostatic energy, written in terms of the electrostatic potential $\varphi$$^{39}$. Here, $\varepsilon_0 = 8.85 \times 10^{-12} \text{ C V}^{-1} \text{ m}^{-1}$ is the vacuum permittivity, and $\varepsilon_b \simeq 10$ is the background dielectric constant of the non-polar ions$^{40}$. The last term of (2) emulates the interaction of the ferroelectric polarization with the material chirality fluctuation, described by the parameter $\Lambda$. In most simulations, we took $\Lambda = 0$. The local fluctuation, defining the direction of the chirality jumps naturally arise due to the random mesh configuration. To calibrate the effect of the mesh fluctuations we swept the value of $\Lambda$ and found that at the threshold value $\Lambda_c \simeq 7 \times 10^{-5} \text{ C}^{-2} \text{ m}^3 \text{ N}$, the direction of the chirality jump 2-3 at Fig. 2d changes to the opposite one.

**Computational techniques.** The simulations were performed using the FERRET package$^{40}$, designed for the multi-physics simulation environment MOOSE$^{41}$. In order to determine the minima of the free energy functional we use the standard technique of letting system to evolve ("relax") into these minima. This effective relaxation process is set by the time-dependent Ginzburg-Landau (TDGL) relaxation equation $-\gamma \partial P_i / \partial t = \delta F / \delta P_i$, where $F$ is the total static free energy functional (2) including the electrostatic potential $\varphi$. The latter is found at each step of the relaxation as a solution of the Poisson equation $\varepsilon_0 \varepsilon_b \nabla^2 \varphi = \partial_i P_i$. The TDGL equation includes the parameter $\gamma$ that sets a time scale for the dissipation of energy hence the rate of the motion of the system to a local free energy minimum. Since we are interested in the rapid computation of the sought local minima, but not in the real dynamics of the system, we choose $\gamma$ that provides the fast relaxation ensuring, at the same time, that the observation time scale remained longer than the relaxation time of the ferroelectric system. Typically, this is the case also in experimental situations since nanoscale ferroelectric systems can relax in a static field on times of the order of a few nanoseconds. Therefore, the TDGL equation gives indeed a good description of quasi-static properties, as has been demonstrated over the past few decades$^{25,31,42-45}$. Hence, when talking about the "relaxation" and "time scale," we refer to the standard so-called Phase Field Relaxation numerical method designed for finding the minima of the free energy of complex systems. In our calculations we thus set the relaxation parameter $\gamma$ equal to unity. For more details on simulations see Supplementary Information. The geometry and conceptual setup of the simulated system are shown in Fig. 1a. We selected the diameter and the thickness of the nanodot as 40 and 20 nm respectively as the optimal geometry for the skyrmion's observation. According to estimates the modification of the geometrical parameters over 15-20% does not change the results much, although the evolution of the Bloch points may somewhat differ marginally and in ways that are inconsequential for our results.

The driving field was controlled by the voltage $U$, applied to the electrodes. The upper electrode was separated from the nanodot by the empty, or low-$\varepsilon$ dielectric space of the thickness of 60 nm. The initial paraelectric state with the small randomly-distributed polarization was used as an initial condition for the quench to the original (virgin) state. Then, the quasi-static field variation protocols were used with the polarization distribution at the previous stage taken as the initial condition. The different finite-element meshes were used to ensure the stability of the process.

## Data availability

Computational scripts are available online at https://github.com/ferroelectrics/skyrmion.

Received: 20 February 2020; Accepted: 30 April 2020;

Published online: 26 May 2020
## References

1. Kelvin, W. T. B. *The molecular tactics of a crystal* (Clarendon Press, 1894).

2. Hegstrom, R. A. & Kondepudi, D. K. The handedness of the universe. *Scientific American* **262**, 108–115 (1990).

3. Wagnière, G. H. *On chirality and the universal asymmetry: reflections on image and mirror image* (John Wiley & Sons, 2008).

4. Berger, R. & Stohner, J. Parity violation. *Wiley Interdisciplinary Reviews: Computational Molecular Science* **9**, e1396 (2019).

5. Neufeld, O. & Cohen, O. Optical chirality in nonlinear optics: Application to high harmonic generation. *Phys. Rev. Lett.* **120**, 133206 (2018).

6. Naaman, R. & Waldeck, D. H. Spintronics and chirality: Spin selectivity in electron transport through chiral molecules. *Annual review of physical chemistry* **66**, 263–281 (2015).

7. Yoo, S. & Park, Q.-H. Metamaterials and chiral sensing: a review of fundamentals and applications. *Nanophotonics* **8**, 249–261 (2019).

8. Inaki, M., Liu, J. & Matsuno, K. Cell chirality: its origin and roles in left-right asymmetric development. *Philosophical Transactions of the Royal Society B: Biological Sciences* **371**, 20150403 (2015).

9. Bahr, C. & Kitzerow, H.-S. *Chirality in liquid crystals* (Springer, 2001).

10. Chandrasekhar, S. *Chirality in liquid crystals* (Springer Science & Business Media, 2006).

11. Seki, S. & Mochizuki, M. *Skyrmions in magnetic materials* (Springer, 2016).

12. Liu, J. P., Zhang, Z. & Zhao, G. *Skyrmions: topological structures, properties, and applications* (CRC Press, 2016).

13. Zhang, S. *Chiral and topological nature of magnetic skyrmions* (Springer, 2018).

14. Phatak, C., Heinonen, O., De Graef, M. & Petford-Long, A. Nanoscale skyrmions in a nonchiral metallic multiferroic: Ni₂MnGa. *Nano letters* **16**, 4141–4148 (2016).

15. Navas, D. *et al*. Route to form skyrmions in soft magnetic films. *APL Materials* **7**, 081114 (2019).

16. Tagantsev, A. K., Courtens, E. & Arzel, L. Prediction of a low-temperature ferroelectric instability in antiphase domain boundaries of strontium titanate. *Physical Review B* **64**, 224107 (2001).

17. Lee, D. *et al*. Mixed Bloch-Néel-Ising character of 180° ferroelectric domain walls. *Phys. Rev. B* **80**, 060102 (2009).

18. Wojdel, J. C. & Iñiguez, J. Ferroelectric transitions at ferroelectric domain walls found from first principles. *Phys. Rev. Lett.* **112**, 247603 (2014).

19. Cherifi-Hertel, S. *et al*. Non-ising and chiral ferroelectric domain walls revealed by nonlinear optical microscopy. *Nature communications* **8**, 15768 (2017).

20. Baudry, L., Luk'yanchuk, I. A. & Sené, A. Inhomogeneous polarization switching in finite-size cubic ferroelectrics. *Ferroelectrics* **427**, 34–40 (2012).

21. Baudry, L., Luk'yanchuk, I. A. & Sené, A. Switching properties of nano-scale multi-axial ferroelectrics: geometry and interface effects. *Integr. Ferroelectr.* **133**, 96–102 (2012).

22. Baudry, L., Sené, A., Luk'yanchuk, I. A., Lahoche, L. & Scott, J. F. Polarization vortex domains induced by switching electric field in ferroelectric films with circular electrodes. *Phys. Rev. B* **90**, 024102 (2014).

23. Nahas, Y. *et al*. Discovery of stable skyrmionic state in ferroelectric nanocomposites. *Nat. Commun.* **6** (2015).

24. Gonçalves, M. P., Escorihuela-Sayalero, C., Garca-Fernández, P., Junquera, J. & Iñiguez, J. Theoretical guidelines to create and tune electric skyrmion bubbles. *Science advances* **5**, eaau7023 (2019).

25. Das, S. *et al*. Observation of room-temperature polar skyrmions. *Nature* **568**, 368–372 (2019).

26. Luk'yanchuk, I., Tikhonov, Y., Razumnaya, A. & Vinokur, V. Hopfions emerge in ferroelectrics. *arXiv preprint arXiv:1907.03866* (2019).

27. Bratkovsky, A. & Levanyuk, A. Abrupt appearance of the domain pattern and fatigue of thin ferroelectric films. *Phys. Rev. Lett.* **84**, 3177 (2000).

28. Stephanovich, V., Luk'yanchuk, I. & Karkut, M. Domain-enhanced interlayer coupling in ferroelectric/paraelectric superlattices. *Phys. Rev. Lett.* **94**, 047601 (2005).

29. De Guerville, F., Luk'yanchuk, I., Lahoche, L. & El Marssi, M. Modeling of ferroelectric domains in thin films and superlattices. *Mater. Sci. Eng. B* **120**, 16–20 (2005).

30. Zubko, P., Stucki, N., Lichtensteiger, C. & Triscone, J.-M. X-ray diffraction studies of 180° ferroelectric domains in PbTiO³/SrTiO³ superlattices under an applied electric field. *Phys. Rev. Lett.* **104**, 187601 (2010).

31. Yadav, A. *et al*. Observation of polar vortices in oxide superlattices. *Nature* **530**, 198–201 (2016).

32. Martelli, P.-W., Mefire, S. M. & Luk'yanchuk, I. A. Multidomain switching in the ferroelectric nanodots. *EPL (Europhysics Letters)* **111**, 50001 (2015).

33. Baudry, L., Luk'yanchuk, I. & Vinokur, V. M. Ferroelectric symmetry-protected multibit memory cell. *Scientific reports* **7**, 42196 (2017).

34. Malozemoff, A. & Slonczewski, J. *J. Magnetic Domain Walls in Bubble Materials: Advances in Materials and Device Research*, vol. 1 (Academic press, 2016).

35. Salje, E. & Scott, J. Ferroelectric bloch-line switching: A paradigm for memory devices? *Appl. Phys. Lett.* **105**, 252904 (2014).

36. Stoica, V. *et al*. Optical creation of a supercrystal with three-dimensional nanoscale periodicity. *Nature materials* **18**, 377–383 (2019).

37. Pertsev, N., Zembilgotov, A. & Tagantsev, A. Effect of mechanical boundary conditions on phase diagrams of epitaxial ferroelectric thin films. *Phys. Rev. Lett.* **80**, 1988 (1998).

38. Wang, J., Shi, S.-Q., Chen, L.-Q., Li, Y. & Zhang, T.-Y. Phase-field simulations of ferroelectric/ferroelastic polarization switching. *Acta Materialia* **52**, 749–764 (2004).

39. Landau, L., Lifshitz, E. & Pitaevskii, L. *Eloctrodynamics of Continuous Media* (Oxford: Pergamon, 1984).

40. Mangeri, J. *et al*. Topological phase transformations and intrinsic size effects in ferroelectric nanoparticles. *Nanoscale* **9**, 1616–1624 (2017).

41. Gaston, D., Newman, C., Hansen, G. & Lebrun-Grandie, D. Moose: A parallel computational framework for coupled systems of nonlinear equations. *Nucl. Eng. and Design* **239**, 1768–1778 (2009).

42. Li, Y., Hu, S., Liu, Z. & Chen, L. Effect of substrate constraint on the stability and evolution of ferroelectric domain structures in thin films. *Acta materialia* **50**, 395–411 (2002).

43. Choudhury, S., Li, Y., Krill Iii, C. & Chen, L. Effect of grain orientation and grain size on ferroelectric domain switching and evolution: Phase field simulations. *Acta materialia* **55**, 1415–1426 (2007).

44. Chen, L.-Q. Phase-field method of phase transitions/domain structures in ferroelectric thin films: a review. *Journal of the American Ceramic Society* **91**, 1835–1844 (2008).

45. Zhang, J., Schlom, D., Chen, L. & Eom, C. Tuning the remanent polarization of epitaxial ferroelectric thin films with strain. *Applied Physics Letters* **95**, 122904 (2009).

## Acknowledgements

This work was supported by the US Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division (O.H. V.M.V. and partly I.L. and A.G.); by the H2020 RISE-ENGIMA, RISE-MELON and ITN-MANIC actions (Y.T., A.R. and I.L.), By Bourse Ostrogradski of French Government (Y.T.), and by the Southern Federal University, Russia (A.R. and Y.T.).
**Author contributions**

I.L., O.H. and V.M.V. conceived the work, S.K., A.S., O.H., I.L., and V.M.V. performed analytical calculations, J.M., O.H., and S.N. developed the Ferret code, Y.T. elaborated on the Ferret code adding new functionalities, Y.T., S.K., J.M., M.P., L.B., A.G., A.R. performed numerical calculations, Y.T. and A.R. visualized numerical data, Y.T., S.N., A.R., O.H., I.L. and V.M. analyzed the results, I.L., and V.M.V. wrote the manuscript, all the authors discussed the manuscript.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41598-020-65291-8.

**Correspondence** and requests for materials should be addressed to V.M.V.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_5.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020