Discovery of stable skyrmionic state in ferroelectric nanocomposites

Non-coplanar swirling field textures, or skyrmions, are now widely recognized as objects of both fundamental interest and technological relevance. So far, skyrmions were amply investigated in magnets, where due to the presence of chiral interactions, these topological objects were found to be intrinsically stabilized. Ferroelectrics on the other hand, lacking such chiral interactions, were somewhat left aside in this quest. Here we demonstrate, via the use of a first-principles-based framework, that skyrmionic configuration of polarization can be extrinsically stabilized in ferroelectric nanocomposites. The interplay between the considered confined geometry and the dipolar interaction underlying the ferroelectric phase instability induces skyrmionic configurations. The topological structure of the obtained electrical skyrmion can be mapped onto the topology of domain-wall junctions. Furthermore, the stabilized electrical skyrmion can be as small as a few nanometers, thus revealing prospective skyrmion-based applications of ferroelectric nanocomposites.

T he topological invariance of nontrivial field structures such as vortices and skyrmionic configurations against fluctuations and deformations triggers considerable interest for device applications 1,2 . Many advances have been achieved in predicting their formation and characterizing the emergent order they subtend in magnetic systems 3,4 . Whereas vortex configurations, or related flux-closure configurations, have already been revealed in ferroelectrics (FE) within nanoscale geometries (see, for example, refs 5-13 and references therein), the search for a spontaneous formation of skyrmionic textures in FEs was somewhat hindered by the absence of intrinsic chiral interactions that are known to stabilize such configurations in non-centrosymmetric magnetic systems [14][15][16] . However, it is now unequivocal that such interactions are not a prerequisite for obtaining skyrmionic topological patterns 17 , as non-coplanar swirling field structures have been amply predicted and observed in magnets that do not meet the conditions typically deemed necessary 18-20 . Recognizing in this lifted restriction a breach for probing skyrmions in ferroelectric materials, we build on the topology of the polarization field in confined geometries and show here, via the use of a first-principles-based technique, that nanoscale structures can readily become the locus of skyrmionic configuration of polarization. We also demonstrate that the field-induced skyrmionic texture is thermodynamically stable against electric and thermal perturbations. We further find that its topological charge density markedly differs from that of axisymmetric magnetic skyrmion textures in that it is broken into equally contributing fractions anchored at domain walls junctions. Our finding therefore unravels the interplay between geometry and topology in stabilizing extrinsically protected skyrmionic configuration, and brings to the fore the possibilities of extending skyrmion-based devices to ferroelectrics.

Results
Confined nanocomposite geometry. We choose to focus on a ferroelectric nanocomposite consisting of a cylindrical BaTiO 3 (BTO) nanowire with a radius of 2.7 nm embedded in a SrTiO 3 (STO) matrix in our search for electrical skyrmion. The choice of such a nanostructure is motivated by the advancement in the controlled growth of composites with tailored functionalities 21 , and the wide variety of novel behaviours they feature (note, however, that the growth of low-dimensional perovskites inside a perovskite matrix remains a difficult task to experimentally achieve). A primarily relevant topological feature of this type of nanocomposites is the occurrence of a size-induced chiral symmetry breaking, which stabilizes a vortex domain structure coexisting with a spontaneous polarization along the axial direction of the nanowire 10 . The investigated ferroelectric structure is schematized in Fig. 1a. This nanocomposite is mimicked by a 36 Â 36 Â 6 supercell that is periodic along the x, y and z axis (which lie along the pseudocubic [100], [010] [001] directions, respectively). Its properties are predicted by performing Monte Carlo simulations of the first-principlesbased effective Hamiltonian scheme of ref. 22.
Electric field treatment. The creation of an electrical skyrmion is presently achieved by a numerical procedure consisting of the following few steps. We first perform a temperature annealing under an external electric field, E [001] ¼ 10 8 V m À 1 , applied along the pseudocubic [001] direction. On reaching 15 K, we set the field to zero and further relax the priorly obtained low-temperature configuration. The resulting relaxed configuration features a spontaneous polarization along the axial direction of the nanowire, co-occurring with a flux-closure four-domain vortex structure of the cross-sectional polarization field, in-plane pattern whereby the strength of the depolarization field is reduced 23 . More precisely, we find that the (1-3) Newnham's connectivity 24 (BTO phase is one-dimensionally connected and the STO phase is three-dimensionally connected) of the considered nanocomposite structure hosts a polarization state possessing translational invariance along the axial direction of the wire. The polarization field thus depends only on x and y spatial coordinates, enabling a visualization of the system as consisting of two-dimensional layers (z-planes). The corresponding crosssectional polarization field is shown in Fig. 1b, which displays the x and y components of the electric dipoles in an arbitrary (001)plane of the shifted periodic supercell. This state will be referred to as (V xy |FE z ) in the following, to emphasize that it has a vortex state in the z-planes but also possesses an electrical polarization along the [001] direction. The distribution of the latter as a function of x and y spatial coordinates is such that that local dipoles within the wire have a z-component around 30% larger in magnitude than that of the local dipoles located in the matrix (B30.5 mC cm À 2 versus B22.4 mC cm À 2 , respectively). Moreover, in this (V xy |FE z ) state, in addition to the central vortex whose core is confined within the nanowire, the matrix exhibits a vortex at mid-way between second-nearest neighbour wires as well as two anti-vortices at mid-way between first-nearest neighbour wires of the periodic supercell. These punctual topological defects are associated with singularities in the two-dimensional cross-sectional polarization field and are found to be anchored at the junction of DW x and DW y domain walls, which interpolate between domains of distinct in-plane (x and y) components of polarization (Fig. 1b). The sum of the O(2) winding numbers, each defined as a line integral measuring the change in angle of the two-dimensional cross-sectional polarization field over a closed path and giving n ¼ þ 1 for vortices and n ¼ À 1 for anti-vortices 25 , yields a zero net topological charge, in accordance with Poincaré-Hopf theorem 26 .
In our search for skyrmionic configuration, we then proceed with further subjecting the (V xy |FE z ) state to an E 00 1 ½ electric field, that is, to a field pointing oppositely to (V xy |FE z )'s polarization. The evolution of the z-component of polarization, , one can see that P z first slightly decreases before exhibiting a steep drop. Microscopical insight into the behaviour of the overall polarization is given by Fig. 2b-d, which show the distribution of P z within an arbitrary z-plane of the supercell for three different values of E 00 1 ½ . It is therein seen that as E 00 1 ½ increases, the matrix transiently develops inhomogeneity and is thus primarily accountable for the observed reduction of P z . After completion of this transient process, at a certain threshold value E* of the field (indicated by a vertical line in Fig. 2a), the dipolar configuration is such that the z-component of the electric dipoles of the matrix is anti-parallel to that of the dipoles comprised within the wire.
As a result of the wire and the matrix having different switching fields, this final electric field treatment has thus prompted the nucleation of a p-wall (DW z ) positioned at the interface between the wire and matrix. Let us note that for values of E 00 1 ½ exceeding E* and rising up to À 3 Â 10 8 V m À 1 , the z-component of the dipoles located in the matrix slightly increases in magnitude, while the anti-aligned z-component of the dipoles within the wire remains quasi constant. Remarkably and as shown in Fig. 4a,b, the ensuing dipolar configuration at E* forms an electrical skyrmion (to be denoted as Sk) with a polarization antiparallel to E* at its center and parallel to E* at its periphery. Such skyrmionic configuration appearing at E* persists up to À 3 Â 10 8 V m À 1 and bears a radial p-modulation of the z-component of polarization, as well as an azimuthal modulation of the x and y components. This double modulation is characteristic of a chiral, or Bloch-like skyrmion 27 . Notably, the removal of the applied field E 00 1 ½ does not alter the morphology of the (Sk) configuration, which indicates the stability of the obtained skyrmionic state (Sk). Such a finding features resemblance with magnetic systems within which occur magnetic-field-induced skyrmions that can persist after releasing the external field 17,28,29 .
Topological characterization. As a result of the radial p-modulation involving the z-component of the polarization field, an O(3) topological invariant, the Pontryagin density q, is required to unambiguously identify the topological nature of the obtained three-dimensional swirling field structure 17 . This invariant is , where the three-component unit vector u denotes the normalized local dipole moment. For a skyrmion, its surface integral over an arbitrary (001)-plane is an integer topological charge Q yielding the number of times u wraps the unit sphere 17 . Let us note that the topological density q can be expressed as a total derivative, and is thus intrinsically insensitive to the polarization texture details. Figure 3 provides the calculated topological density distribution with spatial derivatives being approximated using a fourth-order accurate finite difference scheme on interpolated data.
Remarkably, Fig. 3 reveals the presence of several narrow peaks located at the junctions of (DW x and DW y ), (DW z and DW x ) and (DW z and DW y ) domain walls, and equally contributing to a topological charge Q ¼ 1 (note that, by contrast, this invariant is trivial for (V xy |FE z ), for which Q ¼ 0, as it is seen in Fig. 2e). This finding, which we further cross-checked through the method of topological charge computation of ref. 30, attests to the topological nature of the obtained skyrmionic configuration. It also demonstrates that its broken topological structure can be mapped onto the topology of domain-wall junctions, as a genuine consequence of the considered nanocomposite geometry. A distinctive feature of the electrical (Sk) is thus its four-fold symmetric topological charge density, markedly differing from that of axisymmetric magnetic skyrmion textures and other topological objects. Moreover, the geometry of the system enables obtaining an electrical skyrmion whose size can be as small as a few nanometers whereas magnetic skyrmions non-arising from Dzyaloshinskii-Moriya interaction extend on significantly larger scales 17,18,29 , thereby widening the scope and capabilities of future skyrmion-based applications.
Assessment of the skyrmionic state's stability. Furthermore and as above demonstrated, the (Sk) state is at the least metastable. To bridge the topologically distinct (V xy |FE z ) and (Sk) states and estimate their relative stability, we have performed a Metropolis Monte Carlo sampling of the transition paths connecting them in configurational space. Different Markov chains allowed to gather statistics on the pathways accessible under the application of E* on the (V xy |FE z ) configuration. Averaging the energies associated with different trajectories parametrized by the P z component of polarization yielded the internal energy profile shown in Fig. 4c. It appears clearly that the obtained (Sk) is metastable with a local minimum only 7 mHa higher than that of the (V xy |FE z ) state. While both the (V xy |FE z ) and (Sk) states are stabilized by long-range dipolar interactions, we find that the difference in their energy mainly originates from the p-wall energy formation cost encoded in the short-range interaction. The topological stability of the (Sk) state is reflected in the finite energy barrier which protects it against spontaneous decay. This energy barrier (B0.02 mHa) underpins the first-order nature of the transition associated with the abrupt partial switching process (Fig. 2a), in that the skyrmionic structure cannot be continuously deformed into that of the (V xy |FE z ) state.
We also probe the thermal stability of the (Sk) state by increasing the temperature. At moderate temperatures, thermal fluctuations do not destabilize the (Sk) state, which remains robust up to 65 K. Specifically, thermally activated excitations in the form of confined and self-compensating O(2) punctual defects (vortex-antivortex bound pairs) induce homotopical configurations preserving the topological texture of the (Sk) state. This finding is consistent with the range of stability shown in Fig. 4d, which provides the temperature dependence of the estimated threshold external electric field E* needed to drive the system towards the skyrmion state. It is seen that E* decreases in magnitude as the temperature is increased. The transition line between the (Sk) and the (V xy |FE z ) states terminates at a finite temperature and field (65 K and À 0.24 Â 10 7 V m À 1 ), on which the topological protection is no longer in effect. At this termination point, thermal and electrical excitations endow the system with a sufficiently high energy to overcome the topological barrier. The (Sk) thermal stability range can be enhanced by substituting the BTO material composing the wire with a higher-T C ferroelectric, while conserving a discrepancy between the wire and matrix polarizabilities.
We further appraised the skyrmion stability by subjecting it to a gradually increasing E [111] electrical field oriented along the [111] direction at 15 K. While under a certain threshold of E [111] the excited configuration relaxes back to the initial (Sk) state, for strong enough magnitude, the induced nonequilibrium state converts into an homogeneous polar rhombohedral state (to be denoted by (Rh)) of higher energy ( À 1.45 mHa). This metastable (Rh) state consists of electric dipoles pointing along the [111] direction in both the wire and the matrix. We find that the energy barrier associated with this annihilation pathway (B0.05 mHa) is higher than the one connecting the (Sk) and the (V xy |FE z ) states. The reason behind the difference in height of the barriers can be approached through topological considerations; while a transition from the (Sk) state to the topologically trivial (Rh) state requires disentangling the former by unwinding its two modulations, the transition from the (Sk) state to the (V xy |FE z ) state involves releasing the p-modulation, thereby only affecting the z-component of polarization while conserving the crosssectional configuration of polarization and its domain-wallpinned O(2) punctual defects.

Discussion
Inquiring into the dependence of the topological (Sk) texture on the geometry of the structure, we performed additional simulations where the effect of the radius of the nanowire, its shape and the length n z along the [001] pseudocubic direction of the supercell were assessed. We found that the results are not affected by neither n z (over the range n z ¼ 6 to n z ¼ 10) nor the shape (square versus circle cross-sections) of the nanowire. Varying the size of the wire yielded a non-monotonic dependence of E* on the radius R of the cylindrical wire. With increasing radii of BaTiO 3 nanowires, E* was found to increase, merely indicating the relatively empowered influence of the wire on the matrix as a result of its enlargement. Moreover, we found that considering a rectangular array of cylindrical wires within the supercell endowed the resulting skyrmion lattice with helicity Z 2 ð Þ degrees of freedom (the circulation direction of the cross-sectional polarization around the skyrmion core); depending on the inter-wire distances, the interaction between skyrmions led to either a single-chirality (phase-locking) pattern or to an alternate arrangement of skyrmionic chiralities within this lattice. We also note that skyrmion-like configurations should exist in other geometries than the nanocomposite structure herein investigated, provided that similar physics comes into play 31 . Moreover, free charges and structural defects can screen or modify depolarizing field and thus affect the formation of electrical skyrmions. However, we are confident that these latter topological defects will soon be observed based on the fact that other depolarizing-field-related dipolar patterns have already been experimentally seen in some ferroelectric systems (see, for example, refs 9-13).
In conclusion, we have shown that ferroelectric nanocomposites can accommodate a skyrmionic configuration as a stable state of polarization. Our study brings to the fore the interplay between geometry and topology in providing extrinsic topological protection. Finally, analysing the (Sk) topological properties, we found that its topological charge density is broken into fractions pinned at domain walls junctions, thereby singling out electrical skyrmions from topological objects with axisymmetric topological charge densities.

Methods
Effective Hamiltonian. The effective Hamiltonian employed in this study is that of ref. 22. Its total internal energy has the form: The first energy term, E ave , depends on the local soft mode u i centred on the Ti-sites of the 5-atom unit cell i and directly proportional to the electric dipole moment of the corresponding unit cell, the dimensionless displacement variables {v i } defined at the cell corners and entering in the calculation of the inhomogeneous strain tensor components of the cell i 32,33 , and the homogeneous strain tensor {Z H }, which allows the simulations to account for the change in size and shape of the supercell 32,33 . E ave consists of five parts: a local-mode self-energy, a long-range dipole-dipole interaction, a short-range interaction between soft modes, an elastic energy, and an interaction between the local modes and local strain 32,33 . Its parameters are fitted from first-principles calculations performed on a uniform virtual hAiTiO 3 system that averages the potentials of the pure parent compounds 34,35 to model (Ba 0.5 Sr 0.5 )TiO 3 solid solutions. The second energy term, E loc , involves the {s j } parameters and {Z loc }, which respectively correspond to the set of variables {s j } characterizing the atomic distribution of the mixed A-sublattice 36 , with s j ¼ þ 1 or À 1 corresponding to the presence of either Ba or Sr atom at the A-lattice site j, and the local strain {Z loc } stemming from the difference in ionic radii between Ba and Sr atoms (which is ' 2 % ). This second energy term E loc stands as a perturbative correction to the virtual crystal approximation, as it accounts for the real nature of the A-site atom and its effect on the local soft modes and the inhomogenuous strain tensor, as well as it includes corrections pertaining to the strain induced by the size difference between Ba and Sr ions 22 . The parameters entering in E loc are derived from first principles. The nanocomposite herein considered is stress-free, which allows the matrix, interface and wires to adopt different lattice constants, thereby modelling a (realistic) inhomogeneous strain. The Monte Carlo simulations were performed using up to 50,000 sweeps. It is worth mentioning that the above described first-principles-based effective Hamiltonian scheme 22   (Ba,Sr)TiO 3 materials (BST) 39 . For instance, our effective Hamiltonian gives transition temperatures of 385 K (from cubic to tetragonal), 280 K (from tetragonal to orthorhombic) and 230 K (from orthorhombic to rhombohedral) in BaTiO 3 bulk, which compares rather well with the experimental values of 400, 280 and 180 K 22 . Other examples demonstrating the accuracy of this effective Hamiltonian is the subtle temperature-gradient-induced polarization 40 , and the existence of two modes (rather than a single one as previously believed for a long time) contributing to the GHz-THz dielectric response of BST compounds 39 .
Notably, nanocomposites such as the one considered in the present study and consisting of BTO nanowires embedded in a STO matrix, were predicted to exhibit a spatially modulated field structure from which a variety of striking phenomena [41][42][43] originates. Let us note that the simulated matrix does not rigorously correspond to pure STO but rather coincides with Ba 1 À x Sr x TiO 3 having 85% of strontium because of some limitations related to the employed effective Hamiltonian (in particular, this effective Hamiltonian predicts a paraelectric-to-ferroelectic transition of pure STO at a Curie temperature that experimentally corresponds to the one of Ba 0.15 Sr 0.85 TiO 3 ).
A final remark pertains to the magnitude of electric fields reported in this study. Generally, effective Hamiltonian schemes have a tendency to overestimate electric experimental fields by a factor of about 25 with respect to experiments 44 . However, a field of 10 9 V m À 1 has been applied to some nano-structures in recent measurements 45 .