Thermally assisted Skyrmion drag in a nonuniform electric field

TL;DR

This paper introduces a novel method to manipulate magnetic skyrmions using inhomogeneous electric fields combined with thermal bias, enhancing their velocity for potential spintronics applications.

Contribution

It proposes a new approach for skyrmion control via electric and thermal stimuli, expanding the toolkit beyond current methods like charge and magnonic currents.

Findings

Inhomogeneous electric torque effectively drags skyrmions.

Thermal bias enhances skyrmion velocity.

Method is demonstrated for insulating magnetic materials.

Abstract

Magnetic skyrmions are topologically protected excitations of the magnetization vector field with promising applications in spintronics and spin-caloritronics, particularly due to their high mobility. Skyrmions can be steered by a spin-polarized charge current or by exposure to a magnonic spin current. Here, we propose a further method for driving skyrmions by applying an inhomogeneous electric field and a homogeneous thermal bias. We show that the inhomogeneous electric torque leads to an efficient skyrmionic drag which can be thermally assisted as to enhance the skyrmion velocity. The calculations and analysis are limited to insulating samples; for conducting materials the influence of the inhomogeneous electric field on the charge carriers need to be taken also into account.

Figures

Equations

$$\displaystyle\frac{\partial\mbox{\boldmath$\mathrm{M}$}}{\partial t}=-\gamma\mbox{\boldmath$\mathrm{M}$}\times\bigg{(}\mbox{\boldmath$\mathrm{H}$}_{\rm\mathrm{eff}}+\mbox{\boldmath$\mathrm{h}$}_{\rm l}-\frac{1}{\mu_{\rm 0}M_{\rm s}}\frac{\delta E_{\rm ele}}{\delta\mbox{\boldmath$\mathrm{m}$}}\bigg{)}+\frac{\alpha}{M_{\rm s}}\mbox{\boldmath$\mathrm{M}$}\times\frac{\partial\mbox{\boldmath$\mathrm{M}$}}{\partial t}.$$

$$\displaystyle-\frac{1}{\mu_{\rm 0}M_{\rm s}}\frac{\delta E_{\rm ele}(E_{\rm i})}{\delta\mbox{\boldmath$\mathrm{m}$}}=$$

$$\displaystyle+\sum_{\rm j(j\neq i)}2E_{\rm i}(-\partial_{\rm j}m_{\rm j}\mbox{\boldmath$\mathrm{e}$}_{\rm i}+\partial_{\rm j}m_{\rm i}\mbox{\boldmath$\mathrm{e}$}_{\rm j})].$$

$$\displaystyle-\gamma\mbox{\boldmath$\mathrm{m}$}\times\bigg{(}-\frac{\delta E_{\rm ele}(\partial_{\rm x}E_{\rm i})}{\mu_{\rm 0}M_{\rm s}\delta\mbox{\boldmath$\mathrm{m}$}}\bigg{)}=-\frac{\gamma c_{\rm E}\partial_{\rm x}E_{\rm i}}{\mu_{\rm 0}M_{\rm s}}\mbox{\boldmath$\mathrm{m}$}\times(\mbox{\boldmath$\mathrm{m}$}\times\mbox{\boldmath$\mathrm{p_{\rm E}}$}).$$

$$\displaystyle E_{\rm ele}(\mbox{\boldmath$\mathrm{S}$}^{n})=-\mbox{\boldmath$\mathrm{E}$}^{n-1}\cdot\mbox{\boldmath$\mathrm{P}$}^{n-1}-\mbox{\boldmath$\mathrm{E}$}^{n}\cdot\mbox{\boldmath$\mathrm{P}$}^{n}.$$

$$\displaystyle E_{\rm ele}(\mbox{\boldmath$\mathrm{S}$}^{n},E_{\rm y})$$

$$-\Delta E_{\mathrm{y}}(S_{\mathrm{x}}^n{\bar{S}}_{\mathrm{y}}^n-S_{\mathrm{y}}^n{\bar{S}}_{\mathrm{x}}^n)],$$

$$\displaystyle E_{\rm ele}(\mbox{\boldmath$\mathrm{S}$}^{n},E_{\rm z})$$

$$-\Delta E_{\mathrm{z}}(S_{\mathrm{x}}^n{\bar{S}}_{\mathrm{z}}^n-S_{\mathrm{z}}^n{\bar{S}}_{\mathrm{x}}^n)].$$

$$\displaystyle-\mbox{\boldmath$\mathrm{S}$}^{n}\times\bigg{(}-\frac{\delta E_{\rm ele}^{n}}{\delta\mbox{\boldmath$\mathrm{S}$}^{n}}\bigg{)}=-\frac{eaJ\Delta E_{\rm i}}{E_{\rm so}}\mbox{\boldmath$\mathrm{S}$}^{n}\times[\langle\mbox{\boldmath$\mathrm{S}$}^{n}\rangle\times\mbox{\boldmath$\mathrm{p_{\rm E}}$}],$$

$$\frac{\partial m_{\mathrm{x}}}{\partial x}{|}_{\partial \mathrm{V}}+\frac{c_{\mathrm{E}}{E}_{\mathrm{y}}}{2A_{\mathrm{ex}}}{m}_{\mathrm{y}}+\frac{c_{\mathrm{E}}{E}_{\mathrm{z}}}{2A_{\mathrm{ex}}}{m}_{\mathrm{z}}$$

$$\frac{\partial m_{\mathrm{y}}}{\partial x}{|}_{\partial \mathrm{V}}-\frac{c_{\mathrm{E}}{E}_{\mathrm{y}}}{2A_{\mathrm{ex}}}{m}_{\mathrm{x}}$$

$$\frac{\partial m_{\mathrm{z}}}{\partial x}{|}_{\partial \mathrm{V}}-\frac{c_{\mathrm{E}}{E}_{\mathrm{z}}}{2A_{\mathrm{ex}}}{m}_{\mathrm{x}}$$

$$\frac{\partial m_{\mathrm{x}}}{\partial y}{|}_{\partial \mathrm{V}}-\frac{c_{\mathrm{E}}{E}_{\mathrm{x}}}{2A_{\mathrm{ex}}}{m}_{\mathrm{y}}$$

$$\frac{\partial m_{\mathrm{y}}}{\partial y}{|}_{\partial \mathrm{V}}+\frac{c_{\mathrm{E}}{E}_{\mathrm{x}}}{2A_{\mathrm{ex}}}{m}_{\mathrm{x}}+\frac{c_{\mathrm{E}}{E}_{\mathrm{z}}}{2A_{\mathrm{ex}}}{m}_{\mathrm{z}}$$

$$\frac{\partial m_{\mathrm{z}}}{\partial y}{|}_{\partial \mathrm{V}}-\frac{c_{\mathrm{E}}{E}_{\mathrm{z}}}{2A_{\mathrm{ex}}}{m}_{\mathrm{y}}$$

$$\displaystyle\mbox{\boldmath$\mathrm{T}$}_{\rm(nx=1)}$$

$$\displaystyle\mbox{\boldmath$\mathrm{T}$}_{\rm(nx=n_{\rm 0})}$$

$$\hat{H}={\hat{H}}_{\mathrm{ex}}+{\hat{H}}_{\mathrm{D}},$$

$${\hat{H}}_{\mathrm{ex}}=J\sum _{\mathrm{n}}{\hat{S}}_{\mathrm{n}}{\hat{S}}_{\mathrm{n}+1},$$

$$\displaystyle\hat{H}_{\rm D}=-2D_{\rm z}\sum_{\rm n}\hat{z}\cdot\big{(}\hat{S}_{\rm n}\times\hat{S}_{\rm n+1}\big{)}-$$

$$\displaystyle 2\sum_{\rm n}D_{\rm y,n}\hat{y}\cdot\big{(}\hat{S}_{\rm n}\times\hat{S}_{\rm n+1}\big{)}.$$

$$\displaystyle\frac{d}{dt}\langle\hat{H}_{\rm mag}\rangle=i[\hat{H}_{\rm D},\hat{H}_{\rm mag}]=4S^{2}J\times\big{\langle}\sum_{\rm n}i\big{(}\hat{a}^{{\dagger}}_{\rm n}-\hat{a}_{\rm n}\big{)}\cdot$$

$$\displaystyle\big{(}D_{\rm y,n}\hat{a}^{{\dagger}}_{\rm n+1}\hat{a}_{\rm n+1}-D_{\rm y,n-1}\hat{a}^{{\dagger}}_{\rm n-1}\hat{a}_{\rm n-1}\big{)}\big{\rangle}.$$

Full text

I Introduction

The notion of skyrmion was introduced in high-energy physics by Skyrme Skyrme as a topologically stable vector field configuration. The vector order parameter of a magnetically ordered system also shows such stable excitations, as experimentally observed in the chiral itinerant compound MnSi Binz . Essential for the formation of magnetic skyrmions is the absence of inversion symmetry Iwasaki ; Dai ; HoonHan ; Bogdanov ; Garst ; Gavilano , that can be imposed for instance by the Dzyaloshinskii-Moriya interaction. Skyrmions exhibit a quasi-particle character with mobility higher than a magnetic domain wall Bulaevskii ; Kong ; Mishra ; Batista ; Hoogdalem ; Papanicolaou . Besides, they are less pinned to dislocations and impurities and may have an extension on the nanometer scale, which make them interesting for spintronic applications Saxena ; Rosch ; Finocchio . An important issue thereby is the controlled driving of skyrmions. One may set skyrmion in motion using a spin-polarized charge current Jonietz ; Iwasaki1463 ; yzhou2016 or as a thermally assisted magnonic current Batista ; Kong ; Mochizuki2014 .

The pressure exerted by spin transfer torque (STT) on the surface of skyrmion moves the skyrmion and allows achieving reasonably high speed. The magnonic spin Seebeck current exploits the presence of a thermal gradient which is usually difficult to realize and hard to control swiftly on the nanoscale. To generate STT, various types of methods are employed, including the spin Hall effect in heavy metal/magnet heterostructure Liu2012 ; sinova2015 , spin injection in nonlocal structure Demidov2016 ; Lin5177 , and STT in metallic multilayer Slonczewski ; Berger . Each of these methods provides a unique insight into the spin transfer torque manipulation. Nevertheless, standard methods require demanding preparation of the sample, and their application is usually hampered by Joule heating or interface effect.

In the present work, we show that at a finite uniform temperature, applied inhomogeneous electric field generates an inhomogeneous electric torque (IET) that can be used to drive skyrmions. Such fields are widely available, and their spatiotemporal structure is accurately controllable. The proposed scheme does not rely on dissipative charge currents and hence involves less energy dissipation. In addition, electric fields are more easily spatio-temporally engineered and controlled on the nanoscale than temperature gradients that might be also used to move skyrmions. One should note however, that in principle, an inhomogeneous electric field applied to an electrical conductor induces a local inhomogeneous Joule heating and through the formed inhomogeneous temperature profile can influence the magnon density. Thus, the inhomogeneous Joule heating has an extra effect that can additionally influence the skyrmion drag but further dissipates energy. This effect is absent for insulating materials to which the current study is restricted.

The paper is organized as follows: in Sec. II we specify the model, in Sec. III we explore the mechanisms of the formation of IET, in Sec. IV we study the inhomogeneous damping and the intrinsic frequencies of the system. The thermally assisted magnonic current is addressed in the Sec. V, in the Sec. VI we analyze the skyrmion motion induced by the IET. The final Sec. VII is dedicated to the effect of the time-dependent electric field.

II Model

In spite of the absence of itinerant electrons in single phase multiferroic or magnetic insulators, the virtual hopping of electrons between the $d$ orbitals and the strong spin-orbit interaction leads to a net ferroelectric polarization \mbox{\boldmath\mathrm{P}}=c_{\rm E}[(\mbox{\boldmath\mathrm{m}}\cdot\nabla)\mbox{\boldmath\mathrm{m}}-\mbox{\boldmath\mathrm{m}}(\nabla\cdot\mbox{\boldmath\mathrm{m}})]. Here, $\mathrm{m}$ is the unit magnetization vector of the magnet, and $c_{\rm E}$ is a magneto-electric (ME) coupling constant (see Tliu20 ; VRis for details). The net ferroelectric polarization is coupled to the applied external electric field \mbox{\boldmath\mathrm{E}}=(E_{\rm x},E_{\rm y},E_{\rm z}). The ME coupling energy E_{\rm ele}=-\mbox{\boldmath\mathrm{E}}\cdot\mbox{\boldmath\mathrm{P}} mimics that resulting from a dynamical Dzyaloshinskii-Moriya (DM) interaction and leads to the formation of N$\mathrm{\acute{e}}$el-type skyrmion . The $E_{\rm z}$ component of the field, applied along the whole sample stabilizes the skyrmion Troncoso2016 ; loidl2015 ; Tokura2013 . The same skyrmion structure can be stabilized utilizing the bulk (in the material with broken inversion symmetry in lattices) or interface type DM interaction term (at the interface of magnetic films) Tokura2013 .

For a uniform finite temperature the magnetization dynamics is governed by the stochastic Landau-Lifshitz-Gilbert (LLG) equation supplemented by the ME term

[TABLE]

Here, \mbox{\boldmath\mathrm{M}}=M_{\rm s}\mbox{\boldmath\mathrm{m}} and $M_{\rm s}$ is the saturation magnetization, $\gamma$ is the gyromagnetic ratio, and $\alpha$ is the phenomenological Gilbert damping constant. The effective field \mbox{\boldmath\mathrm{H}}_{\rm\mathrm{eff}} consists of the exchange field and the applied external magnetic filed, \mbox{\boldmath\mathrm{H}}_{\rm\mathrm{eff}}=\frac{2A_{\rm ex}}{\mu_{\rm 0}M_{\rm s}}\nabla^{2}\mbox{\boldmath\mathrm{m}}+H_{\rm z}\mbox{\boldmath\mathrm{z}}, where $A_{\rm ex}$ is the exchange stiffness, and $H_{\rm z}$ is the external magnetic field applied along the z-direction. The thermal random magnetic field is characterized by the correlation function of a white noise JLGar22 , meaning \langle\mbox{\boldmath\mathrm{h}}_{\rm l,i}(x,t)\mbox{\boldmath\mathrm{h}}_{\rm l,j}(x^{\prime},t^{\prime})\rangle=\frac{2k_{\rm B}T\alpha}{\gamma M_{\rm s}V}\delta_{\rm ij}\delta(x-x^{\prime})\delta(t-t^{\prime}) , where $k_{\rm B}$ is the Boltzmann constant and $T$ is the temperature.

Existence of skyrmions can be verified through the topological number $Q=\int dxdy\rho_{\rm\mathrm{sky}}$, with \rho_{\rm\mathrm{sky}}=-\frac{1}{4\pi}\mbox{\boldmath\mathrm{m}}\cdot(\partial_{\rm x}\mbox{\boldmath\mathrm{m}}\times\partial_{\rm y}\mbox{\boldmath\mathrm{m}}). In the case of a single skyrmion, we obtain $Q=1$. The inhomogeneous electric (E) field $\left|E_{\rm y}\right|$ is applied only to a part of the system. This can be achieved for instance a screening metalic layer to block the E-field from the respective region. In what follows we show that the inhomogeneous E-field $E_{\rm y}$ modifies the thermal magnon density profile and the induced magnon flow drives the skyrmion. The numerical simulations based on Eq. (1) are done at zero and finite temperatures for the saturation magnetization $M_{\rm s}=1.4\times 10^{5}$ A/m, the exchange constant $A_{\rm ex}=3\times 10^{-12}$ J/m, the ME coupling strength $c_{\rm E}=0.9$ pC/m, and the Gilbert damping constant $\alpha=0.001$. The N$\mathrm{\acute{e}}$el-type skyrmion is stabilized by the electric and magnetic fields $E_{\rm z}=1.7$ MV/cm, $H_{\rm z}=3.2\times 10^{5}$ A/m. The transversal component of the electric field $\left|E_{\rm y}\right|$ is in the order of (0, 0.15 MV/cm): This value of $\left|E_{\rm y}\right|$ is small enough and cannot induce switching of the equilibrium magnetization.

III Inhomogeneous electric torque

The term \frac{1}{\mu_{\rm 0}M_{\rm s}}\frac{\delta E_{\rm ele}}{\delta\mbox{\boldmath\mathrm{m}}} that enters in the LLG equation Eq.(1) quantifies the influence of the ME coupling on the effective field. Thus, the ME coupling constant $c_{\rm E}$ and the configuration of the applied electric field are important issues to consider. We investigated inhomogeneous electric field \mbox{\boldmath\mathrm{E}}(x) (varying along the $x$ axis) leading to

[TABLE]

Here, $i,j=x,y,z$. The total ME torque term that enters the LLG equation -\gamma\mbox{\boldmath\mathrm{m}}\times[-\frac{1}{\mu_{\rm 0}M_{\rm s}}\frac{\delta E_{\rm ele}}{\delta\mbox{\boldmath\mathrm{m}}}] has two sources: the electric field \mbox{\boldmath\mathrm{E}}=(E_{\rm x},E_{\rm y},E_{\rm z}) and its gradient \partial_{\rm x}\mbox{\boldmath\mathrm{E}}=(\partial_{\rm x}E_{\rm x},\partial_{\rm x}E_{\rm y},\partial_{\rm x}E_{\rm z}). The inhomogeneity of the electric field is manifested in the spatially-inhomogeneous DM interaction and in the additional terms $\partial_{\rm x}E_{\rm i}$, where $i=x,y,z$. After some algebra, we infer the expression for the IET solely induced by the electric field gradient \partial_{\rm x}\mbox{\boldmath\mathrm{E}}

[TABLE]

The vector \mbox{\boldmath\mathrm{p_{\rm E}}}=\mbox{\boldmath\mathrm{x}}\times\mbox{\boldmath\mathrm{e}}_{\rm i} is set by \mbox{\boldmath\mathrm{e}}_{\rm i} which points into the direction of electric field. Obviously the expression Eq.(3) is identical to the standard spin transfer torque -c_{\rm j}\mbox{\boldmath\mathrm{m}}\times(\mbox{\boldmath\mathrm{m}}\times\mbox{\boldmath\mathrm{p}}), in which case $\mathrm{p_{\rm E}}$ mimics the spin polarization direction $\mathrm{p}$. While $c_{\rm j}$ depends on the electric current density, the amplitude of the IET depends on the gradient of an electric field $\partial_{\rm x}E_{\rm i}$ and the ME coupling strength $c_{\rm E}$.

For further insight let us utilize an atomistic model in which the polarization vector $\mathrm{P}$ is expressed as \mbox{\boldmath\mathrm{P}}=-\frac{Jea}{E_{\rm so}}\mbox{\boldmath\mathrm{e}}_{\rm n,n+1}\times(\mbox{\boldmath\mathrm{S}}^{n}\times\mbox{\boldmath\mathrm{S}}^{n+1}) , where \mbox{\boldmath\mathrm{S}}^{n} is the spin localized on the n-th site, the strength of the spin-orbit interaction is quantified by $E_{\rm SO}=\hbar^{2}/(2m_{\rm e}\lambda^{2})$, $m_{\rm e}$ is the electron mass, $\lambda$ is the spin-orbit coupling constant, $J$ is the exchange coefficient, $a$ is the distance between neighboring magnetic ions, $e$ is the electron charge, \mbox{\boldmath\mathrm{e}}_{\rm n,n+1} is the unit vector connecting the ions. The ME coupling term E_{\rm ele}(\mbox{\boldmath\mathrm{S}}^{n}) in the atomistic model can be rewritten in the form

[TABLE]

We adopt the geometry of 1D chain stretched out along the $x$ axis and for the electric field components $E_{\rm y},~{}E_{\rm z}$ we write down the explicit expressions of ME terms:

[TABLE]

Here, we implemented the ansatz: \mbox{\boldmath\mathrm{E}}_{\rm n-1}=\mbox{\boldmath\mathrm{E}}_{\rm 0}-\Delta\mbox{\boldmath\mathrm{E}}/2, \mbox{\boldmath\mathrm{E}}_{\rm n}=\mbox{\boldmath\mathrm{E}}_{\rm 0}+\Delta\mbox{\boldmath\mathrm{E}}/2, and \mbox{\boldmath\mathrm{S}}^{n-1} = \langle\mbox{\boldmath\mathrm{S}}^{n}\rangle-\Delta\mbox{\boldmath\mathrm{S}}/2, \mbox{\boldmath\mathrm{S}}^{n+1} = \langle\mbox{\boldmath\mathrm{S}}^{n}\rangle+\Delta\mbox{\boldmath\mathrm{S}}/2, where \langle\mbox{\boldmath\mathrm{S}}^{n}\rangle means the averaged value of the spin. The atomistic IET has the form:

[TABLE]

where, \mbox{\boldmath\mathrm{p_{\rm E}}}=\mbox{\boldmath\mathrm{e}}_{\rm n,n+1}\times\mbox{\boldmath\mathrm{e}}_{\rm i} and \mbox{\boldmath\mathrm{e}}_{\rm n,n+1} is the unit vector connecting the sites $n,n+1$. In the continuous limit, the atomistic model Eq.(6) goes over into the continuum model upon coarse graining Eq. (3). We use the following boundary conditions for two-dimensional model Rohart2013

[TABLE]

Here, $E_{\rm x}$, $E_{\rm y}$ and $E_{\rm z}$ are boundary values of the electric field components. Derivatives in Eq. (1) are implemented in the sense of central derivatives, i.e., \partial{\mbox{\boldmath\mathrm{m}}}/\partial i=(\mbox{\boldmath\mathrm{m}}_{\rm ni+1}-\mbox{\boldmath\mathrm{m}}_{\rm ni-1})/(2\triangle_{\rm i}), where $i=x,y$, and $\triangle_{\rm i}$ is the step size along the $i$ direction. Thus, beyond the total ME torque (including IET), the boundary conditions impose extra torque \mbox{\boldmath\mathrm{T}}_{\rm\mathrm{boundary}}. However, while ME torque influences the magnetization dynamics in the whole system, the extra torque imposed by boundary conditions is local and acts in the vicinity of boundaries. In what follows the boundary torque we call the ”boundary magnetoelectric torque.”

IV damping

From an energy point of view, the ME interaction provides an additional route for energy transfer. The ME interaction affects the magnetization precession as well as the damping constant $\alpha_{\rm\mathrm{eff}}$. We inspect these effects based on the magnetic resonance analysis. The time dependence of the magnetization \mbox{\boldmath\mathrm{M}}(t) with or without the electric field is calculated based on the LLG equation (Eq. 1). For $H_{\rm z}=1.6\times 10^{5}$ and a zero E field \mbox{\boldmath\mathrm{E}}=(0,0,0), we estimate the resonance frequency and the relaxation time, respectively $\omega_{\rm 0}=\gamma H_{\rm 0}=2\pi\times 6$ GHz, $\tau_{\rm 0}=1/(\alpha\omega_{\rm 0})=26.5$ ns. Note that, for the adopted excitations, only several intrinsic eigenmodes can be excited, for example, the ferromagnetic resonance and standing wave modes.

As is shown in Fig. 1(a), only the single frequency $\omega=\omega_{\rm 0}=6$ GHz magnetization oscillation is excited in the absence of the external electric field and the corresponding effective damping is constant $\alpha_{\rm\mathrm{eff}}=\alpha$. However, applying a linear electric field \mbox{\boldmath\mathrm{E}}=E_{\rm gr}x\mbox{\boldmath\mathrm{e}}_{\rm y} with $E_{\rm gr}=-0.24$ (MV/m)/nm along the $y$-direction (i.e., \mbox{\boldmath\mathrm{p}}_{\rm E}=\mbox{\boldmath\mathrm{e}}_{\rm z}), changes the free-energy landscape allowing the activation of several high frequency modes, as shown by Fig. 1(b). A detailed analysis suggests that four frequencies $\omega_{\rm 1}=2\pi\times 5.99$ GHz, $\omega_{\rm 2}=2\pi\times 6.75$ GHz, $\omega_{\rm 3}=2\pi\times 9$ GHz and $\omega_{\rm 4}=2\pi\times 12.73$ GHz dominate in Fig. 1(b). The spatial distributions of the dominant frequencies are presented in Fig. 2, and one can see that the dominant frequencies correspond to the standing waves with $2n$ nodes ($n=0,1,2$ and 3) along the $x$ axis. Using the fitting equation $M_{\rm x}(t)=\sum_{\rm i=1,2,3,4}A_{\rm i}\cos(\omega_{\rm i}t+\phi_{\rm i})\exp(-t/\tau_{\rm i})$, we determine $\tau_{\rm 1}=26.6$ ns ($\alpha_{\rm\mathrm{eff}1}=0.00099$), $\tau_{\rm 2}=24.9$ ns ($\alpha_{\rm\mathrm{eff}2}=0.00095$), $\tau_{\rm 3}=20.9$ ns ($\alpha_{\rm\mathrm{eff}3}=0.00085$) and $\tau_{\rm 4}=16.5$ ns ($\alpha_{\rm\mathrm{eff}4}=0.00076$). The applied uniform electric field induces these oscillation modes with $2n$ nodes and leaves the effective damping $\alpha_{\rm\mathrm{eff}}$ unchanged. However, the effective damping $\alpha_{\rm\rm{eff}}$ depends linearly on the gradient of electric field $\partial_{\rm x}E_{\rm y}$, see Fig. 3. Therefore, the effective damping can be controlled through an inhomogeneous electric field.

As distinct from the current-induced Slonczewski torque, the ME coupling effects the magnetization precessional damping not only through the magnetoelectric torque term, i.e., Eq. (3) (the torque related to the \partial_{\rm i}\mbox{\boldmath\mathrm{m}} in Eq. (2) is rather weak for the resonant mode). The ME coupling impacts the boundary constraint and indirectly influences the effective damping. We remove the IET part and write out the ME coupling induced magnetic torque under the boundary constraint (Eq. (7)):

[TABLE]

The equation (8) is expressed in the central derivative form and $nx=1$ and $nx=n_{\rm 0}$ represent the cells located at the left ($x=-125$ nm) and right ($x=125$ nm) boundaries, respectively. Besides the IET (Eq. (3)), obviously there is another additional damping-like torque associated with the effect of boundaries. This additional damping-like torque also modifies the magnetization dynamics. The effective damping $\alpha_{\rm\mathrm{eff}}$ is determined by the context of Eqs.(3) and (8). For the electric field profile \mbox{\boldmath\mathrm{E}}=(0,E_{\rm gr}x\mbox{\boldmath\mathrm{e}}_{\rm y},0), adopted in our simulations, the directions of boundary and IET torques in Eqs.(3) and (8) are always opposite. For the $n=0$ resonance mode, two neighboring spins are always parallel ( \mbox{\boldmath\mathrm{m}}_{\rm(nx=1)}\approx\mbox{\boldmath\mathrm{m}}_{\rm(nx=2)} and \mbox{\boldmath\mathrm{m}}_{\rm(nx=n_{\rm 0})}\approx\mbox{\boldmath\mathrm{m}}_{\rm(nx=n_{\rm 0}-1)} ) in Eq. (8). The Fig. 3 demonstrates that for the $n=0$ resonance mode, two opposite torques (boundary torque and IET) totally compensate each other. For higher standing wave nodes and frequencies, the oscillation of non-collinearity between two neighboring spins increases. This relaxes the boundary torque in Eq. (8) and as a consequence the boundary torque does not compensate IET any more. Therefore, the change in $\alpha_{\rm eff}$ is enhanced for higher resonance modes. We plot the variation of the effective damping constant $\alpha_{\rm eff}$ in the presence and absence of the boundary torque, see Fig. 3 and Fig. 4 respectively. Variation of the effective damping constant $\alpha_{\rm eff}$ in the absence of the boundary torques Fig. 4 is caused solely by IET. Comparing Fig. 3 and Fig. 4 we see that the boundary torque reduces the variation of the effective damping constant $\alpha_{\rm eff}$. It is instructive to explore the dependence of the effective damping constant $\alpha_{\rm eff}$ on the inverse frequency $1/\omega$. The result is presented in Fig. 5. Evidently, the effective damping constant $\alpha_{\rm eff}$ depends linearly on the inverse frequency. This result is similar to the case of the current-induced spin-transfer torque Krivorotov2005 . The wavelength and the frequency of the oscillation mode can be steered with the magnet’s length $l$. This also changes the effective damping. As is shown in Fig. 6(a), with the increase of $l$, the oscillation frequency of $n=1$ mode decreases. Similar to Fig. 5(a), the decreasing of the frequency decreases the variation of the effective damping (Fig. 6(b)).

Furthermore, employing a large enough negative $\partial_{\rm x}E_{\rm y}$, we successfully reverse the magnetization direction and drive the auto-oscillation. This fact testifies the feature of the spin transfer like torque. The averaged $M_{\rm z}$ as a function of $\partial_{\rm x}E_{\rm y}$ is shown in Fig. 7(a). The external magnetic field $H_{\rm z}=1.6\times 10^{5}$ is applied along the $z$ axis and the initial magnetization is along the $+z$. When the negative $\partial_{\rm x}E_{\rm y}$ is large enough, the $+z$ local magnetization becomes unstable after 90 ns of evolution. For $\partial_{\rm x}E_{\rm y}=-1.2$ (MV/m)/nm, the averaged $M_{\rm z}$ reveals the nature of magnetization oscillation induced by IET, as demonstrated in Fig. 7(b). Besides, by setting an anisotropy field \mbox{\boldmath\mathrm{H}}_{\rm k}=2K_{\rm z}m_{\rm z}/(\mu_{\rm 0}M_{\rm s})\mbox{\boldmath\mathrm{e}}_{\rm z} with $K_{\rm z}=2.3\times 10^{4}$ J/m3 instead of $H_{\rm z}$, the direction of the equilibrium magnetization can be reversed by IET, as is shown in Fig. 8, and there is no auto-oscillation in this case.

V thermal magnonic spin current

Conjointly with the uniform temperature bias, STT induces a non-equilibrium magnon flow, meaning a thermally assisted magnonic spin current wang2017 . The thermal random magnetic field activates the magnetization oscillations around the equilibrium state and excites the thermal magnons. In this way, the electric spin-polarized current can be converted into a magnonic current employing costly procedure. To assess the costs, one should include ohmic losses in the generation of the electric spin-polarized current and estimate the current conversion efficiency. In the present work, we show that the IET can be used to generate non-equilibrium magnonic flow. The underlying physical mechanism of our method is based on an inhomogeneous electric field $E_{\rm y}$ that modifies the thermal magnon density profile.

Fig. 9 shows the spatial distribution of the magnon density $\rho$ and the magnonic spin current $J_{\rm ex}$. The temperature $T$ is taken to be 25 K, the component of the electric field $E_{\rm y}=\pm 0.15\mathrm{MV/cm}$ is applied only in the region $x<x_{\rm 0}=0$. The magnon density $\rho=M_{\rm x}^{2}+M_{\rm y}^{2}$ is quantified by the squared transversal magnetization components and the total magnonic spin current is equal to $J_{\rm s}=J_{\rm ex}+J_{\rm E}$, where $J_{\rm ex}=\frac{2\gamma A_{\rm ex}}{\mu_{\rm 0}M_{\rm s}^{2}}(M_{\rm x}\partial_{\rm x}M_{\rm y}-M_{\rm y}\partial_{\rm x}M_{\rm x})$ is the exchange spin current and the chiral spin current is defined as follows $J_{\rm E}=-\frac{\gamma c_{\rm E}E_{\rm y}}{\mu_{\rm 0}M_{\rm s}^{2}}(M_{\rm x}^{2}+M_{\rm y}^{2})$ see Wang2016 . In the region ($x>x_{\rm 0}$) without the electric field, the chiral spin current is zero and the total current is equal to the exchange spin current $J_{\rm s}=J_{\rm ex}$. In the absence of the electric field ($E_{\rm y}=0$), the uniform temperature cannot induce a net magnonic spin current. When the electric field $E_{\rm y}=0.15$ MV/cm is applied in the region $x<x_{\rm 0}$, the density of thermal magnons is obviously enhanced at the interface $x=x_{\rm 0}$. Meanwhile, non-equilibrium magnons diffuse away from the $x=x_{\rm 0}$ and generate the negative (positive) magnonic spin current $J_{\rm s}$ in the right (left) side. $J_{\rm s}$ decreases with $x$ due to the attenuation and boundary reflection. Reversal of the electric field ($E_{\rm y}=-0.15$ MV/cm) reduces the thermal magnons at $x=x_{\rm 0}$, and $J_{\rm s}$ in the right (left) side becomes positive (negative). Besides, the effect of the enhancement of the thermal magnons ($E_{\rm y}>0$) is obviously stronger than the reduction effect ($E_{\rm y}<0$), as shown in Fig. 10(a). The variation of the magnon density $\left|\rho-\rho_{\rm 0}\right|$ and the magnonic spin current $\left|J_{\rm s}\right|$ both increases with the increase of $\left|E_{\rm y}\right|$. Here, $\rho_{\rm 0}$ represents the magnon density without electric field ($E_{\rm y}=0$). Moreover, an increase in the uniform temperature $T$ enhances the variation of the magnon density and the magnonic spin current $J_{\rm s}$, as shown in Fig. 10(b).

VI Skyrmion motion driven by inhomogeneous electric field

To contrast our inhomogeneous E-field-based method with standard recipes of skyrmion motion based either on the electric spin torque current or on thermally assisted magnonic current, we embed the skyrmion at the right side of the interface $x=x_{\rm 0}=-35$ nm and let the induced thermal magnonic spin current $J_{\rm s}$ drives it (as demonstrated by Fig. 11). As is shown in Figs. 12 and 13, the skyrmion moves towards the interface $x=x_{\rm 0}$ along the $x$ axis when the thermal magnon density is enhanced and the exchange magnon spin current is negative $J_{\rm ex}<0$ for $E_{\rm y}>0$. Besides, the skyrmion drifts along the $y$ axis due to the skyrmion Hall effect. Switching sign of the electric field ($E_{\rm y}<0)$ leads to an inversion of the skyrmion motion in both $x$ and $y$ directions (Fig. 14(a)), however with smaller velocities. Besides, the motion at the higher temperature is faster (Fig. 14(b)), i.e. thermal effect enhances the skyrmion velocity.

The STT with spin polarization $\mathrm{p}$ flows along the $y$ axis $c_{\rm j}>0$ and drags the skyrmion in both $x$ and $y$ directions simultaneously at zero temperature $T=0$. The corresponding velocities $v_{\rm x}$ and $v_{\rm y}$ satisfy the condition $v_{\rm x}=\alpha v_{\rm y}$ yzhou2016 . Due to the equivalence of IET and the STT, an interesting question is whether IET as well can drag skyrmion at zero temperature, without thermal assistance and for the same electric field. The answer to this question is positive, as inferred from Fig. 15.

Fig. 15 illustrates the motion of the skyrmion driven by an electric field with a uniform gradient \mbox{\boldmath\mathrm{E}}=E_{\rm gr}x\mbox{\boldmath\mathrm{e}}_{\rm z} MV/m, \mbox{\boldmath\mathrm{p}}_{\rm E}=-\mbox{\boldmath\mathrm{e}}_{\rm y} at zero temperature $T=0$, for $E_{\rm gr}=4.7\times 10^{14}>0$. The skyrmion moves in the $+y$ direction with the velocity $v_{\rm y}=4.2$ m/s. Due to the skyrmion Hall effect, the skyrmion is drifted along the $+x$, and the velocity is $v_{\rm x}=0.005\mathrm{m/s}\approx\alpha v_{\rm y}$ yzhou2016 . Reversing $E_{\rm gr}$ leads to an inversion of the skyrmion motion. Besides, the skyrmion speed increases linearly with $E_{\rm gr}=\partial_{\rm x}E_{\rm z}$, as demonstrated by Fig. 16. In this scenario the spatially homogeneous electric field cannot move the skyrmion.

The important issue is the stability of the skyrmion when passing the interface, i.e. the region of the large electric field gradient. The strong inhomogeneity of the electric field may have an impact on the stability of the skyrmion. The calculations in Fig. 17 demonstrate the stability of the skyrmions for the considered electric field parameters when traversing the interface.

VII Time-dependent electric field

The energy balance of the magnonic spin current depends on two main factors: Due to the phenomenological Gilbert damping, the driven magnonic system continuously loses the energy. On the other hand the ferroelectric polarization \mbox{\boldmath\mathrm{P}}=c_{\rm E}[(\mbox{\boldmath\mathrm{m}}\cdot\nabla)\mbox{\boldmath\mathrm{m}}-\mbox{\boldmath\mathrm{m}}(\nabla\cdot\mbox{\boldmath\mathrm{m}})] and the ME coupling with the external electric field E_{\rm ele}=-\mbox{\boldmath\mathrm{E}}\cdot\mbox{\boldmath\mathrm{P}} supplies a ferroelectric energy to the magnonic subsystem (eventually through the energy needed to generate and maintain $\mathrm{E}$). Thus, the Gilbert damping plays the role of a sink, and the ferroelectric polarization plays the role of the energy reservoir that sustains the spin current. To prove this scenario, instead of the constant electric field, we applied a series of rectangular pulses of the electric field, i.e., the electric field is periodically switched on and off. As we see, when the electric field is switched off, the skyrmion stops moving Fig. 18. The same behavior we observe for the magnonic spin current (not shown), the current disappears soon after the electric field is switched off. After the electric field is switched off, due to the Gilbert damping, the ferroelectric energy becomes zero. Switching on the electric field pumps the ferroelectric energy into the system and restores the ferroelectric torque. Thus, an electric field has to continuously supply the ferroelectric energy to the magnetic system in both cases, either for restoring the ferroelectric term or for compensating energy losses due to damping. For exploring the microscopic mechanisms of the energy exchange we utilize the quantum Hamiltonian and apply the external electric field E=\big{(}0,E_{\rm y},E_{\rm z}\big{)}. We find

[TABLE]

Here, the effective DMI constants are defined as follows $D_{\rm z}=E_{\rm z}Jea/2E_{\rm so}$, D_{\rm y}\big{(}x\big{)}=E_{\rm y}\big{(}x\big{)}Jea/2E_{\rm so} and we assumed that the field applied along the $\mathrm{\hat{y}}$ axis is inhomogeneous in the $\mathrm{\hat{x}}$ direction. Utilizing the Holstein-Primakoff transformation we derive the equation that quantifies the energy exchange between the magnonic $\hat{H}_{\rm mag}=2SJ\sum_{\rm n}\hat{a}^{{\dagger}}_{\rm n}\hat{a}_{\rm n}$ and the ferroelectric subsystems \hat{H}_{\rm D}=\frac{2S}{i}\sum_{\rm n}D_{\rm z,n}\big{(}\hat{a}_{\rm n}\hat{a}^{{\dagger}}_{\rm n+1}-\hat{a}_{\rm n}^{{\dagger}}\hat{a}_{\rm n+1}\big{)}+2S\sum_{\rm n}D_{\rm y,n}\big{(}\hat{a}^{{\dagger}}_{\rm n}\hat{a}_{\rm n}\big{(}\hat{a}^{{\dagger}}_{\rm n+1}+\hat{a}_{\rm n+1}\big{)}-\hat{a}^{{\dagger}}_{\rm n+1}\hat{a}_{\rm n+1}\big{(}\hat{a}^{{\dagger}}_{\rm n}+\hat{a}_{\rm n}\big{)}\big{)}, (here $J$ is the exchange constant, $S$ is the spin of the electron and $\hat{a}^{{\dagger}}_{\rm n},~{}~{}\hat{a}_{\rm n}$ are the magnon creation and annihilation operators). After standard algebraic transformations we obtain

[TABLE]

On the right-hand side of the Eq.(VII) we have the product of two operators. The expectation value of the first operator \langle i\big{(}\hat{a}^{{\dagger}}_{\rm n}-\hat{a}_{\rm n}\big{)}\rangle is nonzero at nonzero temperature, while the expectation value of the second operator \langle\big{(}D_{\rm y,n}\hat{a}^{{\dagger}}_{\rm n+1}\hat{a}_{\rm n+1}-D_{\rm y,n-1}\hat{a}^{{\dagger}}_{\rm n-1}\hat{a}_{\rm n-1}\big{)}\big{\rangle} is nonzero if the translation symmetry in the system is broken, meaning when the inhomogeneous electric field is applied $n\mapsto n+2,~{}~{}\hat{a}^{{\dagger}}_{\rm n-1}~{}\hat{a}_{\rm n-1}\mapsto\hat{a}^{{\dagger}}_{\rm n+1}\hat{a}_{\rm n+1}$ but $D_{\rm y,n+1}\neq D_{\rm y,n}$. Thus, in order to have energy flow from the ferroelectric subsystem to the magnonic subsystem and to sustain the magnonic spin current one needs simultaneously two ingredients: applied thermal bias and a nonuniform electric field. We note that in the standard spin Seebeck experiments, translation invariance is broken by the spatially non-uniform temperature profile, while in our case temperature profile is uniform and translation symmetry is broken by the nonuniform electric field.

VIII Conclusions

Our motivation in the present work has been to explore new ways of controlling the motion of skyrmions in a thin ferromagnetic insulator film. We find that a spatially inhomogeneous electric field serves this purpose in that it leads to a specific torque capable of dragging the skyrmion even at a zero temperature, while thermal effects assist the skyrmion drag and enhance skyrmion velocity. Electric fields are advantages in several ways. They can be generated and temporally controlled in a versatile manner. The spatial inhomogeneities can also be well designed, for instance by nanopatterning the sample by a metallic shielding that blocks the electric field on a microscopic scale. Generally, the current findings point to an interesting and rich spin-current physics of skyrmionic systems driven by local THz fields.

IX Acknowledgment

We acknowledge financial support from DFG through SFB 762 and SFB TRR227, the National Natural Science Foundation of China (No. 11704415, 11674400, 11374373), and the Natural Science Foundation of Hunan Province of China (No. 2018JJ3629, 13JJ2004).