Skyrmion Mass from Spin-Phonon Interaction

Magnetic skyrmions are swirls of magnetization in thin films. They proliferated into condensed matter physics BelPolJETP75 ; WriMerRMR89 ; SonKarKivPRB93 ; StonePRB93 ; Bogdanov94 ; YeKimPRL99 ; AlkStoNat01 from field models of atomic nuclei and topologically stable elementary particles SkyrmePRC58 ; Polyakov-book ; Manton-book ; D1 ; D2 . In ferro- and antiferromagnets they are topological defects of the uniform magnetization (Néel vector) that cannot be easily destroyed. Unlike micron-size magnetic bubbles studied in the past MS-bubbles ; ODell , skyrmions can be small compared to the domain wall thickness, making them promising candidates for topologically-protected nanoscale information processing Nagaosa2013 ; Zhang2015 ; Klaui2016 ; Leonov-NJP2016 ; Hoffmann-PhysRep2017 ; Fert-Nature2017 .

Skyrmions can be moved by current-induced spin-orbit torques Yu-NanoLet2016 ; Fert-Nature2017 ; Legrand-Nanolet2017 ; review2020 . The speed of the information processing with skyrmions depends on their inertia. The effort to compute and measure skyrmion inertial mass has been limited so far. The mass of a skyrmion bubble of dipolar origin, similar to the D$\ddot{\text{o}}$ring mass Doring of the domain wall in the thin-wall approximation, has been discussed by Makhfudz et al. Makhfudz-PRL2012 . Large inertia has been reported in experiments on skyrmion breathing modes in the gigahertz frequency range Buttner-Nature2015 . Similar effects have been observed in the breathing and hypocycloid motion of skyrmions by Shiino et al. Shiino2017 . Inertial mass of electromagnetic origin due to excitation of magnons by a moving skyrmion has been studied by Lin Lin-PRB2017 . Psaroudaki et al. Psaroudaki-PRX2017 have demonstrated that translational symmetry makes classical skyrmions massless. They computed the mass arising from defects, non-uniformity of the magnetic field, and confining potentials, and elucidated contribution of thermal and quantum fluctuations to the mass. Kravchuk et al. Kravchuk-PRB2018 used collective coordinates to demonstrate that skyrmion dynamics in a continuous spin-field model is massless even if one accounts for magnon excitations. Massive skyrmions have been reported by Li et al. Li-PRB2018 in simulations of collective magnetic dynamics on a two-dimensional honeycomb lattice. Measurement of the mass of the oscillating skyrmion in a confined geometry of a semicircular nano-ring has been recently proposed by Liu and Liang Liu-MMM2020 .

The range for the skyrmion mass obtained in experiments is rather broad. It is not always clear whether it is associated with the confined geometry or more fundumental effects studied by theorists. The latter hint towards zero skyrmion mass in the presence of full translational invariance. Crystal lattice violates such invariance. In this Letter we show that skyrmions acquire a finite mass due to the spin-phonon interaction even within translationally invariant continuous spin-field and elastic theories. The physics behind the contribution of the atomic lattice to the skyrmion mass is transparent. The time-dependent spin field corresponding to the moving skyrmion induces, through the magnetoelastic coupling, the motion of the atoms whose inertia contributes to the mass of the skyrmion. Materials that host skyrmions are rather complex. In addition to the dominant exchange interaction, they possess various other kinds of magnetic interactions that are important for stabilization of skyrmions, such as Dzyaloshinskii-Moriya, Zeeman, and crystal field interactions Bogdanov94 ; Bogdanov-Nature2006 ; Heinze-Nature2011 ; Boulle-NatNano2016 ; Leonov-NJP2016 . Pertinent to our purpose, we shall consider in this Letter a toy model of Belavin-Polyakov skyrmion BelPolJETP75 interacting with isotropic elastic environment. We shall study two simple forms of the magnetoelastic coupling, the extreme anisotropic and fully isotropic.

The skyrmion is modelled by the dimensionless three-component spin field ${\bf S}$ of unit length given by BelPolJETP75

	$\displaystyle S_{x}$	$\displaystyle=$	$\displaystyle\frac{2\lambda(x\cos\gamma-y\sin\gamma)}{\lambda^{2}+x^{2}+y^{2}},\,\,S_{y}=\frac{2\lambda(x\sin\gamma+y\cos\gamma)}{\lambda^{2}+x^{2}+y^{2}},$
	$\displaystyle S_{z}$	$\displaystyle=$	$\displaystyle\frac{\lambda^{2}-x^{2}-y^{2}}{\lambda^{2}+x^{2}+y^{2}},\quad{\bf S}^{2}=S_{x}^{2}+S_{y}^{2}+S_{z}^{2}=1,$		(1)

confined to the $xy$ layer of thickness $d$, with ${\bf S}$ looking down at infinity. Here $\lambda$ can be viewed as the lateral size of the skyrmion and $\gamma$ describes the rotation of the spin field, with $\gamma=0,\pi$ and $\gamma=\pm\pi/2$ corresponding to the Néel-type and Bloch-type skyrmions respectively, see, e.g., Ref. DCG-PRB2018, . Magnetoelastic interaction would generally be of the form $A_{ikjl}u_{ik}S_{j}S_{l}$, where

$$u_{ik}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial r_{k}}+\frac{\partial u_{k}}{\partial r_{i}}\right)$$

(2)

is the strain tensor, ${\bf u}$ is the phonon displacement field, and tensor $A_{ikjl}$ represents components of the magnetoelastic energy density.

We shall start with the extreme anisotropic form of the magnetoelastic coupling, $Au_{zz}S_{z}^{2}$, that together with the elastic contribution yields for the energy

$$E=\int d^{3}r\left[\frac{1}{2}\rho\left(\frac{\partial{\bf u}}{\partial t}\right)^{2}+\mu\left(u_{ik}^{2}+\frac{\sigma}{1-2\sigma}u_{ll}^{2}\right)+Au_{zz}S_{z}^{2}\right].$$

(3)

To simplify the problem we assume that the magnetic layer is confined between two non-magnetic semi-infinite solids (see Fig. 1) having the same mass density $\rho$, the same shear modulus $\mu>0$ and the same Poisson coefficient $\sigma=E/(2\mu)-1$ (satisfying $-1\leq\sigma\leq 1/2$), with $E$ being the Young modulus Elasticity . If the speed of the moving skyrmion is small compared to the speed of sound, the elastic deformation adiabatically follows the skyrmion via extremal equation for the energy:

$${\bm{\nabla}}^{2}{\bf u}+\frac{1}{1-2\sigma}{\bm{\nabla}}({\bm{\nabla}}\cdot{\bf u})=-\frac{A}{\mu}\frac{\partial S_{z}^{2}}{\partial z}{\bf e}_{z},$$

(4)

where ${\bf e}_{z}$ is the unit vector along the $z$-axis. Its solution is

$$u_{i}({\bf r})=-\frac{A}{\mu}\int d^{3}r^{\prime}G_{iz}({\bf r}-{\bf r}^{\prime})(\partial S_{z}^{2}/\partial z),$$

(5)

where

$$G_{ik}=\frac{1}{4\pi}\left[\frac{\delta_{ik}}{r}-\frac{1}{4(1-\sigma)}\frac{\partial^{2}r}{\partial r_{i}\partial r_{k}}\right]$$

(6)

is the Green function Elasticity of Eq. (4).

If the skyrmion moves along the $x$-axis at a speed $v$ the solution (5) must be replaced with ${\bf u}(x-vt,y,z)$. Its substitution into the first term of Eq. (3) gives for the kinetic energy, $K.E.=\frac{1}{2}M_{S}v^{2}$, where

$$M_{S}=\rho\int d^{3}r\left(\frac{\partial{\bf u}}{\partial x}\right)^{2}$$

(7)

is the mass of the skyrmion due to spin-phonon coupling. Substitution of Eq. (5) in the expression for the mass yields

$$M_{S}=\rho\left(\frac{A}{\mu}\right)^{2}\int d^{3}r^{\prime}\int d^{3}r^{\prime\prime}F_{zz}({\bf r}^{\prime}-{\bf r^{\prime\prime}})\frac{\partial S_{z}^{2}({\bf r}^{\prime})}{\partial z^{\prime}}\frac{\partial S_{z}^{2}({\bf r}^{\prime\prime})}{\partial z^{\prime\prime}},$$

(8)

where

$$F_{zz}({\bf r}^{\prime}-{\bf r^{\prime\prime}})=\frac{\partial}{\partial x^{\prime}}\frac{\partial}{\partial x^{\prime\prime}}\int d^{3}rG_{iz}({\bf r}-{\bf r}^{\prime})G_{iz}({\bf r}-{\bf r}^{\prime\prime}).$$

(9)

Using the Fourier transform of the Green function (6),

$$G_{ik}(k)=\frac{1}{k^{2}}\left[\delta_{ik}-\frac{1}{2(1-\sigma)}\frac{k_{i}k_{k}}{k^{2}}\right],$$

(10)

$F_{zz}$ can be written as

$$F_{kj}({\bf r}^{\prime}-{\bf r^{\prime\prime}})=\int\frac{d^{3}k}{(2\pi)^{3}}e^{-i{\bf k}({\bf r}^{\prime}-{\bf r}^{\prime\prime})}\frac{k_{x}^{2}}{k^{4}}\left[\delta_{kj}-p\frac{k_{k}k_{j}}{k^{2}}\right],$$

(11)

where

$$p=\frac{1}{(1-\sigma)}\left[1-\frac{1}{4(1-\sigma)}\right]$$

(12)

is the parameter of the elastic theory satisfying $7/16\leq p\leq 1$.

At this point it suffices to consider a thin-film approximation, $d\ll\lambda$, when one can write

$$S_{z}^{2}({\bf r})=S_{z}^{2}({\bm{\rho}})d\delta(z)$$

(13)

with ${\bm{\rho}}=x{\bf e}_{x}+y{\bf e}_{y}$ being the radius-vector in the $xy$ plane of the magnetic layer. Integrating by parts in Eq. (8) one obtains

$$M_{S}=\rho\left(\frac{A}{\mu}\right)^{2}d^{2}\int d^{2}\rho^{\prime}\int d^{2}\rho^{\prime\prime}K({\bm{\rho}}^{\prime}-{\bm{\rho}}^{\prime\prime})S_{z}^{2}({\bm{\rho}}^{\prime})S_{z}^{2}({\bm{\rho}}^{\prime\prime}),$$

(14)

with

$$K({\bm{\rho}})=\int\frac{d^{3}k}{(2\pi)^{3}}\frac{k_{x}^{2}k_{z}^{2}}{k^{4}}\left[1-p\frac{k_{z}^{2}}{k^{2}}\right]\exp(-ik_{x}x-ik_{y}y).$$

(15)

Changing the order of integration in Eq. (14) one can write it as

$$M_{S}=\rho\left(\frac{A}{\mu}\right)^{2}d^{2}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{k_{x}^{2}k_{z}^{2}}{k^{4}}\left[1-p\frac{k_{z}^{2}}{k^{2}}\right]f(k_{\perp}),$$

(16)

where $k_{\perp}=\sqrt{k_{x}^{2}+k_{y}^{2}}$ and

$$f(k_{\perp})=\left|\int d^{2}\rho\,S_{z}^{2}(\rho)\exp(i{\bf k}_{\perp}\cdot{\bm{\rho}})\right|^{2}.$$

(17)

Its independence on $k_{z}$ allows one to integrate over $k_{z}$ in Eq. (16). This leads to

$$M_{S}=\frac{1}{16\pi}\rho\left(\frac{A}{\mu}\right)^{2}d^{2}\left(1-\frac{3}{4}p\right)\int_{0}^{\infty}dk_{\perp}k_{\perp}^{2}f(k_{\perp}).$$

(18)

We now have to compute $f(k_{\perp})$. Substituting $S_{z}(\rho)=(\lambda^{2}-\rho^{2})/(\lambda^{2}+\rho^{2})$ from Eq. (1) into Eq. (16), we get

$$f(k_{\perp})=(8\pi)^{2}{\lambda^{4}}u^{2}(k_{\perp}\lambda),\quad u(q)=-\int_{0}^{\infty}dr\,r\frac{r^{3}J_{0}(qr)}{(1+r^{2})^{2}},$$

(19)

where $J_{0}$ is the Bessel function. This gives for the mass

$$M_{S}=4\pi\rho\left(\frac{A}{\mu}\right)^{2}\left(1-\frac{3}{4}p\right)d^{2}\lambda\int_{0}^{\infty}dq\,q^{2}u^{2}(q).$$

(20)

Here $u(q)$ given by Eq. (19) can be expressed via special functions, which facilitates numerical computation of the integral. The answer yields

$$M_{S}=0.787\left(1-\frac{3}{4}p\right)\left(\frac{A}{\mu}\right)^{2}\rho\,d^{2}\lambda.$$

(21)

We have double-checked this result by performing a more tedious integration in real space without replacing the layer of thickness $d$ with a $\delta$-function. It produces the same answer with the numeric factor given by

$$c=\frac{1}{16\pi}\int d^{2}\bar{\rho}^{\prime}\int d^{2}\bar{\rho}\left[S_{z}^{2}({\bar{\bm{\rho}}}^{\prime}+{\bar{\bm{\rho}}})-S_{z}^{2}({\bar{\bm{\rho}}}^{\prime})\right]^{2}\frac{2\bar{x}^{2}-\bar{y}^{2}}{\bar{\rho}^{5}},$$

(22)

where $\bar{\bm{\rho}}={\bm{\rho}}/\lambda$. This four-dimensional integral reduces to a one-dimensional integral of an awkward elementary function that has a numerical value of $0.785$, very close to the factor in Eq. (21). Notice that for the anisotropic magnetoelastic interaction that we have chosen the mass does not depend on the chirality angle $\gamma$.

To see how general this result is, consider now isotropic magnetoelastic coupling of the form $Au_{ik}S_{i}S_{k}$. Repeating the steps of the previous calculation we obtain for the skyrmion mass

$$M_{S}=\rho\left(\frac{A}{\mu}\right)^{2}d^{2}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{k_{x}^{2}}{k^{4}}\left[|\mathbf{G}|^{2}-p\frac{|\mathbf{k}\cdot\mathbf{G}|^{2}}{k^{2}}\right],$$

(23)

where

$$\mathbf{G}=\int dxdy\left(\mathbf{k}\cdot\mathbf{S}\right)\mathbf{S}\exp(-ik_{x}x-ik_{y}y).$$

(24)

This calculation requires more effort as it involves all three components of the skyrmion spin-field. The final answer reads

$$M_{s}=c(p,\gamma)\rho\left(\frac{A}{\mu}\right)^{2}d^{2}\lambda,$$

(25)

with the numerical factor given by,

	$\displaystyle c(p,\gamma)$	$\displaystyle=$	$\displaystyle 4.118+0.727\cos(2\gamma)$
		$\displaystyle-$	$\displaystyle p[1.612+0.795\cos(2\gamma)+0.255\cos(4\gamma)].$

Thus, in general, one should expect the skyrmion mass to depend on both, the elastic properties of the crystal and the chirality of the skyrmion.

The majority of the materials have $p\sim 1$. At $p=1$ one obtains from Eq. (Skyrmion Mass from Spin-Phonon Interaction) $c=2.186$ for the Néel skyrmion ($\gamma=0$) and $c=2.316$ for the Bloch skyrmion ($\gamma=\pi/2$). Notice that this factor for the anisotropic magnetoelastic coupling at $p=1$ is $0.197$, which is an order of magnitude smaller. Up to that factor the proportionality of the skyrmion mass to $\rho\left({A}/{\mu}\right)^{2}d^{2}\lambda$ is robust. The model correctly captures universal scaling of the mass as the square of the strength of the magnetoelastic coupling, the square of the thickness of the ferromagnetic layer, and the first power of the lateral size of the skyrmion. Notice that the proportionality of the skyrmion phonon mass to its size instead of its volume ($V\sim d\lambda^{2}$) is related to the fact that only spin-field derivatives contribute to the effect. If a thin-wall skyrmion bubble (or a cylindrical domain) of radius $R$ were considered instead, the mass would have been proportional to the area of the wall and would scale linearly with $R$.

To estimate the magnitude of the effect, notice that the magnetoelastic energy density $A$ is of the relativistic origin (it often comprises a noticeable part of the magnetocrystalline anisotropy), while the shear modulus $\mu$ is of the electrostatic origin arising from the coupling between atoms in a crystal. This allows one to roughly estimate the ratio $A/\mu$ to be in the ballpark of $10^{-4}$. At $\rho\sim 5\times 10^{3}$kg/m³ and $d\sim 2$nm it gives $M_{S}$ of the order of a few electron masses for a skyrmion of size $\lambda\sim 10$nm. However, in materials with high magnetostriction this mass can be significantly greater as it scales as square of the strength of the magnetoelastic coupling. The framework used in this Letter for computing the mass of a skyrmion allows one to develop of a software package for obtaining the skyrmion mass in materials with arbitrary crystal symmetry and arbitrary structure of the magneto-elastic coupling.

This work has been supported by the grant No. DE-FG02-93ER45487 funded by the U.S. Department of Energy, Office of Science.