Multilayer structure

Ta(4.7 nm)/[Pt(4 nm)/Co(1.6 nm)/Ta(1.9 nm)] × 5 (hereafter Pt/Co/Ta) multilayers were deposited by using ultrahigh vacuum magnetron sputtering on PMN-PT(001) single-crystal substrate (as shown in Fig. 1a), and separately on a 30-nm-thick Si₃N₄ membrane for Lorentz transmission electron microscopy (L-TEM) observation. The stack structure and setup for MFM measurement with in situ electric fields are schematically shown in Supplementary Fig. 1. The bottom Ta layer was used to flat the surface of substrate and promote the interface quality¹⁴, which enables a large PMA in the multilayer and a strong interfacial DMI at the Pt/Co interface³¹. Pt/Co and Ta/Co interfaces generate strong interfacial DMI and very weak interfacial DMI (see ref. ⁸), respectively, leading to a considerable net interfacial DMI, which is beneficial for stabilizing Néel-type skyrmions. All experiments were performed at room temperature. Figure 1b shows the high angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) image of the sample, revealing flat and distinct interfaces. Energy dispersive X-ray (EDX) mapping further confirms the periodical Pt/Co/Ta tri-layer structure. Figure 1c shows the out-of-plane magnetic hysteresis loop (in-plane one in Supplementary Fig. 2).

Fig. 1

Structure and skyrmion characterizations.

a Schematic of the sample configuration. b Cross-sectional HAADF-STEM image of the multilayer and corresponding EDX mapping for element distributions of Pt (green), Co (red) and Ta (blue) with a scale bar of 5 nm. c Out-of-plane magnetic hysteresis. d–f Magnetic-domain evolutions with increasing magnetic field with the same scale bar of 1 μm. g Line-cut profile along the yellow dot line in the inset with Gaussian fitting (red curve), the inset is magnified individual skyrmion with a scale bar of 100 nm. h Representative L-TEM image of skyrmions with a scale bar of 1 μm.

Microscopic imaging of skyrmions

Due to its high-spatial resolution and the ease of operation, MFM has been widely used for studying skyrmions^9,19,32. Figure 1d–f shows the magnetic-domain structures with an out-of-plane bias magnetic fields of B = 0, 48, and 100 mT, respectively, where the blue (red) contrast indicates magnetization along the positive (negative) z direction. At B = 0 mT, the sample shows a labyrinthine domain pattern (Fig. 1d), consistent with the vanishing remnant magnetization of the multilayer (Fig. 1c). Upon applying magnetic fields, the labyrinth domain transforms into skyrmions (e.g., Fig. 1e), which shrink with increasing magnetic field. At the saturation magnetic field (100 mT), the skyrmions transform into an FM single domain state (Fig. 1f). A full set of MFM images are shown in Supplementary Fig. 3. A Gaussian fitting of a typical skyrmion profile (Fig. 1g), extracted from the inset of Fig. 1g along the yellow dot line, gives a skyrmion size/diameter of 214 nm comparable to the previous reports^8,14,32,33. Furthermore, the state of multiple isolated skyrmion was also observed (Fig. 1h) using L-TEM for sample grown on Si₃N₄ membrane (other images in Supplementary Fig. 4). Both the MFM and L-TEM imaging demonstrate the existence of skyrmions with small differences in skyrmion density and size due to different substrates.

Variations of strain, interfacial DMI, and magnetic anisotropy with electric field

Variations of strain, interfacial DMI, and magnetic anisotropy with electric field are then investigated. To obtain the strains in PMN-PT induced by electric fields, we performed in situ X-ray diffraction (XRD) measurements under different electric fields to get lattice constants c(E), calculated using (002) peak by Bragg’s law (raw data in Supplementary Fig. 5). Electric-field-induced out-of-plane strain (Fig. 2a), defined as [c(E) − c(+0 kV/cm)]/c(+0 kV/cm), can be used to calculate the in-plane compressive strain (Supplementary Fig. 6) based on the crystal symmetry of (001)-poled PMN-PT (see ref. ³⁴). The out-of-plane strain is largely volatile when decreasing electric field and asymmetric for positive and negative electric fields. Moreover, the strain is close to 0.4% at −4 kV/cm, much larger than the typical value of 0.1% (see ref. ³⁴). More discussion is shown in Supplementary Fig. 5. The strain and magnetization under different electric fields were measured for thousands of times in previous studies^35,36, which demonstrated the stability and durability of our devices.

Fig. 2

Variations of strain, interfacial DMI, and magnetic anisotropy.

a Out-of-plane strain variations of PMN-PT(001) substrate. b Wave-vector dependence of Δf under different electric fields with the interfacial DMI values in the inset with the error bar obtained from the standard error of Lorentzian fitting. c Angle-dependent FMR resonance field H_r(θ) and corresponding Kittel formula fitting (solid lines). d K_eff versus electric-field curve.

Interfacial DMI is essential for stabilizing skyrmions in Pt/Co/Ta (see ref. ⁸), and we used Brillouin light spectroscopy (BLS)^21,22 to measure it under different electric fields. By measuring the wave-vector dependence of the frequency difference Δf between the anti-Stokes and Stokes peaks, the interfacial DMI strength of the Pt/Co/Ta multilayer (D) can be determined by the following equation:

Besides the interfacial DMI, the effective PMA of the multilayer is also critical to the control of skyrmion^41,42. To determine the variation of the effective PMA with electric fields, we performed angular dependent FM resonance (FMR) measurements under different electric fields (see detailed experimental configuration in Supplementary Fig. 9). For each angle θ, a resonance field H_r(θ) can be determined from the FMR spectrum, as shown in Fig. 2c. The effective PMA (denoted as K_eff) can then be determined by fitting H_r(θ) with the Kittel formula (see text accompanying Supplementary Fig. 9), and the fitting results are plotted by the lines in Fig. 2c. The variation of the K_eff with electric field is shown in Fig. 2d, which resembles the variation of the interfacial DMI strength (inset of Fig. 2b) and is consistent with the strain variation (Fig. 2a). Thus, it can be concluded that strain transferred from the PMN-PT alters both the interfacial DMI strength of the Pt/Co interface and the effective PMA of the Pt/Co/Ta multilayer.

Nonvolatile electric-field-induced creation of multiple skyrmions

For investigating electric-field control of skyrmions using MFM, pretreatment procedures have been done for the sample. First, the sample was poled by varying electric fields between +4 kV/cm and −4 kV/cm for two cycles with a step of 0.2 kV/cm. Electric fields of +4 kV/cm or −4 kV/cm were then removed, yielding two remanent polarization states, denoted as +0 kV/cm and −0 kV/cm, respectively. After then, the sample was magnetically poled to saturation, followed by reducing the bias magnetic field B_bias to a lower value for subsequent in situ MFM measurement.

Figure 3 shows the electric-field-induced creation of skyrmions with B_bias = 60 mT. Initially, at E = +0 kV/cm, the sample is almost an FM single domain (Fig. 3a), and skyrmions appear (Fig. 3b) when applying E = −4 kV/cm. It has been suggested that the stray magnetic field from the scanning magnetic tip can alone lead to the creation of skyrmions³², but this possibility has been ruled out based on our control experiment (Supplementary Fig. 10). Notably, the skyrmion state remains after removing −4 kV/cm (Fig. 3c, at which E = −0 kV/cm). When further applying E = +4 kV/cm, the skyrmion state barely changes (Fig. 3d). Similar behaviors of skyrmion creation were also observed in sample starting from E = −0 kV/cm (Supplementary Fig. 11), that is, the creation of skyrmions can only occur when applying E = −4 kV/cm, under which the electric-field-induced strain from the PMN-PT is much larger than that the case of +4 kV/cm (see Fig. 2a). It was shown that the MFM images in Fig. 3a–d were taken from the same location on the sample (Supplementary Fig. 12). Moreover, the creation of skyrmions under other B_bias are shown in Supplementary Fig. 13.

Fig. 3

Skyrmion creation.

MFM images at E = +0 kV/cm (a), −4 kV/cm (b), −0 kV/cm (c), and +4 kV/cm (d) with B_bias = 60 mT. Corresponding simulation results of strain-mediated skyrmion creation with $$ {\epsilon }_{\left[110\right]}={\epsilon }_{\left[-110\right]}=0$$ (initial state) and D = 0.772 mJ/m² (e), $$ {\epsilon }_{\left[110\right]}={\epsilon }_{\left[-110\right]}=-0.189\%$$ and D = 0.585 mJ/m² (f), $$ {\epsilon }_{\left[110\right]}={\epsilon }_{\left[-110\right]}=-0.034\%$$ and D = 0.685 mJ/m² (g), and $$ {\epsilon }_{\left[110\right]}={\epsilon }_{\left[-110\right]}=-0.010\%$$ and D =  0.727 mJ/m² (h), with the blue and red contrasts corresponding to magnetizations pointing up and down, respectively. The scale bar is 1 μm. Evolutions of in-plane biaxial compressive strain $$ {\epsilon }_{\mathrm{in-plane}}={\epsilon }_{\left[110\right]}={\epsilon }_{\left[-110\right]}$$ (i), magnetoelastic energy density f_mel (j), and intrinsic energy density f_intrin (k). The time stages are t₀ (E = +0 kV/cm), t₁ (E = −4 kV/cm), t₂ (E = −0 kV/cm), and t₃ (E = + 4 kV/cm), respectively. The insets of k are the corresponding individual skyrmion evolutions with a scale bar of 50 nm.

We performed micromagnetic simulations to elucidate the strain-mediated electric-field creation of skyrmions in polycrystalline [Pt/Co/Ta] × 5 multilayers, based on the strain-mediated electric-field modulation of effective PMA (K_eff) and interfacial DMI strength (D). Since D decreases sizably when sample was switched from +0 kV/cm (or −0 kV/cm) to −4 kV/cm, which by itself would impede the stabilization of skyrmions, the observed skyrmion creation/stabilization must be attributed to the modulation of K_eff. Here, $$ K_{\mathrm{eff}}\approx K_U-0.5{\mu }_0{\mathbf{M}}_S^2+{\mathbf{M}}_S{\mathbf{B}}_{\mathrm{bias}}+B_1\left(1+\frac{2c_{12}^m}{c_{11}^m}\right){\epsilon }_{\mathrm{in-plane}}$$ (see ref. ²⁹), where K_U is the uniaxial magnetic anisotropy, μ₀ is vacuum permeability, $$ B_1=-1.5{\lambda }_S\left(c_{11}^m-c_{12}^m\right)$$ is the magnetoelastic coupling coefficient of the magnetic multilayer, $$ c_{11}^m$$ and $$ c_{12}^m$$ are its elastic stiffness coefficient, and ε_in-plane is the isotropic in-plane strain the magnetic multilayer experiences when the underlying PMN-PT substrate deforms under applied electric fields. We use MuMax³ (version: 3.10β, which incorporates the module of magnetoelastic coupling) for simulating the evolution of local magnetization in polycrystalline Pt/Co/Ta multilayers under such isotropic in-plane strains transmitted from PMN-PT (see “Methods”).

Figure 3e–h shows the simulated magnetization states at equilibrium under different in-plane biaxial compressive strains $$ {\epsilon }_{\mathrm{in-plane}}={\epsilon }_{\left[110\right]}={\epsilon }_{\left[-110\right]}$$ obtained from Supplementary Fig. 6, with E = +0, −4, −0, and +4 kV/cm, respectively. The average interfacial DMI strength under these electric fields, D = 0.772, 0.585, 0.685, and 0.727 mJ/m², respectively, are directly extracted from Fig. 2b. Overall, the simulations (Fig. 3e–h) agree well with results from MFM imaging (Fig. 3a–d). Figure 3e shows the equilibrium magnetization state under zero strain and D = 0.772 mJ/m² (+0 kV/cm), relaxed from a uniform upward (+z) magnetization. When applying $$ {\epsilon }_{\mathrm{in-plane}}=-0.189\%$$ and D = 0.585 mJ/m² (−4 kV/cm), multiple skyrmions appear (Fig. 3f). As mentioned above, since the reduction of D does not favor the formation of skyrmions, the observed skyrmion creation is attributed to the reduction of the K_eff by the compressive in-plane compressive strains (note that the Pt/Co/Ta multilayer has a negative magnetostriction²⁴). Being consistent with experiments, the created skyrmions can be retained (Fig. 3g) after turning off the electric field (−0 kV/cm), at which $$ {\epsilon }_{\mathrm{in-plane}}=-0.034\%$$ and D = 0.685 mJ/m². However, neither this small remanent strain nor the slightly enhanced D are the key factors for retaining the skyrmions, because the simulations shows that created skyrmions can be retained even when ε_in-plane is set to 0 and D remains to be 0.585 mJ/m² (Supplementary Fig. 14). This control study indicates that the spatially heterogenous magnetic parameters in a polycrystalline system are the key to the retainment of the created skyrmions. For example, the spatial variation of uniaxial local magnetic anisotropy mimics the presence of local pinning sites. Simulations results studying the influence of such parametric heterogeneity, average grain size, size of computation cell, and room-temperature (298 K) thermal fluctuation on skyrmion creation are provided in Supplementary Figs. 15—19. Finally, when applying $$ {\epsilon }_{\mathrm{in-plane}}=-0.010\%$$ and D = 0.727 mJ/m² (+4 kV/cm), the overall morphology of the skyrmion state (Fig. 3h) remains largely the same compared to that in Fig. 3g, partly because the reduced ε_in-plane and the slightly enhanced D act against and in favor of stabilizing the skyrmion state, respectively.

For demonstrating that strain-mediated control of effective PMA plays a dominant role in such skyrmion creation process, Fig. 3i–k shows the evolution of ε_in-plane, magnetoelastic energy density f_mel and intrinsic energy density f_intrin, respectively. The evolution of the constituent energy terms (f_anis, f_exch, f_stray, f_Zeeman, f_DMI, see their expressions in “Methods”) of f_intrin are shown in Supplementary Fig. 20. At t₀ (E = +0 kV/cm), the system has an almost uniform perpendicular magnetization (Fig. 3e). At t₁ (E = −4 kV/m), f_mel increases to 11.85 kJ/m³ associated with the sudden rise of strain to −0.189%. To minimize the f_mel, perpendicular magnetization would then start to rotate toward the film plane (the magnetic easy plane). Indeed, f_mel drops to 1.19 kJ/m³ at equilibrium. During this process, the interfacial DMI energy f_DMI also needs to be minimized (Supplementary Fig. 20), which favors the formation of a Néel domain wall (inset in Fig. 3k). It is the simultaneous reduction of f_mel and f_DMI that provides the thermodynamic driving force for the creation of skyrmions (Fig. 3f). Furthermore, the energy cost is Δf_intrin = 6.74 kJ/m³, while energy reduction is Δf_mel = −10.66 kJ/m³, yielding Δf_tot = Δf_intrin + Δf_mel = −3.92 kJ/m³ < 0. Therefore, the skyrmion creation is favored.

Upon removing the applied large strain at t₂ (E = −0 kV/cm), f_mel increases by 0.63 kJ/m³ because of the increase in the magnitude of perpendicular magnetization (Fig. 3g) and f_intrin decreases by 6.03 kJ/m³. Moreover, f_intrin increases by 0.71 kJ/m³ compared with the initial state (t₀), suggesting that the system in which skyrmions are maintained (Fig. 3g) is metastable. Similarly, further reducing the compressive strain from −0.034 to −0.01% and increasing D slightly from 0.685 to 0.727 mJ/m² at t₃ (E = +4 kV/cm), the corresponding f_mel and f_DMI both decreases, but the former shows a more sizable change (see Supplementary Fig. 20), indicating a stronger strain effect. Since the reduction of an in-plane compressive strain (hence larger K_eff) reduces the stability of skyrmions, some smaller skyrmions annihilate. From the magnetization vector distributions in the local region of the multilayer (insets in the bottom panel), it can be seen that most of the magnetization vectors near the domain wall are perpendicular to the wall, directing from inside of the skyrmion to the outside. This indicates a DMI-stabilized left-chirality Néel wall, which is further confirmed by our asymmetric domain expansion study using L-TEM (Supplementary Fig. 21).

Deformation and annihilation of isolated skyrmions

We also investigated deformation and annihilation of isolated skyrmions using MFM, with a high-spatial resolution via a smaller scanning size (2 μm), which is likely unperceived in Fig. 3 due to the larger scanning size (5 μm). The sample was in the remanent polarization state after positive electric poling (E = +0 kV/cm) with B_bias = 55 mT. Figure 4 demonstrates the anisotropic skyrmion deformation induced by in situ electric fields. The anisotropic deformation factor f is defined as f = 1 − r₁/r₂, with r₁ and r₂ the minor and major axes of the elliptic skyrmion, respectively. In the initial state (E = +0 kV/cm), the skyrmion is circular (Fig. 4a) with the line-cut profiles shown in Fig. 4g. At −4 kV/cm, the skyrmion exhibits evident anisotropic deformation (Fig. 4b) along the diagonal direction of magnetic multilayer, i.e., the [110] direction of PMN-PT, with f ~ 48% (r₁ = 180 ± 12 nm, r₂ = 346 ± 22 nm) deduced from Fig. 4h. At −0 kV/cm (after removing −4 kV/cm), the skyrmion becomes circular again with a size similar to its initial state. A subsequent application of +4 kV/cm (Fig. 4d) followed by its removal (+0 kV/cm, Fig. 4e), there is no appreciable deformation. However, the deformation occurs again upon applying −4 kV/cm (Fig. 4f), suggesting that the deformation is repeatable and appears only under −4 kV/cm. This is consistent with the earlier-discussed observation that the creation of multiple skyrmions also occurs only under −4 kV/cm, where both the strain and interfacial DMI show the most significant change (Fig. 2a, b). Note that the phase contrast in Fig. 4d is somewhat pale, which may be attributed to the influence of magnetic tip³², but this does not contradict the main finding.

Fig. 4

Skyrmion deformation.

Isolated skyrmion morphology at E = +0 kV/cm (a), −4 kV/cm (b), −0 kV/cm (c), +4 kV/cm (d), +0 kV/cm (e), and −4 kV/cm (f) with B_bias = 55 mT. g, h Skyrmion line-cut profiles along the major and minor axes in a and b with Gaussian fittings (solid lines), respectively. The insets in a–f show the simulation results of strain-mediated deformation of one single skyrmion with $$ {\epsilon }_{\left[-110\right]}\mathrm{=}{\epsilon }_{\left[110\right]}\mathrm{=}0$$ D  = 0.772 mJ/m² for +0 kV/cm (a), $$ {\epsilon }_{\left[-110\right]}=-0.169\%$$, $$ {\epsilon }_{\left[110\right]}=0$$, D  = 0.585 mJ/m² for −4 kV/cm (b), $$ {\epsilon }_{\left[-110\right]}={\epsilon }_{\left[110\right]}=-0.023\%$$ D = 0.685 mJ/m² for −0 kV/cm (c), and $$ {\epsilon }_{\left[-110\right]}={\epsilon }_{\left[110\right]}=-0.012\%$$ D = 0.727 mJ/m² for +4 kV/cm (d). The scale bar is 100 nm.

The experimental work on skyrmion deformation is rather limited⁴³. A deformation factor of 30% was reported for B20 compound FeGe under a uniaxial strain of 0.3% (see ref. ⁴³). This exotic emergent elasticity of skyrmions in FeGe was attributed to either the anisotropic magnetoelastic energy^44,45 or the anisotropic DMI strength⁴³ induced by the uniaxial strain. In magnetic multilayers with interfacial DMI, there is only one computational report on such emergent elasticity of isolated skyrmion⁴⁶, which was attributed to anisotropic magnetoelastic energy.

In our sample, although the average in-plane strains from a [001]-poled (001)PMN-PT crystal are expected to be isotropic based on its crystal symmetry³⁴, the local in-plane strain at the PMN-PT surface, arising from the 109° ferroelectric domain switching^25,36, should be anisotropic. It is such local strain anisotropy that leads to the deformation of a single local skyrmion, which is much smaller than the micron-sized ferroelectric domain at the PMN-PT surface. Specifically, such local anisotropic strains, which are along the [110] and [−110] axes of the PMN-PT, would lead to anisotropic Néel wall energy $$ {\sigma }_{\mathrm{w}}=4\sqrt{A_{\mathrm{ex}}{K}_{\mathrm{eff}}}-\pi ∣D∣$$ (refs. ^8,20) along these two axes by modulating either the local magnetoelastic anisotropy (which contributes to K_eff) or the local DMI strength D, or both. Here, A_ex is the exchange stiffness. Regardless of the details, the skyrmion radius is always smaller along the axis with a lower σ_w (see ref. ⁴⁶). This is analogous to Wulff construction: lower surface energy (wall energy) of a crystal (skyrmion) yields shorted vector length (skyrmion radius) at thermodynamic equilibrium.

With these in mind, we simulate the deformation of one single skyrmion by applying biaxial in-plane anisotropic strains. The anisotropic magnetoelastic energy f_mel was automatically considered (see expression of f_mel in “Methods”), but the local DMI strength D was assumed to be isotropic for simplicity. The simulation results, which agree well with experimental observations, are shown in the insets of Fig. 4.

We also observed repeatable annihilation and creation of isolated skyrmion (Fig. 5). The initial state (Fig. 5a) shows the skyrmion with irregular shape implying presence of local inhomogeneity. At E = −4 kV/cm, there is no phase contrast (Fig. 5b), indicating the annihilation of skyrmion. Upon the removal of this negative poling field, the skyrmion reappears at E = −0 kV/cm (Fig. 5c). Then, applying +4 kV/cm, the skyrmion shrinks in size (Fig. 5d). After removing this positive poling field (now at E = +0 kV/cm), the skyrmion returns to its initial state (Fig. 5e). Applying −4 kV/cm again, the skyrmion annihilates again (Fig. 5f).

Fig. 5

Skyrmion creation/annihilation.

Isolated skyrmion morphology at E = +0 kV/cm (a), −4 kV/cm (b), −0 kV/cm (c), +4 kV/cm (d), +0 kV/cm (e), and −4 kV/cm (f) with B_bias = 55 mT. The insets in a–f show the simulation results of strain-mediated annihilation and reappearance of one single skyrmion with $$ {\epsilon }_{\left[-110\right]}={\epsilon }_{\left[110\right]}=0$$, D = 0.772 mJ/m² for E = +0 kV/cm (a, e); $$ {\epsilon }_{\left[-110\right]}={\epsilon }_{\left[110\right]}=-0.0425\%$$, D = 0.585 mJ/m² for E = −4 kV/cm (b, f); $$ {\epsilon }_{\left[-110\right]}={\epsilon }_{\left[110\right]}=-0.0415\%$$, D = 0.685 mJ/m² for E = −0 kV/cm (c); and $$ {\epsilon }_{\left[-110\right]}={\epsilon }_{\left[110\right]}=-0.012\%$$, D = 0.727 mJ/m² for E = +4 kV/cm (d). A 20-nm diameter pinning site with ~5% lower perpendicular anisotropy was specified, as indicated by the dashed circle in the center. The scale bar is 100 nm.

The insets in Fig. 5 show the simulation results for the processes discussed above. Notably, for simulating the reappearance of the skyrmion from exactly the same location where it annihilates, we found that it is necessary to introduce a small pinning site (20 nm in diameter herein) with moderately (~5%) lower perpendicular anisotropy. The critical role of lower-anisotropy pinning site in achieving localized skyrmion deletion/re-creation has been demonstrated in systems with Rashba-type DMI arising from the FM/dielectric interface¹⁹. As shown in the inset of Fig. 5b, this lower-anisotropy pinning site leads to the stabilization of a tiny (diameter ~10 nm) skyrmion, which is too small to discern with our MFM (with a spatial resolution of ~10 nm as well), yet can function as the nucleus for the reappearance (growth) of skyrmion at exactly the same location after the removal of electric field (the inset of Fig. 5c). Our simulation results suggest the possibility of harnessing pinning sites for realizing a spatially precise control of skyrmions in continuous magnetic layer, which could be useful for device applications.

It is worth noting that there is an existing report on tuning the density and stability of the skyrmions in BaTiO₃/SrRuO₃ bilayer heterostructures by ferroelectric polarization switching⁴⁷. In that work, voltage-controlled topological Hall effect (THE) measurement to provide an indirect evidence for skyrmion manipulation at low temperatures (<100 K), but a direct observation via in situ imaging is missing. Moreover, the origin of the so-called THE signal in SrRuO₃ is still an open question^48,49. After completing our manuscript, we noticed a very recent report on creation of micron-sized skyrmion bubbles (3–6 μm) with in-plane radiofrequency voltages by employing surface acoustic waves in multilayers of Pt/Co/Ir grown on LiNbO₃, where it was found that static electric fields hardly affect skyrmions⁵⁰. By contrast, in our work, we show manipulations of nanometer-scale skyrmions by using static and vertical electric fields. The use of static and vertical electric field is generally more compatible with existing manufacturing platform of spintronics, and is more energy efficient.

Sample preparation

PMN-PT FE single-crystal substrates are (001)-cut and one-side polished with a size of 5 × 3 × 0.5 mm³. Ta(4.7 nm)/[Pt(4 nm)/Co(1.6 nm)/Ta(1.9 nm)] × 5 multilayers were grown by magnetic sputtering at room temperature with a base pressure of 1 × 10⁻⁶ Pa without a magnetic field. The first layer, Ta, was used to flat the surface of PMN-PT and promote the interface quality of Pt/Co. The top Ta layer is a capping layer preventing oxidation of Co and also as the top electrode. Au layer with a thickness of 300 nm was deposited on the bottom of the PMN-PT as the electrode.

Magnetic property measurement system

Out-of-plane and in-plane magnetic hysteresis loops were measured with a magnetic property measurement system (MPMS 7 T; Quantum Design).

Magnetic force microscopy

MFM was performed in ambient conditions with an Infinity Asylum Research AFM. The stray field of magnetic tip will inevitably interact with magnetic domain and distort them³². To reduce the influence of magnetic tip on domain, an optimized magnetic tip (SSS-MFMR)⁵¹ has been used in our work, which has a radii of 30 nm with a low coercivity ~ 12 mT and ultralow magnetic moment ~ 8 × 10⁴ A/m. MFM images were taken in tapping/lift mode with a lift height of 30 nm.

The safety voltage for MFM is 200 V. Therefore, we apply ±4 kV/cm to the sample with a thickness of 0.5 mm.

Lorentz transmission electron microscopy

L-TEM imaging was carried out on Lorentz TEM (JEOL 2100F) under perpendicular magnetic fields by gradually increasing the objective lens current.

Scanning transmission electron microscopy

STEM and EDS mapping were performed on an ARM-200F (JEOL) operated at an acceleration voltage of 200 kV and equipped with double-spherical aberration (Cs) correctors. The TEM sample was prepared by using a focused ion beam with the cross-sectional lamella thinned down to 100 nm.

X-ray diffraction

XRD was operated on a Rigaku D/max-RB X-ray diffractometer with Cu Kα radiation.

Brillouin light spectroscopy

The principle of BLS is that⁵² the wave-vector direction of the magnetostatic surface spin wave in FM layer has a spatial chirality depending on the cross products of the normal direction and the direction of the magnetization, thus a pair of symmetric Stokes and anti-Stokes peaks are present in the BLS spectra corresponding to a fixed wave vector but opposite direction. However, the interfacial DMI favors certain spatial chirality same with only one of the wave vector, which induces an asymmetric modification of the frequency dispersion relation of Stokes and anti-Stokes peaks. The wave-vector dependence of the frequency difference between the Stokes and anti-Stokes peaks, Δf, has a linear relation with the in-plane wave vector k, and the slope is determined by D, see Eq. (1).

BLS was performed by using a single-mode solid-state laser with a wavelength of 532 nm and a power of 30 mW. The 180°-backscattered light was analyzed using a JRS Sandercock-type multipass tandem Fabry–Pérot interferometer. An external magnetic field of 500 mT was applied parallel to the sample surface plane and perpendicular to the in-plane wave vector (DE mode configuration).

FMR measurement

FMR measurements were performed by using the JEA-FS200 ESR spectrometer of JEOL Co., equipped with an X-band microwave generator with a frequency of 9.1 GHz and a cylindrical microwave resonant cavity. The sample was mounted to a rotating stage with the magnetic field direction unchanged.

Micromagnetic simulation

We use MuMax³ (see ref. ⁵³) (version: 3.10β, which incorporates the module of magnetoelastic coupling) for simulating the evolution of local magnetization in polycrystalline [Pt(4 nm)/Co(1.6 nm)/Ta(1.9 nm)] × 5 multilayers under voltage-induced strains transmitted from the underlying piezoelectric PMN-PT substrate. In the simulations, the total free energy density f_tot = f_anis + f_exch + f_stray + f_Zeeman + f_mel + f_DMI is contributed by the following energy densities: the magnetocrystalline uniaxial anisotropy energy density $$ f_{\mathrm{anis}}=-K_{\mathrm{u}}{\left(\mathbf{u}\cdot \mathbf{m}\right)}^2$$, where u is the unit vector of the uniaxial anisotropy and m is the normalized local magnetization, K_u is the uniaxial anisotropy constant; the magnetic exchange energy density $$ f_{\mathrm{exchange}}=A_{\mathrm{ex}}{\left(\nabla \mathbf{m}\right)}^2$$, where A_ex is magnetic exchange coupling coefficient; the magnetostatic stray field energy density $$ f_{\mathrm{stray}}=-\frac{1}{2}{M}_{\mathrm{S}}\mathbf{m}\cdot {\mathbf{B}}_{\mathrm{stray}}$$, where M_S is the saturation magnetization and B_stray is the magnetostatic stray field; the Zeeman energy density associated with the externally applied magnetic field B_ext, i.e., $$ f_{\mathrm{Zeeman}}=-M_{\mathrm{S}}\mathbf{m}\cdot {\mathbf{B}}_{\mathrm{ext}}$$; the magnetoelastic energy density $$ f_{\mathrm{mel}}=B_{\mathrm{1}}\left({\epsilon }_{\mathrm{xx}}{m}_{\mathrm{x}}^2+{\epsilon }_{\mathrm{yy}}{m}_{\mathrm{y}}^2+{\epsilon }_{\mathrm{zz}}{m}_{\mathrm{z}}^2\right)+2B_{\mathrm{2}}\left(m_{\mathrm{x}}{m}_{\mathrm{y}}{\epsilon }_{\mathrm{xy}}+m_{\mathrm{x}}{m}_{\mathrm{z}}{\epsilon }_{\mathrm{xz}}+m_{\mathrm{y}}{m}_{\mathrm{z}}{\epsilon }_{\mathrm{yz}}\right)$$, where B₁ and B₂ are magnetoelastic coupling coefficients, but the B₂-related term can be omitted in our case because average shear strains are set to be zero $$ \left({\epsilon }_{\mathrm{xy}}={\epsilon }_{\mathrm{xz}}={\epsilon }_{\mathrm{yz}}=0\right)$$; the interfacial DMI energy density $$ f_{\mathrm{DMI}}=D\left[m_{\mathrm{Z}}\left(\nabla \cdot \mathbf{m}\right)-\left(\mathbf{m}\cdot \nabla \right)m_{\mathrm{Z}}\right]$$, where D is the interfacial DMI strength. The variation of f_tot with respect to the change in m gives rise to the effective magnetic field $$ {\mathbf{B}}_{\mathrm{eff}}=-\frac{1}{M_{\mathrm{S}}}\frac{\delta F_{\mathrm{tot}}}{\delta \mathbf{m}}$$, which drives the motion of the local magnetization m in a way described by the Landau–Lifshitz–Gilbert equation. The influence of thermal fluctuation at room temperature (298 K) on the switching of m can be considered by adding a thermal fluctuation field B_therm to the total effective magnetic field B_eff, that is, $$ {\mathbf{B}}_{\mathrm{eff}}=-\frac{1}{M_{\mathrm{S}}}\frac{\delta F_{\mathrm{tot}}}{\delta \mathbf{m}}+{\mathbf{B}}_{\mathrm{therm}}$$. Here, $$ {\mathbf{B}}_{\mathrm{therm}}=\eta \sqrt{\frac{2\alpha k_{\mathrm{B}}T}{M_{\mathrm{S}}\gamma \mathrm{\Delta }V\mathrm{\Delta }t}}$$, where α is the damping parameter, k_B is Boltzmann constant, T is temperature, γ is the gyromagnetic ratio, ΔV is the volume of simulation cell, and Δt is time step in real unit, $$ \mathbf{\eta }=\left({\eta }_{\mathrm{x}},{\eta }_{\mathrm{y}},{\eta }_{\mathrm{z}}\right)$$, where η_x, η_y, η_z are random numbers that have a standard normal distribution and a mean of zero, and are uncorrelated in both space and time. It is found that adding B_therm (with T = 298 K) yields overall similar behaviors of strain-mediated skyrmion creation to the case where B_therm = 0 (by setting T = 0 K), as shown in Supplementary Fig. 19. The simulations shown in the main paper are performed at B_therm = 0 for investigating and analyzing the energetics of the skyrmion switching more clearly.

The simulations simultaneously consider strain-induced changes in both the magnetic anisotropy and the interfacial DMI strength D. Unless stated otherwise, the former was considered by using average strains measured using strain gauge under different electric fields (Fig. 2a), which are related to the magnetoelastic energy density f_mel. The latter was considered by using the D values measured using BLS under different electric fields (Fig. 2b), which are relate to the interfacial DMI energy density f_DMI. The influence of strain on the magnetic anisotropy can be interpreted in a way that the strain is modulating the effective PMA of the multilayer as discussed in the main text.

The multilayer is described by 3D discretized computational cells of $$ n_{{\mathrm{x}}^{\prime }}\mathrm{\Delta }{x}^{\prime }\times n_{{\mathrm{y}}^{\prime }}\mathrm{\Delta }{y}^{\prime }\times 9\mathrm{\Delta }{z}^{\prime }$$, where n_x′ and n_y′ are numbers of computational cells, and Δx′, Δy′, and Δz′ are cell sizes. Along the thickness direction (z′ axis), the bottom five layers and top two layers of grids with Δz′ = 0.8 nm are set as nonmagnetic materials to describe the Pt and Ta layer, while the middle two layers describe the magnetic Co layer. In-plane periodic boundary conditions are applied by adding one repetition of the simulation system along both x′ and y′ axes. Five repetitions of the simulation system are set along z′ axis to describe five repetitions of the Pt/Co/Ta tri-layer. The coordinate systems used in the experiment and simulations are different, shown in Supplementary Fig. 1b. Specifically, we only consider the strain along the diagonal direction of PMN-PT(001), and $$ {\epsilon }_{{\mathrm{x}}^{\prime }{\mathrm{x}}^{\prime }}^m={\epsilon }_{\left[110\right]}^{\mathrm{PMN}-\mathrm{PT}}={\epsilon }_{\left[110\right]}$$, $$ {\epsilon }_{{\mathrm{y}}^{\prime }{\mathrm{y}}^{\prime }}^m={\epsilon }_{\left[-110\right]}^{\mathrm{PMN}-\mathrm{PT}}={\epsilon }_{\left[-110\right]}$$ as discussed later. We use $$ {\epsilon }_{\mathrm{in-plane}}={\epsilon }_{\left[110\right]}={\epsilon }_{\left[-110\right]}$$ in the main text.

For simulating strain-mediated electric-field control of multiple skyrmions, we first relax an initially uniform [001] magnetization under a bias external magnetic field B_bias. The obtained equilibrium magnetization distribution is used as the initial state. We then apply a spatially uniform strain to the entire magnetic multilayer. A complete transmission of the voltage-induced strain from the PMN-PT surface to the magnetic multilayer is considered²⁶, yielding $$ {\epsilon }_{{\mathrm{x}}^{\prime }{\mathrm{x}}^{\prime }}^m={\epsilon }_{\left[110\right]}^{\mathrm{PMN}-\mathrm{PT}}$$, $$ {\epsilon }_{{\mathrm{y}}^{\prime }{\mathrm{y}}^{\prime }}^m={\epsilon }_{\left[-110\right]}^{\mathrm{PMN}-\mathrm{PT}}$$. The out-of-plane strain of the magnetic layer $$ {\epsilon }_{{\mathrm{z}}^{\prime }{\mathrm{z}}^{\prime }}^m$$ is calculated by $$ {\epsilon }_{{\mathrm{z}}^{\prime }{\mathrm{z}}^{\prime }}^m=-\frac{c_{12}^m\left({\epsilon }_{{\mathrm{x}}^{\prime }{\mathrm{x}}^{\prime }}^m+{\epsilon }_{{\mathrm{y}}^{\prime }{\mathrm{y}}^{\prime }}^m\right)}{c_{11}^m}$$ according to the thin-film boundary condition.

For the simulation of skyrmions creation in Fig. 3, we describe the system as a polycrystal considering the heterogeneity in a polycrystalline sample, with a mean grain size of 20 nm, $$ \mathrm{\Delta }{x}^{\prime }=\mathrm{\Delta }{y}^{\prime }=3\mathrm{nm}$$, and $$ n_{{\mathrm{x}}^{\prime }}=n_{{\mathrm{y}}^{\prime }}=1000$$ and B_bias = 60 mT. We define grain-like regions with region-specific uniaxial magnetic anisotropy constant (K_U), uniaxial magnetic anisotropy direction, and interfacial DMI strength (D) using Voronoi tessellation. Specifically, we generate 256 different K_U obeying a uniform distribution around K_U0 = 6.5 × 10⁵ J/m³ (based on our own measurement, see Supplementary Fig. 2) within the range of (0.9–1.1K_U0), 256 different θ (the angle between the axis of uniaxial magnetic anisotropy and the z′-axis) obeying a uniform distribution around 0° within the range of (−20°–20°), 256 different φ (the angle between +x′-axis and the projection of axis of uniaxial magnetic anisotropy on the x′y′ plane) obeying uniform distribution within the range of (0°–360°), and 256 different D obeying a uniform distribution around D₀ within the range of (0.9–1.1D₀). For each parameter set (K_U, θ, φ, D), the 256 different numbers are randomly assigned to the grain-like regions. Note that the D is the averaged interfacial DMI strength of the entire sample, which is taken from our own experimental measurement and dependent on the applied electric fields (see inset of Fig. 2b). Furthermore, for a (001)-oriented PMN-PT single crystal, the average in-plane strain induced by electric field is expected to be biaxial isotropic⁵⁴, i.e., $$ {\epsilon }_{\left[110\right]}^{\mathrm{PMN}-\mathrm{PT}}={\epsilon }_{\left[-110\right]}^{\mathrm{PMN}-\mathrm{PT}}={\epsilon }_0^{\mathrm{PMN}-\mathrm{PT}}$$. Assuming volume conservation, $$ {\epsilon }_0^{\mathrm{PMN}-\mathrm{PT}}$$ can be calculated from the experimentally measured $$ {\epsilon }_{\left[001\right]}^{\mathrm{PMN}-\mathrm{PT}}$$, which is the average out-of-plane strain induced by electric field, via $$ {\epsilon }_0^{\mathrm{PMN}-\mathrm{PT}}={\nu }_{\mathrm{PMN}-\mathrm{PT}}\cdot {\epsilon }_{\left[001\right]}^{\mathrm{PMN}-\mathrm{PT}}$$, with Poisson’s ratio $$ {\nu }_{\mathrm{PMN}-\mathrm{PT}}=-0.5$$ (see ref. ³⁴), shown in Supplementary Fig. 6.

However, given the local strain variations on the PMN-PT surface caused by spatially heterogenous ferroelectric domain switching^25,36, we introduce anisotropic in-plane strain to simulate the deformation and annihilation of one single skyrmion as discussed in the Figs. 4 and 5 and Supplementary Fig. 22. In these cases, we use $$ \mathrm{\Delta }{x}^{\prime }=\mathrm{\Delta }{y}^{\prime }=0.5\mathrm{nm}$$ and $$ n_{x^{\prime }}=n_{y^{\prime }}=600$$ as well as spatially homogenous magnetic parameters (K_U = 6.5 × 10⁵ J/m³; D varies with electric fields) for simplicity; the bias magnetic field B_bias = 30 mT. Under this setup, a circular-shaped 68-nm diameter skyrmion can be stabilized under zero strain. Anisotropic in-plane strains are then applied.

The rest of the magnetic parameters are all set to be spatially uniform in both the polycrystal and single-crystal (i.e., spatially uniform K_u and D) model. These parameters include: saturation magnetization M_S = 9.48 × 10⁵ A/m (from our own experimental measurement), magnetic exchange coupling coefficient A_ex = 1.0 × 10⁻¹¹ J/m (see ref. ¹⁴), magnetostriction constant λ_S = −180 ppm (which leads to best fit of the experimentally measured skyrmion morphology), and elastic stiffness coefficients $$ c_{11}^m=290.86\mathrm{GPa}$$ and $$ c_{12}^m=160\mathrm{GPa}$$, which are calculated from thickness-weighted average of the Young’s modulus and Poisson’s ratio of the Pt, Co, and Ta layers⁵⁵. Young’s modulus are $$ Y_{\mathrm{Pt}}=160\mathrm{GPa}$$, $$ Y_{\mathrm{Co}}=210\mathrm{GPa}$$, and $$ Y_{\mathrm{Ta}}=186\mathrm{GPa}$$, and Poisson’s ratio are $$ {\nu }_{\mathrm{Pt}}=0.38$$, $$ {\nu }_{\mathrm{Co}}=0.31$$, and $$ {\nu }_{\mathrm{Ta}}=0.34$$ (https://periodictable.com/Properties/A/YoungModulus.html).

Note added to proof

After submission of this work and during the period of preparation for the response to the reviewers’ comments, an article reporting study on electric-field control of skyrmions in multiferroic ferromagnetic/ferroelectric nanostructures⁵⁶ has appeared. We focus on skyrmions in continuous thin films, while they focus on nanopatterned confined skyrmions and the behaviors of the skyrmions in these two cases are different.