The optical tweezer of skyrmions

In a spin-driven multiferroic system, the magnetoelectric coupling has the form of effective dynamical Dzyaloshinskii–Moriya (DM) interaction. Experimentally, it is confirmed, for instance, for Cu2OSeO3, that the DM interaction has an essential role in the formation of skyrmions, which are topologically protected magnetic structures. Those skyrmions are very robust and can be manipulated through an electric field. The external electric field couples to the spin-driven ferroelectric polarization and the skyrmionic magnetic texture emerged due to the DM interaction. In this work, we demonstrate the effect of optical tweezing. For a particular configuration of the external electric fields it is possible to trap or release the skyrmions in a highly controlled manner. The functionality of the proposed tweezer is visualized by micromagnetic simulations and model analysis.


INTRODUCTION
Optimal dynamical control of a particle motion includes several tasks, such as acceleration, braking, and trapping. In the case of nanoparticles, ions, or atoms, the trapping problem becomes more demanding than the others, except trapping of charged particles which is relatively easy with the use of Pauli trap 1,2 . In the early 90-ties, it was realized that light-atom interaction allows trapping of neutral objects-cesium and sodium atoms in particular 3,4 . In the case of optical trapping of neutral objects, the light does two jobs: (i) it attracts the particles towards the anti-nodal points of maximum intensity of the optical lattice with the spatial period of the order of optical wavelength, and (ii) the light additionally cools down the atoms. The invention of optical tweezers in 1986 by Arthur Ashkin was a triumph for the manipulation of microparticles with laser light 5 . Although trapping of various particles is widely discussed in the literature, the problem of trapping of localized excited modes, especially of topological solitons (skyrmions) has not been studied yet.
The concept of skyrmion traces back to the paper of Skyrme 6 , and to the fundamental paper of Belavin and Polyakov 7 . It is now well known that skyrmion has a topological character. In particular, invariance of the topological action of the field theory, S top n ð Þ ¼ iθ 4π R dx 1 dx 2 n Á ∂ 1 n ∂ 2 n ð Þ , with respect to the infinitesimal transformation n x ð Þ ! n x ð Þ þ ϵ a x ð ÞR a n x ð Þ, (where ϵ a is infinitesimal parameter and R a stands for generators of the O(3) group) defines the specific texture of the vector field n x ð Þ [8][9][10] . The set of different textures of n x ð Þ, obtained from each other by means of the continuous deformation, has the same invariant topological action and the related conserved topological charge W ¼ 1 iθ S top n ð Þ. Thus, one could argue that the topological soliton (skyrmion) is a robust object, stable with respect to small perturbations. Apart from this, skyrmions possess dual fieldparticle properties  . Skyrmions are highly mobile objects. There are several precise recipes on how to drive a skyrmion-either by a spin-polarized electron current or with a magnonic spin current that exerts a magnon pressure on the skyrmion surface. In the recent work 41 , an alternative mechanism of skyrmion drag was proposed, which is based on a combination of uniform temperature profile and non-uniform electric field. Nevertheless, a vital question that arises is whether the particle nature of skyrmions facilitates their trapping. In what follows, we explore trapping of a skyrmion in the laser field e z E ls (x, y, z, t) (with e z being the unit polarization vector of the electric field) and the external electric field E 0 = (0, 0, E z0 ). There exist several methods to manipulate the polarization of the laser beam. Through these methods, the polarization of the electric field can be switched to the desired direction. For example, one can utilize ultrafast timedependent polarization rotation in a magnetophotonic crystal 42 . The colloidal microspheres also can produce dominant E z component 43 .
Skyrmions emerge in materials (e.g., in chiral single-phase multiferroics [44][45][46] ) with a sizeable magnetoelectric (ME) coupling term, is the net ferroelectric polarization, with m denoting the unit vector along the magnetization and c E is the magnetoelectric coupling constant. In chiral multiferroics, the coupling of the external electric field with the ferroelectric polarization mimics the Dzyaloshinskii-Moriya (DM) term and leads to the noncolinear topological magnetic order. The mechanism of trapping of a skyrmion relies on the interaction between the electric component of the laser field and the ferroelectric polarization of the skyrmion texture. As for the specific materials, we focus on two types of materials: spin-driven single-phase multiferroics and Yttrium Iron Garnet (YIG) [11][12][13][14] . In particular, we present in addition to YIG, results for the multiferroic material Cu 2 OSeO 3 , which supports skyrmions. As detailed below, the emergence of the finite (but small) electric polarization due to non-collinearity of the spin allows for the movement of the skyrmions with external electric fields. YIG and single-phase multiferroic are described by free energies densities: Here, M = M s m, where M s is the saturation magnetization, A ex is the exchange stiffness, and H z is the external magnetic field applied along the z-direction. The free energy of the single-phase multiferroic Cu 2 OSeO 3 has the bulk-related DM interaction term The effective magnetic field acting on the magnetization follows from the functional derivative of the free energy functional δm . The interaction energy E me = −E ⋅ P between the external electric fields E and the spin-driven polarization P enters Eq. (1) as a term, which is linear in P. We note that for spindriven multiferroics the spin-induced P is quite small (we recall where c E is related to the spin-orbit coupling and the spatial variations in m are smooth on an atomic scale). Thus, higher-order terms (P) n and the spatial variations (∇P) n , which both account for the energy density of ferroelectric polarization, are negligible and therefore do not appear in the Eq. (1) above.
The laser manipulated skyrmion dynamics is governed by the stochastic Landau-Lifshitz-Gilbert (LLG) equation 47,48 , supplemented by the ME term where γ is the gyromagnetic ratio and α is the phenomenological Gilbert damping constant. The effective field H eff for the singlephase multiferroic consists of the exchange field, DM field, and of the applied external magnetic field, The temperature in the LLG equation, is introduced through the correlation function of the thermal random magnetic field h l , hh l;p ðt; rÞh l;q ðt 0 ; r 0 Þi ¼ 2kBTsimα γμ 0 MsV δ pq δðr À r 0 Þδðt À t 0 Þ, where p, q = x, y, z, k B is the Boltzmann constant, and V is the volume of the single cell, used in numerical simulations. The value of the temperature T sim , we determine from the heat equation (see "Methods" section). We note that the physical temperature and the simulation temperature are related through the equation 48 T sim ¼ Ta sim =a L , where a L is the lattice constant and a sim is the cell length in simulation. Therefore, the physical temperature T = 50 K corresponds to the simulation temperature of T sim % 100 K.
The z component, E z0 , of the external electric field stabilizes the skyrmion structure. Due to the Gaussian profile of the laser field, E ls (x, y, z, t) has the maximum (denoted as E 0 ) in the center of laser spot. The total z component of the electric field, E z = E z0 + E ls (x, y, z, t), is not homogeneous in the (x, y) plane. Depending on the sign of the oscillating laser field E ls (x, y, z, t), the total field E z can be either negative or positive. We note that for an ultrashort laser pulse, the pulse compressor allows control of the spectral phase ϕðωÞ; E ls ðx; y; z; ωÞ ¼ ffiffiffiffiffiffiffiffiffi ffi jE ls j 2 q expðÀiϕðωÞÞ, where ϕðωÞ ¼ À ω c nðωÞd, n(ω) is the index of refraction and d is the film thickness 49 . In what follows, we consider both negative E 0 < 0 and positive E 0 > 0 values of the field. We note that modern laser technologies allow generation of ultrashort single E ls (x, y, z, t) = E ls (x, y, z)f scp (t) and half cycle E ls (x, y, z, t) = E ls (x, y, z)f hcp (t) pulses 50 . The temporal profiles of laser pulses are defined as follows: The ultrashort single pulse has both positive and negative E ls (x, y, z, t), while the negative field part of f hcp (t) is too small. Therefore, for half-cycle pulse E ls (x, y, z, t) can be viewed as positively defined. Before presenting the numerical results, we explain the trapping mechanism. The electric field E z is inhomogeneous only in the x direction. The functional derivative of the ME term with respect to the magnetic moment reads: À 1 Here j = x, y. We focus on the first term fueled by the non-uniform electric field ∂ x E z , while the second term corresponds to the effective DM interaction with a strength tunable by a constant electric field 41 . For tweezing, we suggest using the scanning near-field optical microscopy (SNOM) and advanced nanofabrication procedures. These two methods permit to obtain spots of light 10-20 nm in size; see recent review and references therein 51 . Contribution of the non-uniform electric field will be presented in the form of inhomogeneous electric torque (IET): Àγm À δEmeð∂x EzÞ The vector p E = x × e z is set by e z , which points into the direction of electric field. Obviously, the expression of IET is identical to the standard spin transfer torque-c j m × (m × p), because p E in IET mimics the spin polarization direction p. However, while c j depends on the electric current density, the amplitude of the IET depends on the gradient of the electric field ∂ x E z and on the ME coupling strength c E . In the case of Gaussian laser beam (for more details, we refer to the Supplementary Note 2), the coefficient in the expression for IET, c ¼ γcE ∂r Ez μ 0 Ms , is determined by the gradient of electric field, while is the unit vector. The underlying mechanism of the skyrmion tweezer is as follows: depending on the direction of the laser field, the IET torque is either centripetal (drives the skyrmion to the center of the beam) or countercentripetal (drives the skyrmion out of the beam center).
While the energy is supplied through the laser, skyrmion releases energy to the bulk and SNOM shield. Thus, for the comprehensible study of the skyrmion temperature, one needs to solve the heat equation with source and sink terms included. The maximal temperature of the skyrmion texture can be estimated analytically (see "Methods" section). In our case T max = 50 K and therefore the skyrmion is stable.

Skyrmion motion
We perform numerical simulations based on Eq. (2) for Néel-type skyrmion in YIG, and Bloch-type skyrmion in Cu 2 OSeO 3 stabilized by the constant electric and magnetic fields. In Fig. 1, we illustrate attraction and repulsion mechanisms of the skyrmion tweezer. In the first case, Fig. 1a, the skyrmion is initially embedded at the point (x, y) = (−7.5, 0) nm, and the laser field is positive, E ls (t) > 0. Therefore, p E = −y, c < 0, and the torque winds the skyrmion on clockwise to the laser beam center (0, 0). In the second case, Fig.  1b, the direction of the laser field and IET are reversed, E ls (t) < 0, p E = y, c > 0, and the skyrmion winds out anticlockwise from the laser beam center. In Fig. 1c The strategy for skyrmions trapping is as follows: Focus the laser beam on the center of the skyrmion texture. Steer the center of the beam until the electric field is positive E ls > 0, the skyrmion follows then the center of the beam (Fig. 2). Rotation of skyrmion leads to a weak oscillation of the skyrmion center (q x (t), q y (t)). When the beam velocity v x is below a critical velocity v c x , v x < v c x ¼ 14:22 m s −1 , increase of the beam velocity v x leads to an increase in the velocity of skyrmion drag. When the beam velocity is above v c x , the skyrmion is not able to follow the center of the laser beam (Fig. 3). The critical velocity v c x increases linearly with E 0 , as is demonstrated in Fig. 4a. Thus, one can argue that the skyrmion behaves as a massive object. Changing sign of the laser field from positive to negative, E ls < 0, releases the skyrmion and drives it off the center of the beam (not shown). Skyrmions are topologically protected objects. However, when the ground state is the ferromagnetic state, thermal fluctuations may cause a thermal collapse of the skyrmion. This problem has been widely discussed in the recent literature [52][53][54][55] . The rate of skyrmion collapse follows the Arrhenius law where U is the relevant barrier height. In the vicinity of the critical region, the value of the barrier height can be estimated analytically. At a critical value of the magnetic field, H c % JM 0 D=J ð Þ 4=3 , the radius of the skyrmion starts shrinking. For Bloch-type of skyrmions, the critical field and height of the barrier can be estimated as follows 55 U=JM 2 0 % ðD=JÞ 2=3 1 À H=H c ð Þ 3=2 . The analytical estimation of the barrier is valid only in the critical region. Away from the critical region, the rate of skyrmion collapse can be estimated numerically as 56    For v x = 14.19 m s −1 , the skyrmion center (q x (t), q y (t)) follows the center of the laser beam. When, v x = 14.22 m s −1 , the skyrmion center is able to follow the laser beam only at the beginning of the evolution, t < 5 ns. Other parameters as in Fig. 1. x as a function of the frequency f e plotted for electric fields: E 0 ¼ E l0 þ E l1 sinð2πf e tÞ, with E l0 = 1.2 MV cm −1 and E l1 = 0.07E l0 . c Dependence of the critical velocity on the E l1 /E l0 at the frequency f e = 0.22 GHz. Other parameters as in Fig. 1.
number of collapsed skyrmions at the time t i and N is the total number of skyrmions. Calculations done for T = 50 K, t = 100 ns show that barrier height is U = 386 K and therefore the probability of collapse at T = 50 K is zero. We obtained this information as follows: the time-dependent probability P(t) of skyrmion stability was fitted with the function PðtÞ ¼ expðÀΓtÞ, where Γ is the rate coefficient. The probability P(t) was extracted from the statistics collected through hundreds of repeated simulations. For a given temperature of skyrmion T = 135 K, we simulate the timedependent probability P(t), as demonstrated in Fig. 5a. Through the fitting of statistical P(t) to the formula PðtÞ ¼ expðÀΓtÞ, we estimate the value of the parameter Γ = 0.012 (ns) −1 . In turn, the temperature dependence ΓðTÞ ¼ Γ 0 expðÀU=TÞ follows Arrhenius law, and through the fitting of curves in Fig. 5c, we obtain U = 386 K, Γ 0 = 0.15 (ns) −1 . For these parameters, the probability of stability of skyrmion is P = 0.993 at t = 100 ns and T = 50 K, which confirms that the skyrmion is stable (Fig. 5b). This result is supported by the recent experiment 57 showing the robustness of skyrmions in Cu 2 OSeO 3 .
We also analyzed the influences of oscillating laser electric field E 0 ¼ E l0 þ E l1 sinð2πf e tÞ. It turns out that the oscillating field drags the skyrmion, and the critical velocity v c x as a function of the frequency f e is shown in Fig. 4b. The trapping of the skyrmion depends on the frequency of the field. As we see, the critical velocity v c x drops down at f e = 0.2 GHz. Analyzing the spectrum of the skyrmion oscillation frequency (not shown), we find that the frequency f e = 0.22 GHz coincides with the natural frequency of the laser-induced pinning potential of the skyrmion, i.e., the resonant oscillation frequency of the rigid skyrmion. The resonant amplification of the skyrmion oscillations leads to a release of the skyrmion, and thus reduces v c x . Furthermore, increase of E l1 leads to a decrease of v c x (Fig. 4c). The large E l1 activates nonlinear effects and dependence of the critical velocity on the frequency is not linear anymore, see Fig. 4c for E l1 > 0.05E l0 .
Similar to evanescent Gaussian laser beam, the oscillating laser field also traps the skyrmion. We simulate the laser pulses 70 ps in width and period and steer the center of the laser beam on a distance 14 ffiffi ffi 2 p nm along y = x in 1.3 ns. As we see in Fig. 6, the skyrmion is trapped by the laser beam and follows the center of the laser beam (Supplementary Note 2). The speed of the skyrmion moving along the y = x axis is about 15.5 m s −1 . As we have already mentioned, the obtained results can be also interpreted in terms of the Thiele equation that describes motion of a rigid skyrmion 58,59 , see Fig. 1c, d and the inset to Fig. 2, as well as the Supplementary Fig. 4.
In the single-phase multiferroic Cu 2 OSeO 3 the bulk-type DM interaction 11 stabilizes the Bloch skyrmion (its structure is shown in Fig. 7a). The bulk-type DM interaction contributes to the total magnetic field in LLG equation through the effective field  For v x = 13.09 m s −1 , the skyrmion center (q x (t), q y (t)) follows the center of the laser beam. When, v x = 13.25 m s −1 , the skyrmion center is able to follow the laser beam only at the beginning of the evolution.
In case of YIG, the role of the bulk-type DM interaction is replaced by the constant electric field E z0 = 1.7 MV cm −1 . In this case, the laser electric field has an amplitude of the order of E 0 = 0.2 MV cm −1 . As shown in Fig. 7, when the laser center is static, the IET torque induces the centripetal motion of skyrmion. After steering the center of the laser beam, the skyrmion still follows the center of the beam until the critical velocity 13.09 m s −1 .
Thermal influence of laser radiation A further influence of the laser pulses is heating. The temperature profile T(x, y, t) induced by the laser heating, through the beam with a moving center, is shown in Fig. 8. The temperature T(x, y, t) is calculated from the heat equation (see "Methods" section). The region of the largest temperature (about 50 K) follows the center of the laser beam and the temperature gradually decreases with the distance from the center. The laser heating leads to the inhomogeneous time-dependent temperature profile and affects the skyrmion dynamics. The main effect of laser-induced heating is that the trapping process becomes non-deterministic. We performed a set of calculations with the same initial conditions and collect ensemble statistics in Fig. 9. The trapping probability P decreases for the higher velocity of the center of beam, but remains finite. Even at the speed v x = 14.32 m s −1 and E 0 = 1.2 MV cm −1 , P = 54% as is shown in the inset in Fig. 9. An interesting fact is that below the threshold velocity of v x ≤ 14 m s −1 the probability P = 1.
In summary, we have proposed a method of optical control of skyrmions in magnetoelectric materials. Owing to the magnetoelectric coupling, electric field of a laser beam couples to the magnetic moments of the skyrmion. When the field in the laser beam is pre-designed in an appropriate way, one can trap, shift, and then release the skyrmion. Such an optical tweezer may be very useful in control and manipulation of the skyrmion position. Numerical results have been obtained from micromagnetic simulations based on Landau-Lifshitz-Gilbert equation with a contribution from magnetoelectric coupling and additionally from solution of the Thiele equation describing motion of rigid skyrmions. A very good agreement of the results obtained by these two methods has been achieved.

Micromagnetic modeling
To explore numerically the skyrmion dynamics, we utilize Eq. (2). The simulations have been done for the magnetic film (which is the x − y plane) with the size of 500 nm × 500 nm × 2.5 nm. The film is discretized into the simulation cells of 2.5 nm × 2.5 nm × 2.5 nm (i.e., a sim ¼ 2:5 nm). The cell size is smaller than the typical exchange length ( In the simulation, we first let a single skyrmion to relax to the stationary state and only after apply the laser field. The position of the skyrmion center (q x , q y ) is measured through the formulas q x = (∫xqdxdy)/ (∫qdxdy) and q y = (∫yqdxdy)/(∫qdxdy) and the simulation data. Here qðx; yÞ ¼ 1 4π m Á ð∂ x m ∂ y mÞ is the topological charge density.

Heat equation
The temperature profile T(x, y, t) is the solution of the heat equation: ∂Tðx; y; tÞ ∂t ¼ D∇ 2 Tðx; y; tÞ þ Iðx; y; tÞ À bTðx; y; tÞ: Here D ¼ kph ρC is the thermal diffusivity. In simulation, we considered the parameters: k ph = 6 W (m K) −1 is thermal conductivity, ρ = 5170 kg m −3 is the mass density, and C = 570 J (kg K) −1 is the heat capacity. The source term is where c is the light speed, ϵ 0 is the permittivity of vacuum. We exploit axial symmetry and rewrite source in the form I ¼ I 0 exp À ϱ 0 2 σ 2 0 . The last term bT (x, y, t) has the role of the sink and describes heat exchange of the skyrmion with the rest of the system, i.e., substrate and SNOM shield. Its value can be estimated from the heat diffusivity of the substrate 60 b % 3 ffiffi π p Ds σ 2 0 . In particular, for D s = 0.24 cm 2 s −1 , we deduce b = 1.6 × 10 10 s −1 .
The parameter δ T = 1.5 × 10 6 m −1 describes the laser penetration depth and is proportional to the complex part of the refractive index. The absorption efficiency of laser energy l x is equal to the real part of the refractive index. We note that for YIG, both real and complex parts drastically depend on the frequency and magnetic field and, in certain regimes, are rather small 61 . Thus, parameter l x should be reasonably small.
On the other hand, we cannot precisely define the value of l x and consider it as a phenomenological parameter with particular attention to the value l x = 8%. The result for l x = 5%, with a lower temperature T = 30 K, leads to the slightly larger skyrmion trapping probability, see inset in Fig. 9. In simulations, we see that skyrmion temperature linearly increases with l x and for l x = 17%, the temperature is still below T < 100 K. However, the skyrmion temperature depends on the temperature of the substrate and SNOM shield (they play the role of sink). If their temperature is zero, then the maximal temperature of the skyrmion in simulations is below T < 25 K even for l x = 100%. Exploiting the Green function method, we solve the heat equation analytically 62 : G ρ À ρ 0 ; t À t 0 ð Þ dθ 0 ϱ 0 dϱ 0 dt 0 ; Tðϱ; t ¼ 0Þ ¼ 0; Tðρ ! 1; tÞ ¼ 0; G ρ À ρ 0 ; t À t 0 ð Þ ¼ e ÀbðtÀt 0 Þ 4πDðtÀt 0 Þ exp À jρÀρ 0 j 2 4DðtÀt 0 Þ n o : The exact analytic solution of Eq. (4) can be obtained for the center of the beam, i.e., the region of the trapped skyrmion.
Here Ei(⋅) is the exponential integral. The analytic solution shows the dependence of the skyrmion temperature on the phenomenological coupling constant b. The maximal temperature of the skyrmion reads T max = I 0 /b. In our calculations, T max = 50 K. The temperature in our case conforms to the temperature in the experiment 57 . The temperature effect can be included in the LLG equation Eq. (1) through the random magnetic field h th , and its correlation function hh l;p ðt; rÞh l;q ðt 0 ; r 0 Þi ¼ 2kBTsimα γμ 0 MsV δ pq δðr À r 0 Þδðt À t 0 Þ. Here, k B is the Boltzmann constant, and V is the volume of the single cell, used in numerical simulations. However one should remember the relation between the physical temperature and the temperature used in the simulations 48 T sim ¼ Ta sim =a L , where a L is the lattice constant and a sim is the cell length in simulation, the lattice constant for YIG a L = 12 Å.
In numerical simulation, temperature is calculated from the heat equation Eq. (3) and Forward-Time Central-Space method 63 . The thermal field in the numerical simulation is determined from h l ¼ η where η is a random vector drawn from a standard normal distribution whose value is changed at every time step, and V ¼ a 3 sim is the volume of a single finite-difference cell.