Controlling the helicity of magnetic skyrmions by electrical field in frustrated magnets

The skyrmions generated by frustration in centrosymmetric structures host extra internal degrees of freedom: vorticity and helicity, resulting in distinctive properties and potential functionality, which are not shared by the skyrmions stemming from the Dzyaloshinskii-Moriya interaction in noncentrosymmetric structures. The present work indicates that the magnetism-driven electric polarization carried by skyrmions provides a direct handle for tuning helicity. Especially for the in-plane magnetized skyrmions, the helicity can be continuously rotated and exactly picked by applying an external electric field for both skyrmions and antiskyrmions. The in-plane uniaxial anisotropy is beneficial to this manipulation.


Introduction
In past century, the topological magnetic textures stabilized by complicated interactions were predicted and investigated extensively [1][2][3][4][5][6]. In particular, topological vortex-like magnetic skyrmions, have attracted a great deal of interest from both the academic and technological fields [7,8]. Since the experimental discovery in 2009 [9], magnetic skyrmions have been observed in many different materials including metals [10][11][12], semiconductors [13][14][15], insulators [16,17], and thin film systems [18,19]. Besides the formation of independent skyrmion excitations in the ferromagnetic background, these skyrmions also crystallize into stable lattice, surviving even down to zero-temperature. Experiments revealed that the motion of skyrmion can be driven and controlled by the ultralow electric currents in the metallic system [20,21]. On the other hand, the skyrmions can magnetically induce electric polarization in insulators or semiconductors, which enables the modulation of skyrmions by an external electric field without losses due to Joule heating [22][23][24][25].
Up to now, magnetic skyrmions are mostly observed in experiments on noncentrosymmetric materials or interfacial symmetry-breaking heterostructures. It has been widely confirmed that these skyrmions are primarily induced by the Dzyaloshinskii-Moriya interaction (DMI), stemming from inversion-symmetry-breaking and the spin-orbit interaction of different kinds. For example, Bloch-type skyrmions can be stabilized in chiral magnets with the Dresselhaus spin-orbit interaction, while Néel-type skyrmions may appear in polar magnets with the Rashba spin-orbit interaction or in magnetic ultrathin films with interfacial DMI [30][31][32][33]. Meanwhile, the magnetic antiskyrmions discovered in Heusler compounds are stabilized by the anisotropic DMI with opposite signs along two orthogonal in-plane directions [34][35][36][37][38]. Although different skyrmions may be observed in different materials, their morphology in one compound is always fixed because the DMI is determined by the structure. On the other hand, various alternative 3 mechanisms are proposed, especially for the skyrmions in centrosymmetric materials, such as the long-ranged magnetic dipolar interactions [17], the four-spin exchange interactions [18], the coupling between itinerant electron spins and localized spins [39,40]. In particular, the magnetic frustration has long served as a relatively simple yet rich source of novel magnetic phases [41][42][43][44][45]. In 2012, it was first theoretically proposed that the skyrmion crystal can be stabilized by the frustrated interactions [46], but the following experimental reports on the skyrmions in frustrated systems were very limited [47,48]. It was exciting that very recently the skyrmion crystal was experimentally discovered in the frustrated centrosymmetric triangular-lattice magnet Gd 2 PdSi 3 [49] and breathing kagomé lattice magnet Gd 3 Ru 4 Al 12 [50]. Due to the frustration source, the skyrmions are typically much smaller than the DMI-driven counterpart, and thus contribute to a giant topological Hall effect. Later, nanometric square skyrmion lattice was reported in a centrosymmetric tetragonal magnet GdRu 2 Si 2 [51].
It has been theoretically predicted that the frustration-induced skyrmions in centrosymmetric structure show distinctive properties [52][53][54][55][56]. One fascinating advantage is the extra internal degrees of freedom-vorticity and helicity. The helicity-dependent current responses unveiled by the dynamics simulation imply that these factors not only characterize the morphology of skyrmion, but also can be exploited to control its motion and furthermore to achieve diverse functionality. On the contrary, for the DMI system, both vorticity and helicity are locked by the sign and direction of the DMI vector. The only way to tune the helicity of the DMI-driven skyrmion is the controllable DMI, i.e. to tune the spin-orbit coupling [57,58]. But the helicity is still locked for a certain composition in these previous reports. Compare to the fixed morphology of the DMI-driven skyrmion, the frustration-induced skyrmion holds the intrinsic degrees of freedom to tune its morphology, and therefore shows the related exotic behaviors not shared by the DMI systems. How to take advantage of these freedom degrees to realize an effective and flexible manipulation remains a crucial physical problem to be resolved. Very recently it was reported that the helicity of dipolarly stabilized skyrmions can be tuned by varying the material parameters and 4 geometry of nanostructure [59,60].

Model and methods
The frustration-driven skyrmion crystal was first predicated as metastable state realized in triangular lattice under applied magnetic field perpendicular to the lattice plane at moderate temperatures [46]. The later investigation testified that an easy-axis anisotropy perpendicular to the lattice plane is sufficient to stabilize the skyrmion crystal with the lowest energy at zero-temperature [53,54,61]. Here the frustrated magnetic model on two-dimensional triangular lattice is considered with the Hamiltonian expressed as where S i represents a classic spin of unit length at the i-th site on xy-plane. The first three terms on the right of Eq. 1 are the exchange energies, where J 1 , J 2 and J 3 represent the exchange interactions between spins on the nearest-neighboring (<i,j>), 5 the second-nearest-neighboring (<<i, m>>) and the third-nearest-neighboring (<<<i, n>>>) sites. Ferromagnetic J 1 =1 is fixed as the energy unit, and all the parameters are simplified with reduced units. It has been reported that either of antiferromagnetic J 2 and J 3 considered is enough to produce a stable skyrmion crystal, and the similar results can be obtained in these two cases [46,61]. So J 2 =0 and J 3 =-0.5 are set unless otherwise noted. The fourth and fifth terms describe the energies of uniaxial anisotropy and magnetic field, where the uniaxial anisotropy with the magnitude k along e k direction and the magnetic field with the magnitude h along e h direction.
Due to the frustration source, the present skyrmion is typically much smaller than the DMI-driven counterpart. The noncollinear spin texture in such a scale of the lattice constant usually provides the condition to produce electric polarization according to the spin current mechanism [62], where an electric dipole p ij can be induced by two neighboring canting spins (S i and S j ) in the form of where e ij denotes the unit vector connecting the two sites of neighboring S i and S j .
The total electric polarization (P) is estimated by the summary over all the bonds, which also can be written as the summary of all the local electric polarization on site Considering the ferroelectric energy with external electric field (E), the total Hamiltonian can be written as The simulation is performed on the triangular lattice of sites N=5184 with periodic boundary conditions. Different lattice sizes are checked to confirm the stability of the main results. The Metropolis algorithm combined with the over-relaxation method is applied to find the stable state with the lowest energy [63,64]. On every parameter point, the system is first evolved from a relatively high temperature to a very low temperature gradually, and then the energy is further minimized to approach the limit of zero temperature. The final result is obtained by 6 comparing independent data sets evolving from different initial states. The spin dynamics at zero temperature is discussed by using the fourth-order Runge-Kutta method to numerically solve Landau-Lifshitz-Gilbert (LLG) equation as below.
where the Gilbert damping coefficient α=0.2 to ensure quick relaxation to the equilibrium state, and the time t is measured in units of ћ/J 1 .
The obtained spin configuration is characterized by the spin structure factor, which is evaluated for three spin components (γ=x, y and z) respectively as follows, The topological character is confirmed by the skyrmion number or topological charge (Q), which is defined as Q quantifies the number of times spin vectors wrapping around a unit sphere as the coordinate (x, y) spans the whole planar space, which can be calculated for a spin lattice in the manner as Ref. [65]. To elucidate the detailed nature, the local topological charge density ρ(r) can be expressed as For an isolated skyrmion in the ferromagnetic background in two-dimensional system with a perpendicular magnetic field, the morphology is determined by vorticity and helicity. If the polar coordinate (r, ϕ) is used with the symmetric center of spin texture as the origin, the spin vector can be expressed as S=(sinθcosφ, sinθsinφ, cosθ). θ only depends on r. φ=nϕ+η, where integer n is the vorticity and η is the helicity. Skyrmion and antiskyrmion are distinguished conventionally by the sign of n, i.e. n=1 for skyrmion and n=-1 for antiskyrmion [7]. The helicity has continuous value from -π to π, where η=-π and π are identical. In particular, when n=+1, η=0 or π is the Néel type, η=-π/2 or π/2 is the Bloch type. All the structures of antiskyrmions 7 (n=-1) are equivalent on rotation in the xy-plane [7].

Results and discussion
For the frustrated triangular model with J 1 -J 2 or J 1 -J 3 , rich phase diagrams have been reported [53,54,61,66]. Under the consideration of uniaxial anisotropy and magnetic field both perpendicular to the lattice plane (e k and e h along +z direction), the perpendicularly magnetized (PM) skyrmion crystal emerges at intermediate h and 8 above a small k, and remains down to zero-temperature. Consistent with the previous reports, the spin configuration of skyrmion crystal phase in the present simulation shows the typical spin structure factor of a triple-q magnetic orderings [46], as presented in Figs. 1(a) and 1(b). The topological character is confirmed by the nonzero value of topological charge, and the corresponding topological charge density values, while they give values near zero for Bloch type of η=±π/2. In the case of antiskyrmion (n=-1) as plotted in Fig. 1(j), all the components of p i cancel out each other, i.e. P=0. The P dependence on η for skyrmion of n=1 is plotted in Fig. 2(a).
The magnitude of P (P) varies with η continuously, following P=P m cos(η) where P m is the amplitude. Thus P could be an characteristic index of η, which may be measured more easily. Néel skyrmion crystal with η=0 o r π always presents the maximum value of P (P m ) with opposite orientations. For the Néel skyrmion crystal, increasing h or k reduces P m owing to the suppression of the noncollinear spin parts.
In addition, Ns mainly depends on the frustration tuned by the exchange interaction J 3 .
When the frustration is enhanced by raising |J 3 | with h and k fixed, Ns increases and P m also, but P m per skyrmion (P m /Ns) decreases as plotted in Figs. 2(b-d). Since the states with different η are energetically degenerate, it is easy to align skyrmions with different η to Néel type of η=0 or π by applying a low electric field.
Although the spatially homogeneous electric field cannot move the skyrmion [67], P can be flipped by E without energy loss. When an AC electrical field 10 E z =E m cos(2πt/600) with E m =0.01 and 0.02 is applied along z-axis (Fig. 2(e)), P flips between P m and -P m , and simultaneously η flips between 0 and π, as illustrated in Figs. 2(f) and 2(g). If E m is stronger, the flips occur more promptly. The independent PM skyrmion also exists in the ferromagnetic state near the skyrmion crystal phase in this system. The simulation on the metastable isolated PM skyrmion shows the very similar magnetoelectric behavior. For instance, an isolated skyrmion presents nearly the same E-driven flips of P and η as plotted in Figs. 2(f) and 2(g), where the same parameters are adopted for crystal and isolated states for a better comparison. It is noteworthy that the conventional skyrmions in most investigations are the 11 PM skyrmions with the spin at the core antiparallel to the spins at the perimeter both perpendicular to the lattice plane. As discussed above, the tunability on helicity by electrical field is limited for the PM skyrmions. It is hard to exactly pick one certain helicity for skyrmion, and it is nearly impossible to control the helicity for antiskyrmion without P. It is worth noting that a novel in-plane magnetized (IM) skyrmion, where the spin at the core still antiparallel to the spins at the perimeter but both within the lattice plane, was put forward in 2019 [68,69]. This unconventional skyrmion was also known as one kind of bimeron [70], and the similar spin texture had been observed in experiments [71]. The IM skyrmion can be obtained from the PM skyrmion by a 90° rotation around the y-axis [68], which can be realized by rotating magnetic field. The IM skyrmion is topologically equivalent to the corresponding PM one with the same Q, because they can be associated with a smooth deformation of the spin texture.
Therefore, the topological properties are well retained in this IM version of skyrmion, and a pure topological Hall effect was predicted [68]. The two internal degrees of freedom-vorticity and helicity are also retained, and thus they can be estimated as the PM case. where C: skyrmion crystal and S: isolated skyrmion. Helicity η as functions of t from varied initial η to reach η=0.3π by applying E of (e) E x =0.01cos(-0.3π) and E y =0.01sin(-0.3π) for IM skyrmion crystal, (f) E x =0.01cos(0.3π+π) and E y =0.01sin(0.3π+π) for IM isolated antiskyrmion.
The IM skyrmion can be regarded as a pair of tight-binding meron and antimeron as plotted in Fig. 3(e). The corresponding p i map demonstrates that the IM skyrmion 13 produces the local electric polarization in an asymmetric way. It is attractive that in this case with e k and e h both along x-axis, P shows nonzero value for all the (anti)skyrmions of different η. Furthermore, both orientation and magnitude of P depend on vorticity and helicity. Figs. 4(a) and 4(b) plot the components of P per skyrmion as functions of η, presenting different behaviors for skyrmion and antiskyrmion crystals. P x =P z =P a cos(-η)=P a cos(η) and P y =P a sin(-η)=-P a sin(η) in the skyrmion case, while P x =P a cos(η+π)=-P a cos(η), P y =P a sin(η+π)=-P a sin(η) and It should be mentioned that the IM skyrmion crystal as stable state down to zero temperature can also be realized with easy-plane anisotropy (k<0) along z-axis and in-plane magnetic field. When k is non-negligible, different from the continuous η , , will be preferred in the metastable IM skyrmion crystal with easy-axis anisotropy (k>0) along z-axis and in-plane magnetic field, as plotted in Fig. 5(b). It is means that the degeneracy of helicity is resolved partially by the nonzero k, and six states with lower energy remain. The corresponding P still roughly follows the rules mentioned above, except that P x and P z do not overlap any 15 more due to k along z. When E with the magnitude E m =0.01 is rotated by ϕ E =0.0005πt in xy-plane, the time-dependence of η demonstrates the sloping steps corresponding to the six η values, respectively for k<0 and k>0 as displayed in Figs. 5(c) and 5(d). The stronger anisotropy induces more obvious steps, implying more robust states with the six η values. In this case, the skyrmions are actually deformed by anisotropy. However, as long as the anisotropy is not so strong, the skyrmion crystal can be well preserved due to the topological protection. On the other hand, the isolated skyrmion may also appear. But it is hard to produce isolated skyrmion in ferromagnetic background in this parameter region. Therefore it is hard to achieve an efficient electric control of η because of the complex background.
It should be mentioned that the IM skyrmions in both crystal and isolated forms may survive as metastable state in the case without anisotropy (k=0). Here the helicity can be tuned continuously and exactly for skyrmion and antiskyrmion crystals respectively, just like the case with in-plane uniaxial anisotropy (Fig. 4), whereas the isolated IM skyrmion can not be controlled freely also due to the complex background.
The better tunability can be achieved for crystal form because the contribution from the background is relatively small.
Although the magnetic dipole interaction prefers the Bloch-type texture, the corresponding energy cost of helicity variation is very small, and it is even much smaller for IM skyrmions than that of PM skyrmions. The electric field applied will overcome this energy barrier to realize the helicity control. (See Supplemental Material Note 1 for details.) The magnetic skyrmion is a natural topological texture in two spacial dimensions. In three dimensional space, the situation will be complicated and interesting for frustrated interlayer interaction, where the characterization and tenability of internal freedom degrees will be explored in the future studies (See Supplemental Material Note 2 for details.)

Conclusions
The PM and IM skyrmions driven by frustration are investigated on a triangular lattice with competing interactions in both crystal and isolated forms. Aiming at 16 insulating systems, the magnetoelectric properties are focused on to explore unique electric controllability on the intrinsic helicity degree of freedom. The magnetism-driven electric polarization provides a direct handle for the external electrical field to enable the manipulation of helicity with the topological charge conserved. The flip of η between 0 and π accompanied by the reversible P is observed for the PM skyrmions. The more intriguing tunability can be achieved for the IM skyrmions, where the different dependences of P on η are revealed for skyrmions and antiskyrmions respectively. Take this advantage, the helicity can be continuously rotated and exactly picked by applying E in both skyrmion and antiskyrmion cases.
The in-plane uniaxial anisotropy favors this manipulation, whereas the anisotropy perpendicular to the lattice plane weakens the continuity of the tunability. In addition, the crystal form can be controlled easier than the isolated one due to less contribution from the background.