Quantum skyrmions in frustrated ferromagnets

We develop a quantum theory of magnetic skyrmions and antiskyrmions in a spin-1/2 Heisenberg magnet with frustrating next-nearest neighbor interactions. Using exact diagonalization we show numerically that a quantum skyrmion exists as a stable many-magnon bound state and investigate its quantum numbers. We then derive a phenomenological Schr\"odinger equation for the quantum skyrmion and its internal degrees of freedom. We find that quantum skyrmions have highly unusual properties. Their bandwidth is exponentially small and arises from tunneling processes between skyrmion and antiskyrmion. The bandstructure changes both qualitatively and quantitatively when a single spin is added or removed from the quantum skyrmion, reflecting a locking of angular momentum and spin quantum numbers characteristic for skyrmions. Additionally, while for weak forces the quantum skyrmion is accelerated parallel to the force, it moves in a perpendicular direction for stronger fields.

Magnetic skyrmions are textures in the magnetization which can be characterized by a topological winding number. Magnetic skyrmions were first discovered in the chiral cubic magnet MnSi [1] and subsequently in a wide range of chiral magnets, magnetic monolayers and layered magnetic systems with sizes ranging from nanometers to micrometers [2][3][4][5][6][7][8][9]. Skyrmions can be manipulated by small electric [6][7][8] and heat currents [10] which makes them interesting for future applications, such as data storage [11].
A single skyrmion in a magnetic film can be viewed as a particle. A direct consequence of its topological winding number is that its equation of motion [12] is dominated by a 'gyrocoupling' to an effective magnetic field arising from the Berry phases picked up from the spins during the motion of the texture. Furthermore, its equation of motion as a classical particle can be described by a damping constant, an effective mass and a special gyrodamping [13]. Viewing the skyrmion as a classical particle is justified in most experimental situations: the skyrmions are often large objects involving a large number of spins and the coupling to electrons in a metal or to thermal magnons will destroy effects of quantum coherence.
An interesting fundamental question concerns the quantum nature of magnetic skyrmions. Experimentally, they will mainly be of importance in insulating magnets [4] at temperatures well below the bulk gap of the underlying ferromagnet. Two important questions arise in this context: (i) How can one define and identify a skyrmion in a quantum spin system, and (ii) what are the quantum properties of such a state. In the classical case, the quantized winding number can be used to uniquely identify skyrmions. Due to Heisenberg's uncertainty principle, this is, however, not possible in the quantum case, as has been pointed out e.g. in Ref. [14], where it was also suggested to compare spin-spin correlation functions of classical and quantum spin systems to * vlohani@smail.uni-koeln.de identify skyrmion-like quantum states. While a "topological quantization" does not exist in the quantum case, one still obtains well defined "quantized" particles as stable many-magnon bound states. For the purpose of this paper we therefore define a quantum skyrmion as a stable bound state which has properties that smoothly connect to classical skyrmion states. We will show that such states are stable even in the presence of quantum tunneling and use correlation functions to show their relation to classical skyrmions.
The second question concerns the quantum properties of such a skyrmion state. The ground state properties of a single quantum skyrmion in a chiral magnet are (at least to leading order approximation) rather obvious: as their dynamics is governed by a large magnetic field, the ground state is localized in a Landau level with edge channels at the sample boundary. Corrections to this picture arise from an exponentially small lattice potential which gives rise to a bandstructure [15,16]. Only a few studies have considered quantum properties of skyrmions. Lin and Bulaevskii [17] investigated the role of a defect for skyrmions localized in a Landau band and Psaroudaki et al. [18] calculated the skyrmion mass in a quantum model, while Derras-Chouk et al. investigated the quantum collapse of a skyrmion due to tunneling processes [19]. Diaz and Arovas considered the inverse process − the nucleation of skyrmions by quantum tunneling [20]. In an interesting study, Takashima, Hiroaki and Balents [15] showed that it is possible to obtain a Bose-Einstein condensate of skyrmions. Very recently, Ochoa and Tserkovnyak [16] gave a concise overview of the quantum properties of skyrmions in chiral magnets, including their semiclassical dynamics.
In this paper we will investigate the quantum dynamics of skyrmions with internal degrees of freedom. We will consider the two-dimensional ferromagnetic quantum Heisenberg model where skyrmions arise in the presence of frustrating antiferromagnetic next-nearest neighbor interactions. Frustration stabilized skyrmions are expected to be much smaller than skyrmions stabilized by weak spin-orbit interactions, therefore quantum effects might be more important in this case. In the classical limit, skyrmion states in frustrated magnets have first been investigated by Ivanov et al. [21] and more recently by Okubo et al. [22], Leonov and Mostovoy [23], Lin and Hayami [24], and Zhang et al. [25,26]. They obtained classical skyrmion solutions by minimizing a classical energy functional. Remarkably, due to the absence of spinorbit interactions, skyrmion and antiskyrmion have exactly the same energy. Furthermore, the spins can rotate freely around the z axis. This leads to two internal degrees of freedom, the helicity φ describing rotations of spin and the skyrmion charge σ = ±1, see Fig. 1. Interestingly, the motion of the skyrmion is thereby coupled to a precession of the helicity [23][24][25][27][28][29]. Recently, the Tokura group [30] reported the discovery of skyrmions in a centrosymmetric metallic material magnet, where skyrmion formation seems to be mainly driven by frustrating interactions rather than spin-orbit coupling. A number of other candidate systems have, e.g., been discussed in Ref. [25]. In the following we will first show numerically, using exact diagonalization, that a skyrmion exists as a manybody bound state in a quantum system and that it is a stable quantum excitation. In a second step we will develop a phenomenological theory of skyrmion motion investigating both the coupling to the helicity degree of freedom and the skyrmion-antiskyrmion tunneling.

I. SKYRMIONS IN FRUSTRATED FERROMAGNETS
We consider an XXZ spin-1/2 Heisenberg model on a triangular lattice at zero temperature in a magnetic field: J 1 is a ferromagnetic nearest-neighbor coupling (set to 1 in the following) and J 2 an antiferromagnetic nextnearest neighbor coupling which can destabilize the ferromagnetic state. For J 2 > 1/3 the spin-waves of the ferromagnetic ground state have a maximum rather than a minimum of their dispersion at zero momentum and the corresponding non-linear-σ model obtains a negative spin-stiffness favouring non-trivial magnetic textures [23,24]. K > 0 effectively leads to an easy-axis anisotropy which helps to stabilize skyrmion solutions and B is an external magnetic field. The model is spinrotation invariant around the z axis.

A. Classical solutions
The classical variant of our model (1) (with a local anisotropy term K(S z i ) 2 instead of the nearest neighbor term KS z i S z j ) has been shown [23,24] to support magnetic skyrmions and antiskyrmions which have by symmetry exactly the same energy. A classical field configuration of a single skyrmion embedded in a ferromagnetic background can be described by two polar angles θ s (r − R) and φ s (r − R). Here (sin θ s cos φ s , sin θ s sin φ s , cos θ s ) describes the orientation of a classical spin and R is the position of the skyrmion. In the continuum limit, θ s (x) depends only on the distance from the skyrmion center, x = |r − R|, smoothly interpolating between a central spin pointing opposite to the ferromagnetic state, θ s (0) = π, and the ferromagnetic state, θ s (x → ∞) = 0. The in-plane angle φ s takes the simple form where φ is the polar angle in real space. For σ = 1 one obtains a skyrmion where the spin rotation follows the spatial rotation. For an antiskyrmion, σ = −1, the spin rotation and spatial rotation occur in opposite directions, see Fig. 1. Changing the parameter φ 0 induces a rotation of spins. For φ 0 = 0 or φ 0 = π/2 one obtains, for example, the so-called Néel skyrmions or Bloch skyrmions, respectively. The energy of the classical solutions is independent of φ 0 due to spin-rotation symmetry about the where Sī is the spin located at position −ri opposite to the spin Si located at position ri. Note that S x i = S y i = 0 as the quantum skyrmion is an eigenstate of Sz. The left panel displays a linear superposition of skyrmion and antiskyrmion (lz = 0) with vanishing spin currents (c.f. lowest panel in Fig. 1). The right panel displays a many-body state with finite angular momentum (lz = 2) which we identify as a quantum skyrmion. The width of the black and grey arrows indicates the size of the z-component of the spin current (largest value: 0.1 J1) on nearest-neighbor and next-nearest neighbor links, respectively. z axis. Also the energies of skyrmion and antiskyrmion, σ = ±1, are identical as one can map the skyrmion to the antiskyrmion by several symmetry transformations, for example by the product of time reversal and a rotations of all spins (but not of space) by π around the x axis.

B. Quantum skyrmion
In the following we want to investigate numerically whether the full quantum model (1), made from quantum spins with s = 1/2, also supports stable skyrmion configurations. We will use exact diagonalization results of small systems embedded in a ferromagnetic background. We have two goals: (i) to show rigorously that the pure quantum model is characterized by stable, many-particle bound states and (ii) to identify the skyrmion nature of those bound states using various correlation functions.
We consider an approximately round flake of 31 sites embedded in a ferromagnetic background, see Fig. 2. Eigenstates are characterized by the number of flipped spins, N f = S fm z − S z , where S fm z = 31/2 is the total magnetization of the ferromagnetic state. Furthermore, the flake has a 6-fold rotation symmetry with group elements {exp[iL z 2π 6 j] | j = 0, . . . , 5}, which allows us to use the angular momentum, l z = 0, 1, . . . , 5, defined modulo 6 as a second quantum number. Eigenenergies relative to the ferromagnetic state of the flake are denoted by E n (N f , l z ) (the index n refers to the enumeration scheme for eigenstates within a definite N f -l z sector). E 0 (N f , l z ) is the energy of the ground state in a given N f -l z sector. Sotnikov et al. [14] have also used exact diagonalization of a quantum magnet to search for skyrmion-like ground states in a small flake. However, in contrast to our study they use much smaller flakes, open boundary conditions and a Hamiltonian dominated by Dzyaloshinskii Moriya interactions.
From a quantum mechanical point of view, a skyrmion in a ferromagnetic background is a bound state comprising of a fixed number of flipped spins, N f . To demonstrate that such a bound state exists, we have to show that it is has a lower energy compared to a bound state with N f − N e flipped spins, where N e flipped spins have 'evaporated' and are located at the minimum E m min of the magnon band of an infinitely large ferromagnet with where δ 1 i and δ 2 i are vectors connecting the nearest and next-nearest neighbors, respectively.
On this account we demand that is the binding energy of N f spins, which can be viewed as the energy gained when N f spins come together to form a bound state instead of dispersing to infinity at the bottom of the magnon band (E B 0 (0, l z ) = 0 by definition). The binding energy is independent of the external magnetic field as we consider the stability of the skyrmion with respect to spin-conserving processes and thus compare only states with the same value of S z .
Note that the energy E 0 (N f , l z ) for fixed N f will always decrease when the size of the flake is increased. Therefore, the numerically determined value for E B 0 is a rigorous upper bound for the binding energy in the infinite system. If we numerically find negative values for E B 0 in our finite size system, this would then rigorously establish the existence of multi-particle bound states on the infinite lattice. Fig. 3 shows that for sufficiently large values of K and J 2 , binding energies become negative and have a minimum as function of N f at N f = N min f . This proves the existence of multi-spin bound states in our model. These states are our candidates for quantum skyrmions, the quantum counterparts of the classical skyrmion solutions, as discussed in more detail below.
The phase diagram in Fig. 4 gives an overview for which values of J 2 and K one can obtain skyrmion-like (b) A magnetic field is needed to stabilize the ferromagnetic and the skyrmion phase. Below the red dot-dashed line the ferromagnet (FM) is energetically unstable with respect to single spin flips. In the light-gray region the ferromagnet is the groundstate but skyrmions exist as metastable excitations. In the dark-shaded region the system can minimize its energy by the proliferation of skyrmions and a skyrmion lattice is expected to form. The black dashed line is determined from the condition that the skyrmion with the lowest energy reaches the largest number of spin flips, N f = 9, in our simulation. bound states. First, a sufficiently large frustrating interaction, J 2 0.45, is required, which leads to a negative stiffness of the ferromagnet which changes sign for J 2 = 1/3. For fixed J 2 , a small uniaxial anisotropy K is needed: it favours states where the down-spins stay together rather than fly apart. For too large an anisotropy, however, it is energetically favourable to form a single down-spin domain instead of a skyrmion of well-defined size. In Fig. 4(a) the points mark parameter values for which we obtain numerically a clear minimum when plotting the binding energy per flipped spin as function of N f (with N f ≤ 8 as our numerics is restricted to N f ≤ 9). For these parameters our numerics indicates that for a fixed small magnetization (or a fixed, sufficiently large external magnetic field) it is energetically favorable to form separate skyrmions in a ferromagnetic background rather than a single spin-down domain without internal structure. The dashed line shows that the upper phase boundary approximately follows the line E m min = 0. We find stable, finite-size, multi-spin bound states when at zero magnetic field the underlying ferromagnet is intrinsically unstable, E m min < 0, with respect to spin flips. This implies that a magnetic field is needed to stabilize both the ferromagnetic state and possible skyrmion phases thermodynamically, as is known from the classical case [22][23][24]. Fig. 4(b) investigates which fields are needed (for J 2 = 0.5). The ferromagnet is energetically stable with respect to spin flips above the red dot-dashed line. Below the lower solid black line, a single skyrmion has a negative energy compared to the ferromagnetic state and, consequently, the ferromagnet becomes unstable with respect to proliferation of skyrmions. Typically, a skyrmion lattice will form in this regime (in a tiny region of the phase diagram also a Bose-Einstein condensate of skyrmions may be realized [15]). Note that also other ordered phases (e.g. helical states) may compete with the skyrmion lattice − we did not try to investigate such phases as the focus of our investigation are the properties of single quantum skyrmions.
Therefore, we are mainly interested in the question whether single skyrmions may exist as (meta-) stable excitations above the ferromagnetic ground state. Above, we have shown that quantum skyrmions exist as stable many-magnon bound states in a sector of fixed magnetization. In a real material magnetization is not fixed due to the presence of weak spin-orbit coupling terms and dipolar interactions. Furthermore, acoustic phonons can absorb energy. As a result, any state with energy higher than the ground state will ultimately decay. Nevertheless, metastable skyrmions (light grey area in Fig. 4(b)) exist which we define by two conditions: (i) a metastable skyrmion cannot decay by spin-conserving processes, i.e., by terms included in our Hamiltonian, and (ii) it cannot decay by an (incoherent) sequence of processes where a single spin flips and the energy is lowered. The first condition is independent of the external field B and is fulfilled for skyrmions obeying the inequality (4); the second condition requires that at least one of the stable skyrmions (according to the first condition) has a lower energy than the groundstate with one flipped spin less. The second condition compares the energy of states with different magnetization and thus depends on B. Roughly, our results are generally consistent with the phase diagram obtained for a classical model by Leonov and Mostovoy [23] (for a quantitative comparison with their phase diagram, K has to be multiplied by a fac-tor of 6 and B with a factor of 2); the only main difference seems to be that our metastability regime is smaller, which may be traced back to tunneling processes that do not exist in the classical limit.
For a classical skyrmion (antiskyrmion), a clockwise rotation in space is accompanied by clockwise (anticlockwise) rotation of spins, respectively. This suggests that the spin and angular momentum quantum numbers are not independent. For the following discussion it will be useful to translate the classical wave function into a quantum mechanical one. A quantum state of a skyrmion at position R with helicity φ 0 can be approximated by the spin-coherent state with σ = ±1 for skyrmion and antiskyrmion, respectively. The angle (2) and θ i = θ s (r i − R) characterizes the tilt of the spin. This wave function is expected to become more and more accurate for larger and larger skyrmions as the magnetic texture can be locally approximated by a non-fluctuating ferromagnet. As the wave function is not expected to be accurate on a quantitative level for small skyrmions, we will use it below only for qualitative arguments. The operator which shifts φ 0 is the total spin minus the spin of the ferromagnet and an eigenstate of ∆S z is obtained from For the triangular lattice also a rotation of the position of the atoms by the angle 2π 6 n is a symmetry transformation, whereR 2π 6 n is the 2 × 2 matrix rotating the 2d skyrmion coordinate. Applying this to the state (8) we can compensate the shift of φ 0 by a rotation of the spins and find To compare to our numerical result, we consider a situation where the skyrmion is localized close to the origin R ≈ 0 (this part of the analysis will have to be modified when we consider below skyrmions without a confining flake). As for R = 0 the phase of the wave function (6) is ill-defined, we consider the limit R → 0, where the rotation of R leads to an extra phase factor e −iσ 2π 6 arising from the central spin. Taking this into account, we obtain the following selection rule for the eigenvalues l z of L z of skyrmions localized in space On the classical level, the locking of angular momentum and spin quantum number arises naturally as a skyrmion is invariant under simultaneous, identical rotations of space and spin about the z axis, and is thus an eigenstate of L z + S z , see Fig. 1. In contrast, one has to rotate space and spin in opposite directions to obtain the same configuration for an antiskyrmion, which is therefore an eigenstate of L z − S z . The unexpected shift by 1 in the quantum number is related to the fact that a simultaneous rotation does not affect the central spin, thereby rendering the relevant number of flipped spins equal to N f − 1 rather than N f . We have checked this physical picture using a Schrödinger equation for the coordinate R developed in section II B, from which, in the concluding section, we will argue that this shift arises from the confinement of the skyrmion in the finite size system considered here. Both skyrmions and antiskyrmions share the property that the xy-component of spins located on opposite sites of the skyrmion center are antiparallel. We hence cal- where Sī is the spin located at position −r i opposite to the the location r i of spin S i . Classically one finds that the maximal value of C ⊥ = 1 is obtained for spins with vanishing S z component. For the quantum skyrmions, we also obtain strong antiferromagnetic correlations, see Fig. 2. In Fig. 5 we plot the maximal anticorrelation, C ⊥ = max i C ⊥i , for groundstate spin configurations with different angular momentum quantum numbers, l z , and number of flipped spins, N f . Large anticorrelations are found whenever the selection rule given in (11) is obeyed, while the correlations are much weaker when it is violated. This confirms our analytical arguments: for a 'proper' quantum skyrmion spin and angular momentum quantum numbers are locked by our selection rule. Changing the spin also changes the angular momentum. Whether the wave function with the lowest energy (per flipped spin) obeys the selection rule or not depends on parameters as shown in Fig. 4. If the selection rule is not obeyed, the groundstate can be viewed as a 'doped' skyrmion, wherein an extra spin has been added or removed. The fact that the quantum skyrmions do not always obey simple groundstate selection rules should not be too surprising from the point of view that also for atoms or nuclei simple rules determining ground-state quantum numbers often fail.
We have also calculated a parameter related to the winding number of the skyrmion [31], where the sum is evaluated over a triangulation of the lattice (i, j, k being the sites in a triangle ∆ in the tri- angulation) such that each triangle is oriented counterclockwise. In the classical limit, W as defined above is quantized and obtains an integer value. For quantum spins this will not be the case. In the lower panel of Fig. 5 we show W as a function of l z for different values of N f . For l z = 0, 3 the groundstate is a linear superposition of skyrmion and antiskyrmion and W vanishes by symmetry. For the other values of l z we find again that |W | is maximal when the selection rule (11) is obeyed, confirming our interpretation of the bound-magnon state as a skyrmion.

II. MOBILE QUANTUM SKYRMIONS
A. Interplay of motion and helicity rotation We will now investigate the low-energy quantum dynamics of a skyrmion first ignoring the possibility of skyrmion-antiskyrmion tunneling, which will be the focus of the next subsection.
As the classical skyrmion solution is parametrized by two variables, the position R and the helicity φ 0 of the skyrmion, the low-energy Hilbert space will be spanned by these two variables and the corresponding conjugate momenta, P and S z . Our goal is to discuss a phenomenological Hamiltonian H s valid at energies well below the spin gap of the bulk phase and below the energy of possible high-energy excitation of the skyrmion. For a skyrmion with a radius much larger than the lattice constant, effects of the underlying lattice potential are exponentially small and will be ignored in this section (but are discussed briefly in subsection II B and in more detail in Appendix A). Using the fact that the skyrmion is a large object, its dynamics is expected to be governed by small values of the momentum P and small deviations of S z from its ground state value. Hence, the effective lowenergy theory can be obtained from the first few terms of a Taylor series in the momentum and the deviation of S z from its ground-state value, Here σ z is ±1 for skyrmions and antiskyrmions, respectively, M is the effective mass [13,18], F is an external force pointing in the x direction, κ and κ describe that the skyrmion mass and the effective force depend on the size of the skyrmion, and Θ parametrizes how the energy depends on the deviation of S z from the real number S 0 z (see Fig. 3). Experimentally, the force can, e.g., be realized by a small gradient in the external magnetic field. In this case, F = ∂Bz ∂Rx (S 0 z −N/2) and κ = 1 Forces can also arise, e.g., from the proximity to a sample boundary. Finally, A(R) is an effective vector potential, arising from the Berry phase of the spins which rotate when the skyrmion moves. The classical and quantum equations of motion of a skyrmion are identical to those of a particle in a huge orbital magnetic field [12,13,16], where n s is the spin density and the numerical value is given for a triangular lattice with lattice constant a. Its strength corresponds exactly to one flux quantum per area of the unit cell (for a spin 1/2), i.e. about 400, 000 T if the size of the unit cell is 1Å 2 . We will focus our analysis only on the lowest Landau level as we expect that the next Landau level has an energy larger than the spin gap of the ferromagnet, implying that the effective Hamiltonian (13) is not valid for the second Landau level [15,16]. The Hamiltonian (13) can easily be solved analytically as S z is conserved and the remaining Hamiltonian corresponds to the text-book Landau level problem. In the gauge where A = (0, −B d x, 0), the momentum k y perpendicular to the force is a good quantum number and the exact energies of the eigenstates are given by where n is the Landau level index, N f parametrizes, as above, the number of flipped spins, N 0 f = N/2 − S 0 z , and ±1 describes the solution for a skyrmion or antiskyrmion, respectively. As in the classical case, the drift velocity v s of the quantum skyrmion in real space is perpendicular to the force and simply given by the ratio of external force and magnetic field It is opposite for skyrmions and antiskyrmions, which is perhaps the easiest way to distinguish them experimentally in cases where a direct measurement of the spin configuration is not possible.
FIG. 6. Snapshots of the probability distribution, |ψ(φ0, Ry)| 2 , of the helicity φ0 and the y-coordinate of the quantum mechanical wave function, for times t = 0, 20, 40, 150, 300, of an antiskyrmion driven by a magnetic field gradient (initial condition: ). For short times the angle φ0 grows linearly in time, but its wave function also spreads. For long times, the wave function splits into distinct wave packets. This does not reflect an interference effect but arises because the velocity depends on the number of flipped spins, N f (see Eq. (16)). It therefore indicates a perfect entanglement of position and N f . For a skyrmion the same result is obtained with Ry → −Ry.
A measurement of the helicity φ 0 of the skyrmion (e.g., by an electron microscope which is sensitive to the inplane orientation of spins [2]) will result in a collapse of the wave function to a state with fixed φ 0 (within measurement accuracy) described by a superposition of states with different values of N f . As −N f is the conjugate momentum to φ 0 , one will subsequently observe a precession of φ 0 with the group velocity ∂ t φ 0 ≈ − ∂ ∂N f E ± N f ,ky,n . Remarkably, the motion of the skyrmion induces an additional precession of the helicity and therefore of the inplane spins on top of the quantum mechanical precession for F = 0, where we used that ∓ ky B d = R x for a wave function in the lowest Landau level. More precisely, the result shown above holds only when the force is turned on adiabatically. Switching the force suddenly excites higher quantum numbers n (thus possibly leaving the range of applicability of our low-energy Hamiltonian).
In several studies, e.g. [23][24][25][27][28][29], it has pre-viously been observed that in classical models the helicity dynamics is coupled to the translational motion of skyrmions and antiskyrmions when skyrmion motion is induced by various forces, most notably spin-orbit torques [24,25,28]. These dissipative forces are, however, associated with extra channels of decoherence not captured in our effective model. A direct comparison can hence only be made to the dynamics induced by field gradients studied in Ref. [29], where, however, no detailed analysis has been given.
We have checked that straightforward classical simulations (not shown) of our model reproduce a precession of the helicity proportional to F 2 if a small field gradient is applied. For J 1 = 1, J 2 = 0.5, K = 0.05, and an average magnetic field chosen to describe about 7 flipped spins we find both for skyrmions and anti- , which is consistent with (2κ + κ)M ≈ 12 J 1 /a 2 (or M ∼ 30 J 1 /a 2 using that κ = 1/N f and assuming κ ∼ κ ). The simulation result has been obtained using the standard Landau-Lifshitz-Gilbert equation in the limit of weak damping α, and we find thatΦ 0 is approximately independent of α in the long-time limit. Note, however, that the long-time and vanishing-damping limits do not necessarily commute. A value of κ ∼ 1/N f is also consistent with the assumption that the mass is proportional to N β f which implies κ = β/N f (Ref. [13] obtains β = 2 from a classical calculation). The parameter Θ can be obtained by a straightforward fit to the S z dependence of the skyrmion energies available to us from the exact diagonalization results. For the parameters quoted above, we find Θ ≈ 20/J 1 . The value of Θ has also been estimated previously for classical models in Ref. [25,32]. In our units their formulae translate to Θ = 1/(3K) and Θ = 1/(6K), respectively, which differs from our result. Nevertheless, this concludes our estimation of all parameters of the effective model defined in Eq. (13).
Beyond the classical effects, the quantum mechanical model predicts new phenomena. First, due to the Heisenberg uncertainty relation a measurement of the helicity with precision δφ 0 leads to an uncertainty in the conjugate variable N f − N 0 f and a subsequent quantummechanical spread of the wave function. This effect rapidly washes out the precession of the helicity, see Fig. 6, where less than half of a 2π precession is observed before the wave function covers all angles. To observe at least a full rotation of φ 0 , the condition |κ + 2κ |δφ 0 1 must be fulfilled. Using the parameters of Fig. 6, characteristic of a skyrmion of O(10) flipped spins, the required force would be very large, of the order of 0.1 J 1 /a (corresponding to a field gradient of 0.01 J 1 /a). Such a large force cannot be realized by a field gradient in a bulk sample (it can, however, arise from confining forces at the edge of the sample). We conclude that for small skyrmions consisting only of O(10) flipped spins, the quantum mechanical spread of the wave function is likely to dominate the drift of the angle obtained classically. The second effect is that the position of the particle R y gets entangled with the magnetization of the skyrmion. Similar to a Stern-Gerlach setup, where the trajectories of particles depend on S z , the velocity of the skyrmion, Eq. (16), depends on the discrete variable N f , and consequently the wave function eventually splits into separate wave packets distinguished by their local magnetization, see Fig. 6. This is also a purely quantum mechanical effect arising from the quantization of the magnetization. As the magnetization for each of these wavepackages is fixed, the conjugate variable, the helicity, shows maximal uncertainty. To observe the effect, we require that δt∆v a, where ∆v = κ F B d is the velocity difference of two skyrmion states differing by ∆N f = ±1 and δt is the time scale on which the skyrmion propagates in a quantum coherent way. From this condition and the estimate κ ∼ 1/N f , we obtain the requirement δt N f F a .

B. Skyrmion-Antiskyrmion tunneling
We will now consider the consequences of skyrmionantiskyrmion tunneling. Using that 1/Θ in Eq. (13) is expected to be much larger than the exponentially small tunneling rate, we assume in the following that the magnetization of the skyrmion is fixed to its ground state value. The effective Hamiltonian for the tunneling problem is therefore given by (18) where ∆ R is an operator encoding the position dependence of the tunneling amplitude (specified in more details below in Eq. (29)), σ ± = 1 2 (σ x ± iσ y ) are operators which induce transitions from skyrmion to antiskyrmion and back, and V (R) is the periodic potential generated by the underlying lattice of the spins. Even in the absence of tunneling, such a periodic potential delocalizes quantum particles in the lowest Landau level and leads to a finite dispersion. V (R) is, however, exponentially small in R s /a, the ratio of skyrmion radius and lattice constant [15]. For our skyrmions we show in Appendix A 2 that V (R) is indeed tiny (∼ 10 −4 J 1 for one set of parameters) and much smaller than estimates of the tunneling rates. We therefore set V (R) to zero in this section and discuss effects of a finite V (R) only in Appendix A.
It is tempting to use in Eq. (18) a vector potential, e.g., of the form A = B d 2 R ×ẑ with ∇ × A = B dẑ in combination with a constant ∆. While this choice of the vector potential is completely appropriate in the absence of tunneling, it leads to unphysical results (a single localized state at the origin of the coordinate system) in the presence of tunneling. To understand the problem it is useful to realize that the tunneling event can be viewed as a sudden sign change of the vector potential. Such a sign change creates an unphysical electric field spike E =Ḃ d 2 R × z growing linear in distance from the origin. It is thus imperative that we rederive more carefully the vector potential of the skyrmion. It originates from the Berry phases of the underlying spin-1/2 system. Parametrizing each spin with a unit vectorn i with angles θ i and φ i , the Berry phase action of the spins can conveniently be computed using a singular vector potential a s (n). We use a gauge choice where a s = 1−cos θ sin θφ withφ being the unit vector in φ direction (this gauge choice is also compatible with the wave function (6)). Using that n i =n(r i − R), we obtain for the Berry phase action of the spin-1/2 system [33] with the vector potential A(R) is singular when the center of the skyrmion R where the spin points down is located exactly at the location of a spin, R = r i , because we used a gauge choice for which a s is singular at θ = π. Our choice of the gauge has the advantage that the vector potential of skyrmion and antisykrmion (σ z = ±1) are opposite to each other. When calculating A and the magnetic field corresponding to A one has to take into account the singularity and one finds We obtain a constant uniform magnetic field B d and for R = r i an extra contribution localized in δ functions and carrying exactly one flux quantum per lattice site. The magnetic field integrated over a unit cell, U C B d 2 r vanishes, as expected for a vector potential which is periodic in space. Note that A is periodic, A(R + r i ) = A(R), as the linear term in R in the first term in Eq. (21) is exactly canceled by a similar contribution from the second term (assuming that the lattice approximately has the shape of a disc). We have checked this property numerically.
In the absence of tunneling, one can gauge away the singular part of the magnetic field using a singular gauge transformation. In the presence of tunneling, however, such a gauge transformation will modify the tunneling term in a singular way.
To obtain eigenstates and low-energy spectrum of the Hamiltonian (18), we first construct for ∆ R = 0 eigenstates of momentum k both for skyrmions and antiskyrmions. Starting from the well-known Landau levels in a constant magnetic field, ψ ∼ e − B d 4 , we first perform a singular gauge transformation to obtain the corresponding eigenstate for the vector potential (21), where the only difference between skyrmions and antiskyrmions is the sign of the phase factors. As the vector potential A is a periodic function of R, one can simply translate the wave function by a lattice vector r i to obtain another (in general not orthonormal) eigenstate localized around r i . An approximation for the corresponding many-body wave function localized around the site r i is Note that in the many-body wave function the singular terms e −σiφ(R−ri−rj ) at lattice points r i + r j cancel exactly with a singular contribution in the definition of |σ, R, N f in Eq. (6). Importantly, the angular momentum of the wave function is given by Eq. (11), where the origin of the −1 term can be traced back to the singular gauge transformation. The momentum eigenstates (without normalization) are given by This eigenstate is unique within the lowest Landau level which contains exactly one state per flux quantum. Thus, for skyrmions and antiskyrmions, there is precisely one state per unit cell of the lattice, each. As the tunneling matrix element is much smaller than any Landau level spacing, we can ignore Landau level mixing and compute the tunneling matrix elements projected onto the lowest Landau level directly, The denominator is needed as we used non-orthonormal wave functions. Within the lowest Landau level of skyrmions and antiskyrmions, the Hamiltonian in momentum space is then simply described by with eigenvalues The tunneling process is expected to be local and also the overlaps α ri decay rapidly with distance (α ri = e −B d r i 2 4 e iπ r i 2 a 2 , and numerically equal to 1, −0.16, −0.004, 0.0007, −3 · 10 −6 onsite and for nearest, next-nearest, third-and forth next neighbors). Note that ri α ri = 0 reflecting the fact that the wave function ψ ±,k (R) carries angular momentum ±1 and therefore has to vanish for k = 0.
Tunneling is constrained by the crystalline symmetry and the relative angular momentum of skyrmion and antiskyrmion encoded in the spin-wave function |σ, R, N f . The difference of the angular momentum of the spin-part of the wave function is 2N f mod 6. As a phenomenological ansatz we expand the tunneling matrix element in lattice harmonics, keeping the lowest-order term allowed by symmetry,  where G n are six reciprocal lattice vectors obtained by rotating the first one, G 1 , by the angle 2π 6 (n − 1). The tunneling rate δ is expected to be exponentially small in the skyrmion size. For N f = 1 mod 3 we can obtain an estimate of the tunneling rate using the splitting of the lowest two energy level within our exact diagonalization result. For N f = 7, J 2 = 0.5, K = 0.05 we find, for example, in the l z = 0 sector a sizeable tunneling splitting ∆E t ≈ 0.05 J 1 . Since this splitting arises from tunneling of localized skyrmion and antiskyrmion, it is approximately equal to 2δ 0 0 ≈ 3.1δ, from which we estimate δ = ∆E t /2δ 0 0 ≈ 0.015 J 1 . Note that this tunneling rate is much larger than our estimate for the amplitude of the periodic potential V 0 ≈ 7 · 10 −5 J 1 derived in Appendix A 2 for the same parameters.
In Fig. 7 the resulting bandstructure is shown N f = 0 mod 3 and N f = 1 mod 3, while Fig. 8 displays the bandstructure for N f = 2 mod 3. The qualitative difference between the three bandstructures can be traced back to the angular momentum of skyrmions and antiskyrmions.
The bandstructure for N f = 0 mod 3 in the left panel of Fig. 7 is regular and non-singular. This is, perhaps, surprising as our analysis of localized skyrmions and antiskyrmions revealed that in this case skyrmion and antiskyrmion are in different angular momentum channels l z = ∓1 and the localized skyrmions do not tunnel into each other. In contrast tunneling for k = 0 is possible and non-singular. Apparently, the mobile skyrmions transfer the angular momentum to the emergent magnetic field when tunneling. Mathematically, we find that tunneling is dominated by the nearest-neighbor tunneling matrix element δ 0 ri while local tunneling vanishes, δ 0 0 = 0. In  Fig. 7, the bandwidth is not linear in the tunneling rate δ but proportional to δ 2 /ωc. The splitting of the two bands is proportional to δ 3 /ω 2 c . The right panel shows the trajectory of a classical particle in a magnetic field with random changes of the sign of the charge to mimic tunneling skyrmion-antiskyrmion tunneling events.
the limit k → 0, the k dependent tunneling is highly singular, This implies that the eigenfunction of the effective Hamiltonian (27) obtains a Berry phase 2π when circling in momentum space around k = 0. This Berry phase cancels, however, exactly, a corresponding singular Berry phase arising from the definition of the momentum eigenstates, Eq. (25) and the resulting wave function and bandstructure are both smooth and non-singular for k → 0. For N f = 1 mod 3 the bandstructure is completely changed. As shown in the right panel of Fig. 7, we obtain a parabolic band-touching at the Γ point. This is a direct consequence of the finite angular momentum of the wave function. The quadratic band touching is thereby associated with a Berry phase of 2π for an adiabatic path around the Γ point. We find numerically that the bandminima are now located at the two K points. Thus, at low energies, the skyrmion obtains a new quantum number describing in which band-minimum the quantum skyrmion is located.
For N f = 2 mod 3 we obtain numerically that the tunneling from the lowest Landau level of the skyrmion and the antiskyrmion is not possible, δ k = 0. We have not been able to find an analytic argument proving this. To calculate the band-structure in this case, we therefore have to consider Landau-level mixing induced by the tunneling process. We have used a discretized version of the Hamiltonian (18) to calculate numerically the resulting bandstructure. As is shown in Fig. 8, the bandwidth is proportional to the square of the (exponentially small) tunneling rate δ in this case and thus much smaller than for N f = 0, 1 mod 3 where a dispersion linear in δ was obtained. The energy splitting of the two bands is even smaller, of order δ 3 . Due to the finite angular momentum we obtain again a quadratic band touching at the Γ point where also the minimum of the dispersion is located. As the effects of tunneling are strongly suppressed in this case, one has to reconsider the effects of a tiny periodic potential V (R), which can also induce a finite bandwidth. We discuss this case in Appendix A and find quantitative but no qualitative changes of the bandstructure for N f = 2 mod 3 .
All results in this section have been derived under the assumption that the Hamiltonian (18) is valid. One assumption which me made is that the tunneling ∆ is local and just a function of the coordinate R but not of the momentum. While this is an ad-hoc assumption, we expect that all qualitative features of the bandstructure, which rely purely on symmetry, will remain the same if more complicated tunneling terms are considered.

III. CONCLUSIONS AND OUTLOOK
We have shown that a skyrmion in a frustrated magnet is a quantum particle with a list of rather unusual properties. Most importantly, the motion of the skyrmion and the internal degrees of freedom are directly coupled. As was already known from the classical theory, the helicity couples to the motion of the skyrmion which leads to a characteristic precession of the spin when the skyrmion is moving in the presence of an external force. For small skyrmions this effect is difficult to observe due to a combination of Heisenberg's uncertainty principle and the quantum mechanical spread of the (helicity-) wave function. Furthermore, position and spin become strongly entangled during time evolution which leads to a characteristic quantization of the skyrmion velocity in the presence of a force.
The helicity and simultaneously the position of a moving skyrmion can be measured using, e.g., an electron microscope [2]. It would, for example, be interesting to study how the quantum mechanical spread of the helicity and the entanglement of spin and position is affected by the presence of thermal magnons or by continuous weak measurements due to the electron microscope itself.
In the absence of an external force and of tunneling, a skyrmion is localized in the lowest Landau level. In this case skyrmion-antiskyrmion tunneling can delocalize the particle. A semiclassical explanation of this effect is shown in the right panel of Fig. 8: During a tunneling event the effective charge of the skyrmion changes, allowing for skyrmion motion not confined by cyclotron orbits. The corresponding bandwidth is thereby naturally set by the tunneling rate. The internal angular momentum of skyrmion and antiskyrmion imposes, however, strong constraints on possible tunneling events. As for skyrmions, spin and angular momentum are locked to each other. This implies that the bandstructure changes drastically when a single spin is added or removed from a skyrmion, see Figs. 7 and 8. When the number of flipped spins N f is 0 or 1 mod 3, the bandwidth is proportional to the tunneling rate, while it is quadratic in the tunneling rate for N f = 2 mod 3.
It is interesting to consider the response of the quantum skyrmion to an external force or a confining potential. For simplicity, we consider the simple case, where N f is a multiple of 3, where the bandstructure has a unique minimum at the Γ point. A weak force, F δ/a, therefore simply leads to an acceleration of the particle parallel to the external force and to Bloch oscillations. As the particle is a superposition of skyrmion and antiskyrmion states, the quantum skyrmion effectively carries a topological charge 0 and does not see the emergent magnetic field. For large forces, F δ/a, the picture changes completely. In this case, tunneling is suppressed. Skyrmion and antiskyrmion move with velocity ± F B d , see Eq. (16), perpendicular to the external force in opposite directions. The crossover between the two regimes is driven by the force-induced Zener tunneling between the two states. Using the splitting of skyrmion and antiskyrmion trajectories, one can use an external force to implement Stern-Gerlach type of experiments using field gradients. Close to a sample boundary, which will likely acts as a repulsive force for skyrmion and antiskyrmions, one can expect chiral edge channels for skyrmions and antiskyrmions running in opposite directions and its an interesting open problem how these edge states merge with the bulk bands and how this affects the scattering of quantum skyrmions from sample boundaries. Similarly, in a weak confining potential, the ground state is unique with vanishing angular momentum, while a strong confining potential, such that V 0 a 2 δ, leads to the doubly degenerate groundstate with angular momentum ±1, as observed for the exact diagonalization of small finite systems.
Our study has focused on single skyrmions -the next step is to consider pairs of skyrmions and their mutual interactions. The interaction potential has an oscillating sign [24] as the magnons in the ferromagnetic state have a minimum at finite momentum. This will lead to the formation of boundstate and, for a finite skyrmion density, to crystalline phases.
Appendix A: Effect of lattice potential on the bandstructure The discrete atomic lattice breaks the continuous translational symmetry and thereby induces a periodic potential V (R) for the skyrmion position R. Such a periodic potential is potentially important as in the presence of both magnetic field and periodic potential, Landau bands are not completely flat but obtain a dispersion. In order to estimate the impact of this effect, we derive an approximation for V (R) and include it in our bandstructure calculation.
The atomic triangular lattice or our model (see Fig. 2) is characterized by the lattice constant a. A skyrmion is a smooth texture with radius R s a. It effectively averages over many lattice sites and consequently Fourier components of the resulting periodic potential are exponentially small in |q i |R s [15,16], where q i are the reciprocal lattice vectors. Ergo, the periodic lattice potential of the skyrmion is with exponential precision described by only the six shortest reciprocal lattice vectors. Its shape is completely fixed by symmetry, The prefactor V 0 is the only free parameter, expected to be exponentially small in R s /a. The periodic lattice potential V (R) has to be added to to our phenomenlogical Hamiltonian, where σ 0 is the 2×2 identity matrix. In the following, we will first discuss how a finite V 0 affects the bandstructure of the quantum skyrmions both for small and large values of V 0 , then we will show that V 0 is very tiny for the quantum skyrmions discussed in our paper.

Effect on the bandstructure for arbitrary potential strength
The bandstructure of quantum skyrmions described by Eq. (A2))can be obtained from a straightforward exact diagonalization of a discretized version of the Hamiltonian. Our discussion will focus on the case V 0 > 0, the sign obtained in subsection A 2.
In Fig. 9 we show the resulting bandstructure for N f = 0, 1 mod 3 and three values of V 0 /δ, where δ is the skyrmion-antiskyrmion tunneling rate δ defined in Eq. 29. For small δ and V 0 , the dispersion in units of the tunneling rate, E k /δ, depends only on the ratio V 0 /δ. As expected, for V 0 δ (the relevant limit within V0 is the strength of the periodic potential and δ the skyrmionantiskyrmion tunneling rate. In our model system, we estimate that V0 δ. FIG. 10. Bandstructure for different values of V0ωc/δ 2 for N f = 2 mod 3. In this case E k /δ 2 is just a function of V0ωc/δ 2 . Qualitatively, the shape of the bandstructure does not change significantly upon increasing V0ωc/δ 2 ; the quantitative changes can be significant, though. For high V0ωc/δ 2 , the bandstructure becomes independent of N f . our model, see below) the effect of the periodic potential can be neglected, small corrections linear in V 0 /δ do not change any qualitative features of the bandstructure. For V 0 δ the dispersion is mainly determined by the periodic potential. The splitting of the two bands, however, is governed by the skyrmion-antiskyrmion tunneling. For N f = 0 mod 3 the minimum of the bands is always at the Γ point, while for N f = 1 mod 3 it moves from the K points towards the Γ point.
For N f = 2 mod 3, see Fig. 10, the dispersion in the absence of a periodic potential is not proportional to δ but to δ 2 /ω c , instead, see also Fig. 8. The lattice potential is expected to create a first order perturbative correction proportional to V 0 but not lift the skyrmionantiskyrmion degeneracy. For a combination of weak tunneling and lattice potential, therefore, the bandstructure E k ω c /δ 2 is expected to be a function of V 0 ω c /δ 2 . The numerical results for the bandstructure are shown in Fig. 10 for various values of V 0 ω c /δ 2 . The effect of the lattice potential V (R) is less pronounced in the qualitative sense, although, there can be significant quantitative changes when V 0 is varied for a fixed δ.

Classical approximation for the lattice potential
In order to estimate the prefactor V 0 of the lattice potential we can employ a classical approximation of the magnetization by replacing the quantum mechanical spin operators S in the Hamiltonian, Eq. (1), with classical Heisenberg spins m with m = 1/2. We can exploit that the skyrmion position is fixed by symmetry if it is initialized on a highly symmetric point. Thus we use standard relaxation algorithms without artificially fixing any spins to calculate the energy of a skyrmion which is centered (i) on a lattice site, (ii) on a plaquette, and (iii) on a bond between lattice sites, respectively. Although we only require the energies of two positions for fitting V 0 , we can use the third position to validate the simplified ansatz for the potential, Eq. (A1).
The skyrmion that we consider in Sec. II B is stabilized for J 1 = 1, J 2 = 0.5, K = 0.05, a = 1, and 7 flipped spins. In the classical system, we have to tune the external magnetic field such that the latter condition is fulfilled. Hence, we choose it such that the skyrmion centered on a lattice site has a difference in the total magnetization of ∆m z = 7 with respect to a fully polarized state. Keeping the magnetic field fixed to this value, we obtain from the energy difference between the skyrmion centered on a bond and a plaquette V 0 = 8.15 · 10 −5 and for centered on a site and a plaquette V 0 = 7.90 · 10 −5 . The good agreement between the two estimates with a precision of 3% shows that Eq. (A1) is well justified. By performing atomistic simulations with a small applied electric current, we could furthermore drive the skyrmion to other intermediate positions and evaluate its energy which is in excellent agreement with the simple ansatz in Eq. (A1). ∆m z varies slightly as function of position within the classical approximation while it is integer valued and conserved in the quantum theory. We therefore performed a second calculation, where we adjusted the external magnetic field for each position of the skyrmion such that ∆m z = 7 is fixed. This yields V 0 = 7.20 · 10 −5 (V 0 = 7.41 · 10 −5 ) when comparing bond-and plaquette centered skyrmions (from the difference of site-centered and plaquette-centered skyrmions). Both numbers agree up to an error of 3% and do not differ significantly from the results obtained for constant magnetic field.
As we estimate from the exact diagonalization result that the tunneling rate for the same set of parameters is 0.015, our results indicate that the potential is much smaller than the tunneling rate, V 0 δ.