Quantum Dynamics of Skyrmions in Chiral Magnets

We study the quantum propagation of a Skyrmion in chiral magnetic insulators by generalizing the micromagnetic equations of motion to a finite-temperature path integral formalism, using field theoretic tools. Promoting the center of the Skyrmion to a dynamic quantity, the fluctuations around the Skyrmionic configuration give rise to a time-dependent damping of the Skyrmion motion. From the frequency dependence of the damping kernel, we are able to identify the Skyrmion mass, thus providing a microscopic description of the kinematic properties of Skyrmions. When defects are present or a magnetic trap is applied, the Skyrmion mass acquires a finite value proportional to the effective spin, even at vanishingly small temperature. We demonstrate that a Skyrmion in a confined geometry provided by a magnetic trap behaves as a massive particle owing to its quasi-one-dimensional confinement. An additional quantum mass term is predicted, independent of the effective spin, with an explicit temperature dependence which remains finite even at zero temperature.


I. INTRODUCTION
Although skyrmions have been proposed [1,2] long ago, there has been a strong rise in interest in recent years spurred by experimental observations of skyrmionic phases in various magnetic thin films [3][4][5]. Magnetic skyrmions are characterized by a topologically nontrivial mapping from a twodimensional (2D) magnetic system in real space into threedimensional spin space. This mapping leads to a topological charge Q which characterizes the skyrmion and is given by where m is the normalized magnetization vector field and x and y are the spatial coordinates of the 2D magnetic layer [6,7]. Magnetic skyrmions are attractive candidates for magnetic storage of classical information [8,9] because they are topologically stable in the sense that no continuous local deformation in the magnetic texture can remove a skyrmion (i.e., change Q), and they can be manipulated at high speed with relatively low current densities [8,10,11]. Furthermore, because skyrmions coupled to conventional superconductors are known to support Majorana fermions [12], they may be a route to topological quantum computation [13]. The dynamics of any classical magnetic texture with fixed magnitude, m, is governed by the Landau-Lifshitz-Gilbert (LLG) equation [15,16]. The dynamic properties of skyrmions can be considerably simplified by considering the motion of the average magnetic texture [17], which reduces to an equation of motion for the skyrmion's center-of-mass coordinate, known as Thiele's equation. This skyrmion coordinate behaves as a massless particle under a Magnus force proportional to Q [18], and a possible damping is parameterized by a phenomenological velocity-dependent term that is induced by the coupling of the skyrmion motion to other degrees of freedom in the system such as electrons, phonons, magnons, etc.; the microscopic details of this coupling are typically left unspecified. Using this phenomenological approach, several generalized equations have been proposed [19][20][21][22] for a rigid magnetic structure in the presence of external static as well as rapidly changing forces. The thermal diffusion of a driven skyrmion has been studied within the classical framework by a numerical solution of the stochastic LLG equation [19], where random fields are included to model the effects of thermal fluctuations.
However, for small skyrmions at the nanometer scale, the size of skyrmions in state-of-the-art experiments [23][24][25][26], quantum fluctuations are expected to become relevant [27] and the classical LLG equation is expected to break down. In a fully quantum mechanical treatment starting from the spin and other degrees of freedom present in the system the quantum dynamics of the skyrmion should emerge naturally, giving rise for instance to possible mass and dissipation terms linked to the underlying microscopics.
The goal of this work is to provide such a quantum treatment of skyrmions in insulating magnetic films. Thereby we restrict ourselves to systems where only spin degrees of freedom are present and other sources such as phonons and itinerant electrons are assumed to be negligibly small (due to low enough temperatures) or entirely absent, respectively. By employing quantum field theory methods and the Faddeev-Popov techniques for collective coordinates [28][29][30][31][32], applied to magnetic systems [33,34], we derive the quantum dynamics of the skyrmion in a systematic way, thereby going beyond the classical limit.
Within this formalism the skyrmion position, which parametrizes the classical static solution of the Euler-Lagrange equations, is promoted to a time-dependent dynamical variable and a finite perturbation theory in terms of fluctuations around the skyrmion configuration is performed. The interaction of the skyrmion with these magnon modes gives rise to dissipation and to a nonlocal damping kernel whose form we obtain in a closed form and can evaluate analytically in some regimes of physical interest. In particular, we are able to microscopically calculate the contribution to the damping in the Thiele equation from the incoherent magnon modes which in general depends on the magnetic dispersion relation and thus on the details of the free energy of the system.
For asymptotic (imaginary) times, we find that the damping kernel becomes local in time and reduces to a simple mass term which we can calculate explicitly in the presence of nonuniform magnetic traps and spatially dependent exchange interactions for arbitrary temperatures. Specifically, we find that local defects and magnetic traps which break translational symmetry induce a finite mass, even at vanishingly small temperature. As a striking consequence, we find that in quasi-one-dimensional wire geometries the skyrmion acquires a finite mass. Such geometries are provided for instance by anisotropic magnetic traps, similar to linear tracks used for magnetic memory and logic devices [8,35,36]. Surprisingly, this shows that skyrmions in quantum confined geometries exhibit a fundamentally different dynamical behavior compared to the one in unconfined geometries.
The paper is organized as follows. In Sec. II we derive a coherent-state path integral for a generic magnetic texture, while in Sec. III, we utilize this functional integral approach to calculate the damping experienced by such a magnetic texture as a result of the coupling to magnon fluctuations, and we also demonstrate that this damping reduces to a simple mass term. In Sec. IV, we apply our formalism to a chiral magnet and show in Sec.V, that only a spatially nonuniform perturbation to the magnetic texture induces a skyrmion mass at zero temperature. To complete the description we examine in Sec.VI the skyrmion mass at finite temperatures. We conclude with a discussion and some general remarks in Secs. VII and Sec. VIII. Technical details are deferred to the Appendices.

II. COLLECTIVE COORDINATE QUANTIZATION
We consider a thin magnetic insulator with normalized magnetization m(r) = [sin Θ(r) cos Φ(r), sin Θ(r) sin Φ(r, cos Θ(r)] in polar coordinates, where the z axis points along the smallest dimension of the magnet (for the connection to spins in a 2D lattice model, see below). For sufficiently thin magnets, i.e., thinner than the exchange length, the magnetization is constant along the z axis andr = (ρ cos φ,ρ sin φ) is a vector in the xy plane. The Lagrangian of such a magnetic system is given by where S is the magnitude of the spin, N A is the number of magnetic layers along the z axis [37], α is the lattice spacing,Φ is the time derivative of Φ(r, t) and Π(r, t) ≡ cos Θ is canonically conjugate to Φ. The first term in Eq. (2), i.e., the Wess-Zumino or Berry phase contribution (using the north pole parametrization [33]), describes the dynamics of the magnetization while the second term is the free energy which, in this and the following section, we keep general. For reasons of simplicity we assume that the energy density functional F(Φ, Π) is translationally invariant whereas the rotational invariance is broken. This is supported by the fact that, for the majority of the known models which can stabilize magnetic skyrmions, U(1) rotational symmetry is broken by either the chiral Dzyaloshinskii-Moriya interaction or by long-range magnetic dipolar interactions. Upon minimizing L, one recovers the Euler-Lagrange or, equivalently, the Landau-Lifshitz (LL) equations, which govern the classical motion of the magnetization m(r, t).
To investigate the quantum effects of such a magnetic texture, we employ a functional integral formulation in which the partition function is given by Z = DΠDΦ e −S E , where S E = β 0 dτ L is the Euclidean action. Here, L is the imaginary time Lagrangian, which is obtained upon replacing time t with imaginary timeτ = it, and whereβ = 1/k BT for the system temperatureT with k B the Boltzmann constant. Throughout this work we assume that = 1. We also introduce dimensionless variables r =r/lα, τ =τ ε Λ , and β =βε Λ . Here l, a dimensionless constant, and ε Λ , a characteristic energy scale, depend on the system under study and are specified later in Sec. IV. With increasing number of layers the effective spin N A S of the texture increases, and thereby m(r, τ ) approaches a macroscopic description of a classical magnetization field. We therefore consider the limit N A 1 to capture the leading semiclassical effects.
Henceforth, we consider a specific class of free energies which furnish a static stable or metastable solution, m 0 (r) = [sin Θ 0 (r) cos Φ 0 (r), sin Θ 0 (r) sin Φ 0 (r), cos Θ 0 (r)], in which m 0 is characterized by a topological number Q and the magnetization profile Θ 0 approaches a constant value at spatial infinity |r| → ∞. These textures, known as Skyrmions, are illustrated in Fig.1 which depicts a two-dimensional magnetic layer hosting a Skyrmion with Q = 1 embedded in a bulk magnetic insulator. We first develop a formalism that is valid for general Skyrmionic solutions, while an explicit solution of a Skyrmion is discussed in detail later.
The center of mass of such solutions m 0 (r) is denoted by the collective coordinate R(τ ), which is energy independent owing to the translational symmetry of the system. The dynamics of the skyrmion is then described by R(τ ) and the magnetization fluctuations that arise when the skyrmion moves. Because of the nonlinear character of the LL equation these fluctuations couple back to R(τ ) and thus affect the skyrmion motion. To describe this effect we perform a canonical transformation of the path integral variables which separates the collective coordinate R(τ ) from the fluctuations around the skyrmion [28,29], In the presence of translational symmetry there is a degenerate pair of zero-frequency modes ξ ∝ ∂ i Φ and η ∝ ∂ i Π, associated with translation of the skyrmion position in the i = x, y direction. Instead of fluctuations ξ and η, it is more convenient to work with the spinors χ and χ † defined by We denote as Y x and Y y the degenerate pair of translational zero modes expressed in the spinor notation of Eq. (4). We now construct a finite perturbation theory in terms of fluctuations around m 0 (r − R), where the system is described from a translated spatial framer = r − R(τ ). To second order in χ, the action S E = S cl + S fl + S I separates into three parts. The first part describes the translational motion of m 0 , whereS = Sl 2 and where we ignore an overall constant from the configuration energy (lα) 2 drF(Φ 0 , Π 0 ). Note that S WZ cl contains a gauge-dependent boundary term proportional to −iSN A dτ drṘ · ∇Φ 0 = −iSN A π dτṘ · C, where C = (1/π) dr ∇Φ 0 is the skyrmion chirality vector. This topological term does not affect the classical equation of motion, but could in general have observable effects, for example, in the tunneling probability [33,40]. The most common class of solutions of magnetic skyrmions has vanishing topological term with C = 0, thus becoming manifestly a gauge invariant theory. The second part consists of the fluctuations χ(r, τ ) around m 0 (r), which we can consider as spin waves or magnons, given by Here we introduce the compact scalar product notation for operators and functions, where the scalar product in spinor space is left implicit. Further we define where δ χ is the functional derivative with respect to the field χ(r, τ ), and G is the magnon Green's function. Finally, the coupling between the center of mass R(τ ) and the spin waves χ(r, τ ) is given by where Here, repeated indices, i, j = x, y, are summed over and we introduce the abbreviations into the partition function, where δ(·) is the Dirac delta function. We then obtain where J G = dG(τ )/dR(τ ) is the Jacobian of the coordinate transformation and is treated as an additional term in the action, adaptable for perturbative calculation. We should stress that the scalar product χ † (r)Y i (r) (in the 2-spinor space) is determined by the properties of the Hamiltonian H. In most cases H is not Hermitian for the standard metric and orthogonality conditions are enforced by an inner product of the form χ † (r)σ z Y i (r) [39]. We next integrate out the fluctuations χ and χ † from the partition function neglecting terms O(1) in N A which originate from the FP determinant det(J G ).
The integral over χ and χ † can be reduced to a Gaussian form by completing the square,χ = χ + (G + K) −1 J . The inverse of an operator, O −1 , is defined by dr dτ O(r, τ ;r , τ )O −1 (r , τ ;r , τ ) = δ(r −r )δ(τ − τ ). After integration, the partition function reduces tõ where the prime on the determinant excludes the zero modes, and where Upon exponentiating the determinant and expanding to lowest nonvanishing order inṘ, we get where Tr O = dτ drO(r, τ ;r, τ ) denotes the trace and where we use that Tr (G −1 K) = 0 due to periodic boundary conditions in time. The approximations in Eqs. (14) and Eq. (15) neglect terms O(Ṙ 3 ).
The quantization of a skyrmionic field presented in this section, with the dynamical part of the action proportional toΦΠ, is significantly different from the canonical quantization of theories with a dynamical part proportional toΦ 2 , for example the quantization of a kink soliton of a (1+1) nonlinear field theory [28] or a real two-component scalar field in (2+1) dimensions [41].

III. SKYRMION MASS AND DISSIPATION
BecauseZ[R] is, in general, nonlocal in time, it naturally encodes dissipation of the skyrmion that is captured by the action Let Ψ s n (r)e iων τ / √ β be the normalized eigenvectors of the operator G, where ω ν = 2πν/β are the Matsubara frequencies, with ν integer, and where ε s n is the corresponding eigenenergy of the magnons, given as solutions of the eigenvalue problem (EVP) HΨ s n = 2ε s n σ z Ψ s n . Here we introduce index s = ±1 to distinguish between particle states Ψ 1 n (r) with eigenfrequency ε 1 n = +ε n and antiparticle states Ψ −1 n (r) with eigenfrequency ε −1 n = −ε n . We refer the reader to Appendix B for a detailed discussion. Using this complete basis we get (18) Insertion of the series expression (18) into Eq. (16) yields where summation over repeated indices i, j = x, y is implied. The damping kernel γ ij = γ 0 ij + γ T ij consists of a term proportional to the effective spin N A S, whereε s n = 2ε s n /S withS = Sl 2 , the prime in the sum denotes omission of the zero modes, and the elements A n,s ij are given by The second part of the kernel is independent of the effective spin N AS and is calculated by evaluating the trace in Eq. (16) in the eigenbasis of the operator G, with matrix elements The sums (20) and (22) over Matsubara frequencies ω ν can be explicitly performed using the following exact relation, where the RHS is understood to be periodically extended beyond |τ | ≤ β/2. The action Eq. (19) takes the typical form of a system coupled to a thermal reservoir governed by a linear dissipative process [42]. The expression for S d contains contributions that are nonlocal in imaginary time and describe effective correlations of R i with itself at different times. These correlations are a consequence of an emission of a virtual magnon at time τ generated by the motion of the skyrmion which is then reabsorbed by the skyrmion at a later time τ . The contributions at equal times, if finite, will give rise to a mass of the skyrmion, which otherwise is absent. The form of the kernel γ ij obviously depends on the microscopic details of the system and determines the amplitude of these effects as see below.
For the case of a translationally invariant free energy, the translational modes with zero energy are proportional to the functions f i . Because the summation over the magnon modes excludes translation, A n,s ij = 0, we conclude that there is no mass contribution from the damping kernel γ 0 ij . This is in agreement with the well-known result that the skyrmion behaves like a massless particle in the presence of a magnetic field [17]. However, the situation changes if we allow for a perturbation in our free energy, denoted by V(r), which explicitly breaks translational symmetry. Treating V(r) in standard perturbation theory, the magnon states become in lowest order, where W = δ χ † δ χ V| χ=χ † =0 . Using the orthogonality properties of the unperturbed magnon modes, we find to lowest order in V(r), which in general is nonzero (see below). In contrast to γ 0 ij , the second contribution to dissipation, γ T ij , assumes a finite value Jε gap =JS 2ε gap in general, i.e., we find nonzero matrix elements B ns,n s ij = 0, even for a translationally invariant system.
The first variation of the action S cl +S d vanishes for the extremal (classical) path R i (τ ) obeying the following equation of motion in bosonic Matsubara ω ν frequency space, where ij is the Levi-Civita tensor and summation over repeated indices i, j is implied. Here, the β-periodic functions in imaginary time are expanded into Fourier series as and we assume the presence of a potential term V which breaks translational symmetry which is incorporated in the equation of motion through The first term in (27) originates from the Berry phase and results in a Magnus force acting on the skyrmion proportional to the winding numberQ = 4πSN A Q, recovering well-known results [17,18]. The Fourier coefficients γ ων ij are found from Eqs. (20) and (22), The quantum dynamics of the skyrmion is described by a kernel γ ων ij with a strong frequency dependence, in contrast to the typical Gilbert damping. However, to be consistent with the approximation that terms O(Ṙ 3 ) are negligible, we note that the damping kernel γ ων ij includes processes occurring at frequencies ω ν γ ων ij and the behavior for higher frequencies is beyond the range of validity of our derivation. For low frequencies, i.e. asymptotic imaginary time τ , the damping kernel γ ων ij in Eq. (27) can be expanded, , which holds for all characteristic frequencies in the range ω ν <ε gap , and consequently for temperatures in the range T <ε gap , whereε gap is the lowest magnon energy gap. We therefore observe that the effect of the damping can be reduced to a mass term defined as where the first term, which is temperature-independent and proportional to the effective spin N AS , is given by while the second mass term is explicitly temperature dependent and independent of N AS , Note that there is a nonvanishing particle-antiparticle contribution to the effective mass at zero temperature, which reflects the presence of quantum fluctuations at zero temperature.
The definition of the mass as the zero-frequency limit of the real part of the damping kernel (γ ων ij + γ −ων ji ) is valid under the assumption that the real-time dynamics of the skyrmion is described by circular modesω that are small compared to the magnon gapω <ε gap . For example, in the presence of a parabolic potential well V = K(R 2 x + R 2 y ), the characteristic frequency of the circular motion in the absence of a mass is ω = K/Q. Here, V is a perturbation to the magnon Hamiltonian H andQ is proportional to the effective spin N A S, thus the iω ν → 0 limit in the mass definition is well justified. In the opposite case of a large dominant frequency, the mass tensor is refined aroundω as ]. Finally, note that ordinary friction, of the form of a damping proportional to the velocity of the skyrmion, is also included in our formulas and can be identified by the frequency-dependent co- The equation of motion Eq. (27) and the mass tensor Eq. (31) are the main results of this section. They are valid for a general magnetic texture m(r, τ ) satisfying the specified requirements, and govern the dynamics of its center-of-mass position R(τ ), taking into account quantum and thermal fluctuations and allowing for a potential V that weakly breaks the translational symmetry. This equation includes only dissipation produced by the magnon modes, while other sources giving rise to damping such as phonons and itinerant electrons are not included. Thus, we also must assume that we work at low enough temperatures such that phonon effects can be safely neglected. It is worth mentioning that additional forces Magnetic field h should be included in Eq. (27) in the presence of a constant magnon current generated by a source [43] or by a magnon current induced by a temperature gradient [44]. The presence of a non-negligible mass term leads to additional oscillatory modes in the real-time dynamics of the skyrmion [20,45]. From the structure of the elements A n ij of Eq. (21) and B n,n ij of Eq. (23), it is easy to see that in the presence of a nonuniform but isotropic potential V or at finite T , the diagonal mass terms are equal, M xx = M yy , while off-diagonal terms vanish, M xy = M yx = 0, because the real part of matrix elements A n,s ij and B ns,n s ij is antisymmetric in indices i, j = x, y. Therefore, in the absence of a driving force and assuming that the confinement is given by a parabolic well V = K(R 2 x + R 2 y )/2, the real-time dynamics of the skyrmion is oscillatory with two characteristic circular modes given bỹ where we introduce M = M xx = M yy . Finally, the most general case of an anisotropic potential V gives rise to four distinct circular modes, originating from the fact that the mass tensor elements are no longer equal, a result that could possibly have interesting experimental implications.

IV. CHIRAL MAGNETS
In this section, we apply the formalism derived in the previous section explicitly to a skyrmionic magnetic texture, focusing on the low-temperature contribution to the mass, in the presence of a perturbation in the free energy that breaks translational symmetry. Skyrmions, known to be stable or metastable solutions in chiral magnets with an applied magnetic field, are described on a square lattice by the spin Hamiltonian, Here, S a = (S x a , S y a , S z a ) is the spin of size S at lattice site a and e = (e x , e y ) is the lattice unit vector. We assume the onsite anisotropy constant to be positive,K > 0, so that there is an easy axis along z and fix the chirality of the Dzyaloshinskii-Moriya (DM) interaction,D. The strength of the magnetic field along the z axis isH and the exchange interaction constant is assumed to be ferromagnetic,J > 0.
Using coherent state representation for the spins (assuming S 1) and passing to the continuum limit [33], we find that Eq. (36) corresponds to the free energy of the form where J =JS 2 , K =KS 2 /α 2 , H = gµ BH S 2 /α 2 , and D =DS 2 /α. We now introduce dimensionless spatial variables r =r/(lα) where l =J/D = J/(Dα), while energies are measured in units of ε Λ =JS 2 = J. The energy density in reduced units is defined as w(Φ, Π) = (lα) 2 εΛ F(Φ, Π) (see Eq. (A1)). We consider the regime of isolated skyrmions, which exist for a wide range of parameters K and H as a stable or metastable state. A skyrmion with Q = −1 is parametrized by (see Appendix A) while the helicity, defined by the sign in Φ 0 , is commensurate with the DM interaction. The structure of the stationary  skyrmion Θ 0 (ρ) is determined by the Euler-Lagrange equation derived from skyrmion energy w(Φ 0 , Π 0 ) and is given in Eq. (A4). Because there is no known analytic solution for Θ 0 (ρ), the following function can be used as an approximate solution, Roughly speaking, ρ 0 = 2/(2κ + h) determines the amount of spins that are noncollinear with the easy axis, while λ, which we obtain numerically using an algorithm based on Runge-Kutta formulas [46], is the skyrmion radius of the area for which the spins are parallel to the easy axis. Here, κ =KJ/D 2 and h = gµ BHJ /D 2 are the dimensionless parameters describing easy axis anisotropy and magnetic field, respectively (see Table. I). In the opposite limit of a largeradius skyrmion the magnetization profile is described by where the parameters λ and ∆ 0 are calculated numerically by fitting the approximate function (40) to the numerical solution of the Euler-Lagrange equation (A4). Magnetization profiles of skyrmions in the small-and large-radius limit obtained numerically are depicted in Fig. 2. A detailed discussion of the magnon eigenstates Ψ n is provided in Appendix B. Apart from the magnon scattering states Ψ sc n that lie above the magnon gap ε MS = κ+h/2, induced by the anisotropy and magnetic field, there exist massive internal modes Ψ bs n that are found for energies 0 < ε n ≤ ε MS . These bound states correspond to deformations of the skyrmion into breathing modes and were discovered numerically in Refs. 43 and 47.

V. INDUCED SKYRMION MASS
The mass tensor M 0 ij (32) of a magnetic skyrmion proportional to the effective spin N AS is now investigated, which is nonzero only in the presence of a perturbation in the free energy that breaks translational symmetry. We should stress that once a perturbation is considered, the skyrmion acquires a mass equal to M 0 ij at any temperature. We consider three isotropic mechanisms for the generation of a skyrmion mass, all breaking translational invariance. (1) a magnetic trap, generated by a nonuniform magnetic field, (2) a local defect that alters the exchange constant, and (3) a periodic variation in the exchange constant. In all three cases, we assume the translationally noninvariant terms to be much smaller than those that are translationally symmetric, | drV(Φ 0 , Π 0 )| S E , and they can thus be treated perturbatively. In all cases where an isotropic potential in x and y direction is applied, we find where λ x and λ y characterize the size of the trap along the x and y directions, respectively. See Fig. 3 for a visualization of the radially symmetric, λ x = λ y , and asymmetric, λ x = λ y , magnetic traps. Because the propagating magnon modes are suppressed, we numerically calculate the mass terms by including only the lowest energy bound states, which are themselves found variationally, in the sum of Eq. (31). In Fig. 4 we plot these results as a function of λ x for a radially symmetric trap, λ y = λ smaller than a critical value λ cr x ∼ λ, the mass increases with increasing λ x , while for λ x λ cr x the mass decreases monotonically.
When there is a local defect in the crystal structure of the chiral magnet, e.g. at the origin, the exchange interaction is a function of position around that point which we model as where J D and λ D parametrize the strength and size of the defect, respectively. In Fig. 5 we plot the skyrmion mass due to such a defect as a function of applied magnetic field for several values of λ D and J D = 0.1. The behavior of Fig. 5 implies that as long as λ D < λ, the mass M increases with increasing (decreasing) magnetic field (Skyrmion size), peaking around λ D ∼ λ and then for λ D > λ decreases with decreasing skyrmion size (see Fig. 5(b) and (c)). We see no peak when λ D = 2, Fig.5 (a), because the skyrmion radius for all positive magnetic fields is smaller than the size of the defect. In addition, the dependence of M on λ D is presented in the inset of Fig. 5, where a peaked behavior is suggested with a maximum around λ. We note that although it has been shown that a model with spatially dependent exchange interactions similar to Eq. (42) acts as a pinning potential [48], i.e., ∂V/∂R i in Eq. 27, our findings show that additionally a finite mass is generated. Finally, because of the discrete underlying lattice the exchange interaction can obtain a periodic modulation along the lattice vectors, where J P and λ P are the strength and period, respectively, by which the exchange deviates from J 0 . In Fig. 6 we plot the magnetic field dependence of M for λ P = 1.5 and J P = 0.1.
Focusing on the dependence of M on the skyrmion radius λ, illustrated in the inset of Fig. 6 (a), it becomes apparent that M grows with the area λ 2 . Further results are provided in the inset of Fig. 6 (b), where we examine the dependence of the mass on the period λ P for h = 0.1, λ = 1.57, and J P = 0.01. Figure 6 (b) reveals a peaked behavior around λ/2.

A. Massive Skyrmion in quasi-1D confinement
The case of a magnetic skyrmion confined in a quasi-onedimensional space can be realized by applying an anisotropic magnetic trap whose size is much larger in one spatial direction than in the other, i.e., λ x λ y . A pictorial representation of a skyrmion in the presence of an anisotropic magnetic trap is provided in Fig. 3. When λ x = λ y , the elements of the mass tensor M 0 ij are no longer equal, but acquire different values depending on the size of the magnetic trap. In Fig. 7 we plot the dependence of M 0 xx and M 0 yy on λ x for a fixed value of λ y = 2 and a skyrmion with κ = 0.7, h = 0.1 and λ = 1.57. We find that M 0 xx is maximized at λ; conversely, M 0 yy grows until λ y ∼ λ x upon which it approaches a constant value.
A corollary of the preceding discussion is that the highly anisotropic magnetic trap acts as a magnetic quantum wire, along which the skyrmion moves as a massive particle due to its one-dimensional confinement. Our investigation demonstrates that skyrmions observed in confined geometries such as magnetic nanowires [50] exhibit fundamentally different dynamical behavior compared to the one in unconfined 2D geometries. It would be interesting to test this prediction.

VI. SKYRMION MASS AT FINITE TEMPERATURE
When no translational symmetry breaking term is present, the only contribution to the mass is given by the term M T ij , which is independent of N AS and finite even at vanishingly small temperatures. The series expression Eq. (33) requires the summation over the full magnon spectrum, that is over scattering states Ψ sc m,k which are classified by the azimuthal quantum number m as well as the radial momentum k ≥ 0 with energy ε k = ε MS + k 2 , with ε MS = κ + h/2. In addition to scattering states, one needs to take into account localized modes Ψ bs m , classified by m with energies in the range 0 < ε m ≤ ε MS . Therefore we consider the following terms, where M sc ij denotes the contribution from scattering states, M bs ij the contribution from bound states, and M bs-sc ij takes into account combinations of bound-scattering states. Here, we focus on large-radius skyrmions for low magnetic fields and easy axis anisotropy, where the skyrmion is expected to support more internal modes with energy below ε MS than for small-radius skyrmions. As the magnetic field or easy-axis anisotropy is increased, the radius decreases and the modes sequentially leave the gap region, pass to the continuous spectrum and transform into quasilocalized modes. Therefore, there is a bounded parameter region of h and κ where the requirement ε m ≤ ε MS is fulfilled.  Fig. 9. Explicit expressions on summations over scattering states are given in Appendix C. Summations over the quantum number m converge rapidly and are bounded between −7 m 7. In Fig. 9, we calculate M T ij for a discrete set of h and T and then we adopt a process to obtain the mass as a smooth function of (h, T ).

VII. DISCUSSION
Equation (31) is the first closed formula for the effective mass of a skyrmion, obtained microscopically in the presence of arbitrary external perturbations arising from defects, nonuniform magnetic fields, and external potentials and at finite temperature. Skyrmion phases have been observed in several metallic ferromagnets, including MnSi [3] and FeGe [51], as well as in the insulator Cu 2 OSeO 3 [52,53]. Typical parameter values for Cu 2 OSeO 3 crystals areJ = 33.36 × 10 −4 eV,D = 7.47 × 10 −4 eV and α = 8.911Å (J/D = 4.46) [54]. The induced mass in physical units is then given by M 0 ij × N A · 4 · 10 −26 kg, which corresponds to 5 × 10 4 electron masses. The scale that separates the quantum from the finite-temperature regime is given by the magnon gapε gap ∼ 2 K. A typical scale for the quantum mass M T ij is ∼ 10 −28 kg with a temperature dependence given in Eq. (33) and plotted in Fig. 9. For high temperatures, k BT >ε gap , our approximation breaks down since the reduction of the damping kernel to an effective mass in Eq. (31) is not valid anymore. We remark that the skyrmion dynamics in this classical regime, ε gap < k BT J, was studied in Ref. [19] via a stochastic LLG equation giving a mass that is independent of temperature. This suggests that the weak temperature dependence of the mass found here in the quantum regime is flattened off by some second-order processes activated at higher temperatures, such as magnon-magnon interactions (neglected in this work). It would be interesting to include such higher-order terms in our microscopic approach systematically. We leave such a study for future work.
Further, we note that the authors of Ref. [19] predict a strong frequency dependence of the dynamics of the skyrmion motion, as a result of the presence of time-dependent forces, as well as the coupling of the skyrmion to magnon excitations. Among the results so obtained is a large numerical value temperature independent effective mass m(ω) as well as a frequency-dependent Gilbert damping D(ω). To make contact with these results, a frequency-dependent mass is identified as m(ω) → 2 Reγ ων ii and Gilbert damping is D(ω) → −2ω ν Imγ ων ii , where γ ων ii due to magnons is known exactly and given in Eq. (30). Since eigenenergies and frequencies scale as 1/λ 2 , where λ is the skyrmion radius, from our explicit expressions we obtain that, m(ω) ∼ λ 2 (see inset of Fig. 6) while D(ω) is independent of λ, a result which is consistent with Ref. [19]. Moreover, a large inertial mass is found for the case of thermal diffusion which induces a large skyrmion distortion (strong translational symmetry breaking term), while an almost vanishing mass was found for the case of an electric current-driven motion, when the skyrmion flows with the current with little distortion (weak translational symmetry breaking term). This result is consistent with the interpretation of skyrmion mass generation provided by our calculation; any term which explicitly breaks translational symmetry will induce a finite mass even at vanishingly small temperature.
Quite generally, the skyrmion dynamics obtained by the LLG equation and the (quantum) damping obtained exactly in closed form in the present work are complementary methods.
The former provides a reliable framework especially when a microscopic description is difficult to obtain, for example the inclusion of a large number of microscopic degrees of freedom such as electrons. However, the latter can provide a detailed microscopic origin of dissipation beyond the classical limit and beyond the predictability of the LLG equation.
In Ref. [55] an inertial term was estimated that appears due to a confining potential, independent of the skyrmion size λ and proportional to the stiffness of the harmonic confining potential. Finally, it is also worth mentioning that an approximate evaluation of the mass term of a skyrmion magnetic bubble, a circular domain wall stabilized by long-range interactions has been found in Ref. [20], derived by generalizing Thiele's approach beyond steady motion. We emphasize that the mass tensor depends on the inverse of the magnon energy ε n . We therefore expect that bound states, which correspond to deformations of the skyrmion with energies below the gap, have a larger contribution. This finding is consistent with Ref. [20] where an approximate value for the mass of a magnetic bubble domain has been calculated by considering contributions from distortions of the circular domain wall. In addition, a large value of inertia mass of a magnetic bubble has been experimentally reported in Ref. [45]. The bubbles were subjected to a magnetic disk and the presence of an inertial mass was attributed to a breathing mode associated with a change of the bubble size.
While several studies demonstrate the association of inertia mass to deformations of the skyrmion, an explicit analytical formula has been missing so far and is provided here. Our detailed calculations shed light on the mechanism of mass generation, and we also provide an analytic formula, valid in the presence of arbitrary external perturbations which can serve as a basis for a variety of further studies.

VIII. CONCLUDING REMARKS
The main task of the paper is to describe dissipation of a magnetic skyrmion using the canonical scheme of quantization. Within this formalism the skyrmion position is promoted to a time-dependent dynamical variable and a finite perturbation theory in terms of fluctuations around the skyrmionic configuration is performed. The interaction of the skyrmion with these magnon modes gives rise to a damping term that is found to be time dependent, in contrast to the typical Gilbert damping. The time-dependence of damping kernel indicates that damping depends on the velocity of the skyrmion of the past.
Perhaps, the most interesting feature is that the effect of damping is reduced to a mass term in some limits. We demonstrate that a massless skyrmion at the classical level is a consequence of the translational symmetry assumed for the system. Our investigation suggests that, even at vanishingly small temperature, the skyrmion mass is nonzero in the presence of an external perturbation arising from defects, nonuniform magnetic fields and external potentials. Among the translational symmetry breaking terms we consider here, we emphasize the case of a quasi-one dimensional wire provided by an anisotropic magnetic trap, along which the skyrmion moves as a massive particle owing to its confinement. To complete the description we also examine damping in the absence of perturbations, and we calculate an effective spin-independent quantum term with an explicit temperature dependence and a nonvanishing particle-antiparticle contribution at zero temperature. This result reflects the fact that not only thermal but also quantum fluctuations contribute to a mass term.
This picture suggests that in the absence of perturbations, massive magnon modes are activated at finite temperatures. A utilization of the improved skyrmion dynamics calculated in the present paper could be towards the direction of investigating the possibility of depinning via quantum tunneling, a behavior that has been studied both theoretically [62] and experimentally [63] in similar mesoscopic magnetic structures, such as magnetic domain walls.
In our present study, we exclusively focus on damping caused by intrinsic mechanisms originating from the interaction of the skyrmion with magnon modes which it generates and couples back to, giving rise to damping. Additional sources of damping could arise if one takes into account timedependent forces or the interaction of the skyrmion to extrinsic degrees of freedom, like a current of electrons or phonons. However, this is beyond the scope of this work, and we leave it as a motivation for further studies making use of the general formalism developed here.

Non-Hermiticity
Magnon modes are solution of the eigenvalue problem HΨ R n,a = 2ε n σ z Ψ R n,a , where the explicit form of the Hermitian Hamiltonian H is given in Eq. (B13) below. Here, the index R denotes the fact that Ψ R n,a are right eigenvectors of the non-Hermitian matrix σ z H. Further, n is the quantum number and a counts for any possible degeneracy. Note that the pseudo-Hermitian operator σ z H satisfies the relation In addition to the right eigenstates Ψ R n,a of operator σ z H, we introduce left eigenstates of the Hermitian adjoint matrix (σ z H) † Ψ L n,a = 2ε * n Ψ L n,a . Hermiticity of the operator σ z H is assured with respect to a new inner product between left and right eigenstates Ψ L n,a |Ψ R m,a = δ n,m [14]. Using Eq.(B1), it can be shown that left and right eigenvectors are connected through the relationship Ψ L n,a | = Ψ R n,a |σ z . Therefore, biorthogonality conditions are enforced as Ψ n,a |σ z |Ψ m,a = δ n,m , where we have now dropped the index R/L which distinguish between right and left eigenfunctions. The identity operator is given by while the trace of an operator A is computed as

Particle-Antiparticle symmetry
In this section we discuss a particular symmetry for the magnon Hamiltonian H which is known for the case of magnetic excitations in ferromagnetic 2D skyrmions [19] and for vortices in 2D easy-plane ferromagnets [39]. In particular, there exist a transformation matrix C = Kσ x , where K denotes complex conjugation, under which Hamiltonian H is invariant CHC † = H, while the Pauli matrix σ z changes sign Cσ z C † = −σ z . This symmetry generates an additional class of solutions of the same eigenvalue problem, but with negative eigenfrequencies.
More specifically, it can be shown that if states Ψ n,a are solutions of the eigenvalue problem (EVP) HΨ n,a = 2ε n σ z Ψ n,a with eigenfrequencies ε n ≥ 0, then states CΨ n,a are solutions of the same EVP (B5) with eigenfrequency −ε n . Such a symmetry is called particle-antiparticle. To describe these two classes of solutions, we introduce the index s = ±1 such that Ψ 1 n with eigenfrequency ε 1 n = +ε n corresponds to particles, while Ψ −1 n,a = CΨ 1 n,a with ε −1 n = −ε n to antiparticles. The biorthogonality conditions for the solutions Ψ s n,a are Ψ s n,a |σ z |Ψ s m,a = sδ n,m Ψ s n,a |σ z |Ψ −s m,a = 0 .
From the conditions (B6), it follows that the energy levels are not necessarily linked to the eigenfrequencies ε n . The measure of the energy of a mode Ψ n,a is given by the quantity E n = Ψ s n,a |H|Ψ s n,a = ε s n Ψ s n,a |σ z |Ψ s n,a .
It is easily demonstrated that modes with eigenenergies ε n and −ε n have their energies E n equal in sign (positive) and value.
The unity operator is redefined as and the trace of an operator changes accordingly,

Eigenvalue problem
Here we consider the spectrum of the magnon modes on the skyrmion background. For reasons of convenience, in this section we use the notation Ψ n to denote particle states, and antiparticles states are recovered using the particle-antiparticle symmetry described in the preceding subsection. Magnon scattering states are obtained for energies ε n ≥ ε MS , with ε MS = κ + h/2. In addition to scattering states, we expect localized modes that correspond to deformations of the skyrmion into polygons (breathing modes) in the range 0 < ε n < ε MS . The existence of such modes was found numerically in Refs. [47] and [43]. An extensive analysis of the magnon spectrum has been provided in Ref. [43], where the authors numerically determine the magnon bound states and provide an analytical and numerical description of magnon scattering states. Here we repeat some steps of the analysis adapted to our puprose and introduce a variational approach for the calculation of the bound states in the small and large radius limit of skyrmions.
To begin with, we seek for solutions of the EVP (B5) , where the Hamiltonian is given by Here, the potential terms V (ρ), W (ρ), and U 0 (ρ) are defined by and where the quantum number n labels both localized modes and scattering states. Next, we represent solutions in terms of wave expansions Ψ n = e imφ ψ n,m (ρ), and the EVP is written as H m ψ n,m (ρ) = ε n,m σ z ψ n,m (ρ) with One of the translational modes with zero energy is given by and the other is found using the particle-antiparticle symmetry Here, the upper index corresponds to particle/antiparticle while the lower to the quantum number m.

a. Small-Radius Skyrmion
For magnetic solitons in 2D easy-axis ferromagnets in the small-radius limit there is a bound state with m = −1 that remains localized and is associated with the soliton displacement [60]. We extend these results to the skyrmion field, and we seek for a localized mode with m = −1 by performing a variational calculation. We choose trial functions Ψ −1 = e −iφ ψ −1 (ρ) as [60] ψ −1 (ρ) = Θ 0 − sin Θ 0 /ρ + aρ 2 Θ 0 − bρ sin Θ 0 Θ 0 − sin Θ 0 /ρ + aρ 2 Θ 0 + bρ sin Θ 0 , (B17) and minimize the functional with respect to the variational parameters a and b in order to calculate the energy ε −1 . Then, the parameters a and b can be found from the normalization conditions (B14). In Fig. 10 (a) we illustrate the energy ε −1 as a function of magnetic field h for κ = 0.7. Within the limits of the variational approach, we find ε −1 < ε MS for a wide range of magnetic fields up to h = 8.

b. Large-Radius Skyrmion
In the limit of large radius, the magnetization profile of the skyrmion is described by Eq. (A6) and the localized modes Ψ bs m with energy ε m are again found variationally, but now using the trial functions [49] Ψ bs m (r) = where f 0 (ρ) = A/ cosh( ρ−λ ∆0 ) and A is chosen such that dρρf 2 0 = 1. The condition of minimizing the energy functional U = dr(Ψ bs m ) † (H m − σ z ε m )Ψ bs m , along with normalization conditions satisfied by the functions Ψ bs m (ρ) will specify the eigenenergies ε m as well as the variational parameters a m and b m . In Fig. 10 (b) we display the energies ε 0 and ε 2 as a function of magnetic field h for κ = 0.1.

c. Scattering States
To complete the description of the magnon spectrum we need to include the scattering states Ψ sc m,k (r) classified by m as well as the radial momentum k 0. Here we repeat some steps of the calculation of Ψ sc m,k (r) originally presented in Ref. [43] for reasons of a complete discussion. In the absence of a skyrmion, magnon scattering states are described by Next, we consider the particle-particle expression for M bs-sc xx , with integrals D m 1 = (f 0 J m ), D m 2 = (f 0 J m /ρ) and D m 3 = (f 0 J m Θ 0 cot Θ 0 ). We use (· · · ) = ∞ 0 · · · dρρ. The normalization constantc m is defined as for a finite system size of area πL 2 . Note that we have neglected terms proportional to sin(δ m ) under the assumption cos(δ m ) sin(δ m ). In Fig. 9-(a) we present the magnetic field dependence of M bs-sc ii for T = 0.05 for a skyrmion with κ = 0.1, S = 1,J/D = 4 and L → ∞. For this choice, we take into account two localized modes with energy less than ε MS in the magnetic field region between h = 0.34, where the skyrmion energy becomes positive, up to h = 0.42 where the bound state energy passes into to the continuous spectrum (see Fig. 10 (b)). To calculate the particle-particle contribution from the scattering states, we use the states Ψ sc m,k given in Eq. (B23), and under the assumption that cos 2 (δ m (k)) sin 2 (δ m (k)) we arrive at the following result, sinh 2 (βε k /2) with g k,k m,m =c m (k) 2c m (k ) 2 [− cos 2 (δ m (k)) cos 2 (δ m (k ))+ 1], C k,k (β) = [coth(βε k /2) − coth(βε k /2)] and G k,k m,m = g k,k m,m k (m+m +1) (k ) −(m+m +3) . In Fig. 9-(b) we plot the temperature and magnetic field dependence of M sc ii for a skyrmion with κ = 0.1, S = 1,J/D = 4 and L → ∞. Summations over quantum number m converge and are bounded between −7 m 7. The particle-particle and antiparticle-antiparticle contribution to the effective mass are equal, since operator Γ i is invariant under the particleantiparticle symmetry.
Finally, the particle-antiparticle contribution to the effective mass is calculated in a similar way.