{'fraction_non_alphanumeric': 0.041652224147891395, 'fraction_numerical': 0.013922588099364528, 'mean_word_length': 4.188998800959233, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Magnon band structure of skyrmion crystals and stereographic projection approach

Using semiclassical method combined with stereographic projection approach, we investigate the magnetic dynamics of the skyrmion crystal (SkX), formed in planar ferromagnet with both Dzyaloshinskii-Moriya interaction and uniform magnetic field. The topologically non-trivial ground state of SkX is described in stereographic projection by the complex valued function with simple poles at skyrmions' positions. We use the earlier proposed ansatz for this ground state function in the form of the sum of individual skyrmions. The dynamics follows from the second variation of the classical action. Numerical analysis yields the magnon band structure of tight-binding form in accordance with previously known results. There are two sets of bands, one set with a flat dispersion, topologically trivial and rapidly evolving with magnetic field. Another set is robust to magnetic field, characterized by pronounced dispersion and with the Berry curvature which may be sign-reversal in the Brillouin zone. The developed theory can be straightforwardly generalized for the analysis of magnetic dynamics in topological spin structures of other types.


I. INTRODUCTION
Magnetic skyrmions are one of the most popular examples of topologically protected configuration of local magnetization [1][2][3] . The topological protection and relatively small size of skyrmions make them promising candidates for development of new types of magnetic random access memory 4 and programmable logic 5 devices.
The theoretical concept of magnetic skyrmions appeared in the work by Belavin and Polyakov (BP) 6 , where skyrmions were found as exact solutions of the continuous model of two-dimensional ferromagnet at T = 0. It was also shown there that these solutions with non-zero topological charge are metastable. It was found later that Dzyaloshinskii-Moriya (DM) interaction and external magnetic field can stabilize skyrmions 7,8 , and they tend to form the skyrmion crystal (SkX).
The experimental observations of SkX were reported more than a decade ago by means of neutron scattering 9,10 and Lorentz transmission electron microscopy 11 . Nowadays there is a large number of compounds and systems, where different types of skyrmionic lattices with different crystal symmetries were confirmed, see the review 12 .
The dynamics of skyrmionic systems is both interesting and difficult question, still under active investigation. A single skyrmion can be viewed as a circular domain wall of small size, i.e. a magnetic bubble. One can then use the Thiele equation 13 for the description of single skyrmion motion 14 . At the same time, the skyrmion dynamics is not limited by the motion of its center. Full equations of motion for a single skyrmion should also take into account internal excitations [15][16][17] . These excitations correspond to distortions of the form of skyrmions: dilatation, elliptical deformation etc.
Just as the physics of solids is not reduced to the physics of the atoms in the lattice, so too the SkX dynamics is not reduced to the motion of individual skyrmions. It was shown that the low energy dynamics of SkX in-cludes the Goldstone mode 18 associated with displacement of SkX as a whole. Three types of SkX excitations, breathing, clockwise and counterclockwise modes, have been derived numerically on a lattice 19 . Another interesting way to describe the low energy SkX dynamics in terms of phason excitations was suggested in 20 .
The spatial periodicity of SkX configuration translates into the notion of the Brillouin zone (BZ) in reciprocal space, which in turn poses the question of dispersion of SkX excitations. This question was addressed numerically on a discrete lattice for atomic scale Néel-type skyrmions in 21 . Another calculation in the framework of linear sigma model was done in 22 for some realistic compounds, taking into account also dipole-dipole interaction. The dependence of the band structure on external magnetic field was studied in 23 .
One of the most interesting features of the magnon band structure in SkX is its nontrivial topological properties. It was demonstrated 21-23 that certain low-energy bands have non-zero Chern numbers. It was further shown that this property results in the appearance of magnon edge states 21,23 , that could be useful for magnonics. Besides, nontrivial Berry curvature and Chern numbers should lead to thermal Hall effect 24,25 .
In this paper we combine the semiclassical method and the stereographic projection approach for the investigation of SkX dynamics. Such combination is primarily useful at low temperatures, when the local magnetization saturates to its maximum value. The stereographic projection describes the direction of the equilibrium magnetization in terms of complex-valued function, f . The magnetic energy acquires then the form of highly nonlinear functional, whose extremum is well reached by ansatz borrowed from the original BP paper. 6 Having obtained a suitable description of the classical configuration of SkX, we find the dynamics from the second variation of the classical action. 26 The proposed semiclassical analysis leads to equations of motion structurally identical to those found previously in linear spin-wave theory formalism. 15 Our numerical analysis of these equations yield the dispersion laws of low-lying excitations in good agreement with previous findings. [21][22][23] The subsequent analysis of topological properties of the low-lying bands by link-variable method 27 provides a good accuracy in evaluation of the Berry curvature and the Chern numbers.
Among the advantages of the proposed method we can name a few. One is the possibility to work consistently in the continuum model, without digressing to lattice formulation in order to find the equilibrium magnetization configuration. Another is the principal possibility to investigate the dynamics of irregular skyrmion structures, e.g. liquids, which is described by the same equations of motion with a different stereographic projection function. The third one is better control over the so-called zero mode, describing the motion of skyrmions system as a whole. We plan to explore the two latter features in subsequent publications.
The rest of the paper is organized as follows. In Section II we describe our model, the general formalism of semiclassical approach and stereographic projection, and derive the equations of motion. The discussion of several important technical issues is found here as well. In Section III we briefly describe the construction of static SkX configuration in stereographic projection by the ansatz proposed earlier in 28,29 . Here we also discuss the way to calculate corrections to our trial ground state, and provide corresponding formulas. The magnon dispersion of SkX is calculated in Section IV for a number of values of external magnetic field in our model, the analysis of band dispersion in terms of tight-binding approximation is made here. In Section V we calculate the Berry curvature and Chern numbers for low-lying magnon bands and qualitatively discuss the results. Our final remarks and conclusions are presented in Section VI. Further technical details are discussed in Appendices A and B.

A. Energy and Lagrangian
We consider a minimal model of the two-dimensional (2D) non-centrosymmetric ferromagnet possessing skyrmions, with the energy of the system E = d 2 r E and the energy density here C and D are exchange and DM constants, B is external magnetic field perpendicular to the plane; µ = 1, 2 and i = 1, 2, 3; totally antisymmetric tensor µij and summation over the repeated indices is assumed. This simple form of DM interaction implies cubic crystals, and in general one considers DM interaction via Lifshitz invariants depending on the particular symmetry of the system. For zero temperature considered here we have a constraint S i S i = 1 for local magnetization S(r). We assume D, B > 0 and added unity in the last term of (1) in order to have zero energy for the anticipated uniform ground state, S 3 = 1, at D = 0. The model (1) has two characteristic length scales, C/D and C/B. Taking values of C and D as basic parameters of the material, we measure length in units of l = C/D. The energy density, measured in units of Cl −2 = D 2 /C, can be reduced to the dimensionless form : The dynamics of the magnetization without damping effects is usually described by the Landau-Lifshitz (LL) equation:Ṡ where γ 0 is a gyromagnetic ratio and H = δE/δS is an effective magnetic field. The LL equation (3) is a convenient tool to numerical calculation and micromagnetic simulations, but we prefer below the analysis based on the Lagrangian formalism. The Lagrangian density L = d 2 r (T − E) includes the kinetic term: where ϕ and θ define the magnetization direction S = (cos ϕ sin θ, sin ϕ sin θ, cos θ) and a 0 is the distance between localized magnetic moments, arising here after passing to continuum model.

B. Stereographic projection approach
We parametrize the magnetization in the stereographic projection approach as where the complex-valued function f = tan(θ/2) e iϕ and its complex conjugate,f , depend in our 2D case on z = x + iy andz = x − iy, where x and y are spatial coordinates. We note in passing that in the alternative representation 15 of (4), T = − a 2 0 A ·Ṡ, the gauge field A coincides up to a prefactor with f , namely The representation (5) is constraint-free, |f | ∈ [0, ∞). In the following we also use another field variable ψ = sin(θ/2)e iϕ , whose Jacobian of transformation from the initial spherical angles is especially simple, dS = sin θ dθ dϕ = 2i df df /(1 + ff ) 2 = 2i dψ dψ. Notice however that, due to the inconvenient constraint |ψ| ≤ 1, we use below only the differential form, dψ = df /(1 + f 0f0 ), for description of small fluctuations around the static unconstrained configuration, f 0 .
The complex function f , with the additional condition of homogeneity, f → cst at r → ∞, represents a map of one two-dimensional sphere onto another S 2 → S 2 . Such maps can belong to different homotopy classes π 2 (S 2 ) = Z, with an integer number classifying the degree of the map called topological charge. In terms of the complex function f (z,z) the topological charge is given by : where ∂ z = (∂ x − i∂ y )/2 and ∂z = (∂ x + i∂ y )/2. The kinetic part of the Lagrangian is given by whereas the energy (2) is Here and below we use the curly brackets in order to designate the terms originated from the DM term in Eq. (1).
The state corresponding to the (local) energy minimum is found by variation of (8), leading to the equation: which is time-independent. Generally, nonlinear equations may have multiple solutions whose search is difficult and analytically unavailable in most cases. Equation (9) is not an exception in this sense. We postpone the discussion of our trial function to (9) until Section III A, and now discuss a general formalism of the dynamical fluctuations around the static configuration.

C. Dynamics
In this subsection we will describe the well known formalism of semiclassical quantisation 26 and its application to our situation. We assume that we found the static configuration f 0 (z,z), obeying the Eq. (9).
A simplistic way to consider the dynamics is to allow small time-dependent fluctuations around the classical solution f (t, z,z) = f 0 (z,z) + δf (t, z,z). The Lagrangian then takes the form of functional Taylor series: where L n is the nth variational derivative in f and dots represent higher order terms beyond the scope of our paper. The first term vanishes, L 1 [f 0 ] = 0, for f 0 satisfying Eq. (9).
The kinetic part (7), in view of T [f 0 ] = 0, reads This expression has an unconvenient denominator, which can be eliminated by local redefinition, ψ = δf /(1+f 0f0 ). This choice of new variable is corroborated by the trivial Jacobian for ψ,ψ, mentioned above. The kinetic part takes now a simple form: It is known however that the above simplistic way of treating fluctuations is ill-suited for description of the socalled zero modes. In our case this mode corresponds to the translational symmetry of the skyrmion configuration on the infinite plane, f 0 (r) → f 0 (r − R). The appropriate treatment of fluctuations including the zero mode is discussed in 26 and in our case is given by here f 0 ≡ f 0 (r−R(t)) acquires the dynamics of its own. 15 The velocity ∂ t R(t) is also assumed to be small, although different by construction from ψ. We provide further details in Appendix A, and set ∂ t R(t) ≡ 0 from now on. Proceeding along these conventions, we represent the second order Lagrangian as with the Hamiltonian operatorĤ of the form Here U , V and A = e x A x + e y A y are functions of f = f 0 (r) and its gradients: The Lagrangian (14) results in the Euler-Lagrange equation: with σ 3 the third Pauli matrix. The energy of normal modes, n , is found from The similar equation for excitations above the skyrmionic ground state has previously appeared in the literature, see e.g. 15 . Our approach generalizes this equation to the multi-skyrmion ground state by describing it in terms of the function f . A particular case of 15 is reproduced by using in (15) the single skyrmion solution for f 0 , as discussed in the next Section.

D. Qualitative discussion of the Hamiltonian
Before proceeding further with our calculation, let us qualitatively discuss the Hamiltonian. The form of the Hamiltonian (15) resembles one of Bogoliubov-de Gennes in the theory of superconductivity, although it describes bosonic excitations instead of fermions. After standard bosonic representation of spin operators, the Hamiltonian similar to (15) appears in the linear spin-wave theory of antiferromagnets and, generally, of ordered magnets with non-parallel alignment of spins.
What differs SkX from other non-collinear magnets is the existence of the gauge field A. The definition (18) of A includes two terms, arising from exchange interaction and DM interaction. One can easily see that the part, stemming from DM interaction, is regular over the whole plane. The first term in (18) is more complicated and can be rewritten as: We show below that near the center of one skyrmion placed at the origin the profile function f contains singular and regular parts, f SkX i/z + f reg . It leads to an appearance of singular part in (21): This in turn leads to delta-function singularity for gauge field intensity, B = ∇ × A: For a single skyrmion we have f → 0 as r → ∞. The application of the Stokes' theorem gives then the zero total flux 15 For multi-skyrmion configuration the field intensity B has a set of δ-functions placed at centers of skyrmions, and this set is infinite and regular in SkX case. One can show that the zero flux property should hold also for multi-skyrmion configuration, when the integration in (24) is performed over polygons, referring to individual skyrmions (presumably after Voronoi partition). Numerically we checked that the zero flux is obtained in our approach after the integration over the unit cell in SkX case.
Let us ignore for a moment the existence of anomalous terms, i.e. set V = 0 in (15) The simplified Schrödinger equation reads which describes the motion of the particle in the nonuniform magnetic field and electrostatic potential, both with the same spatial periodicity. The δ-function contributions in (23) can be removed by the unitary transformation with z j corresponding to centers of skyrmions, see (32) below. The regular part of the pseudo-magnetic field in (23) is smooth and could be well approximated by the uniform field, with two flux quanta per unit cell. At first glance such commensurate field can be ignored, following the conclusions by Azbel and Hofstadter. However, a closer look reveals that the derivation in 30,31 started from the tight-binding form of the spectrum which implies the limit U → ∞ first and then projecting the Hamiltonian onto the lowest level. In our situation all quantities, A, U and V , are defined by the same f and hence are of the same order, so that A cannot be eliminated.

E. Orthogonality and second quantization
For each Ψ n with the energy n , we observe that the Let us define the matrix Υ n = u n , v * n v n , u * n , obey-ingĤΥ n = n σ 3 Υ n σ 3 . We choose the normalization d 2 r |u n | 2 − |v n | 2 = 1, which leads to d 2 r Υ † nĤ Υ n = n 1. It can be further shown that the orthogonality condition holds, d 2 r (u * m u n − v * m v n ) = δ mn , which indicates the possibility of expansion with some complex valued coefficients, c n . It follows that the Hamiltonian in Eq. (14) takes the form with the latter expression acquiring familiar form n n (c † n c n + 1 2 ) after promoting the coefficients c n into second quantization operators.
We can write the orthogonality relations in the form which allows us to write the Green's function The sums over n in the formulas above are performed for all n > 0.

F. Zero modes
Two modes possessing zero energy can be constructed as follows. The classical energy (8) is obviously invariant upon the (infinitesimal) shift of coordinates, z → z + δz,z →z + δz, the function f transforms upon this as f (z,z) → f (z,z) + δz ∂ z f + δz ∂zf . It follows that the quadratic form (14) vanishes for Ψ ø ,Ψ ø of the form Notice that the form (29) ensures the zero energy even for those f 0 which do not obey the extremum condition (9). It means that in the latter case Ψ ø ,Ψ ø are orthogonal to the first variational derivative L 1 [f 0 ], see below.

A. Construction of skyrmion crystal
We showed previously 28,29 that the stereographic projection approach is an effective instrument for description of static skyrmionic configurations including SkX. Below we sketch the way to construct the trial function f , adopted for SkX configuration f SkX .
Our starting point is a simple form of f (r) in the Belavin-Polyakov (BP) solution of the model of planar ferromagnet with D = B = 0. The principal observation comes from (9) showing that any holomorphic or antiholomorphic function provides the local energy minimum in this case. The single skyrmion configuration of topological charge Q = 1 is given, e.g., by f = z 0 /z, with |z 0 | being a skyrmion radius. The multiskyrmion configuration can be represented as a sum of individual stereographic functions, f = z j /(z − Z j ), with arbitrary residues, z j , and positions, Z j .
The inclusion of DM interaction and b = 0 into the model, Eq. (8), obviously complicates the problem. However, this inclusion has continuous character and leads to the smooth deformation of the BP solution, whereas topological properties such as Q should not change. Our choice b > 0 leads to a decrease f → 0 at the infinity, i.e. we should restore the function of the type f = z 0 /z in the BP limit. The DM interaction defines a helicity of skyrmions (Néel or Bloch type), which means that z 0 should have a certain phase. We propose the ansatz for a single BP-like skyrmion where the shape of skyrmion is defined by the real-valued function κ of distance, r = √ zz, from the center of skyrmion; we fix κ(0) = 1. The real-valued z 0 gives the size of skyrmion, and the factor i in the numerator corresponds to the Bloch type of skyrmion stemming from (1), (9). The substitution of (30) into (9), gives the equation for κ in the form: with x = zz/z 2 0 , κ (x) = dκ(x)/dx etc. This nonlinear equation can be solved numerically for given b, z 0 .
It was shown 28 that the ansatz (30) allows to take into account the interaction between skyrmions and find optimal parameters of SkX. We model the SkX by the combination with a 1 = (0, a) and a 2 = (− √ 3a/2, a/2), where a is a cell parameter of SkX. Determining the profile of individual skyrmion for given b and some z 0 from Eq. (31), we place the skyrmions on the triangular lattice with parameter a. Calculating the average energy density from (8), we then find the optimal values of z 0 and a. We show the results of this calculation in Fig. 1. Interestingly, the optimal radius of skyrmions approximately follows the relation z 0 ≈ 1.5/b in the most relevant region, b ∈ (0.3, 0.7), which makes the profile function nearly gaussian one, κ(x) exp(−x/2b). 28 The primitive vectors of the reciprocal lattice space are b 1 = 2π a (1/

B. Correction to the trial ground state
A legitimate question arises, what is the consequence of the approximate character of our trial ground state function (32) which does not satisfy the equation (9). To answer this question, we return to (10) and represent the static part of it in the form where we use Eq. (9) and definẽ We see that the action starts from the quadratic terms of ψ , ψ new , which means the correction towards the true solution of (9) is −Ĥ −1L 1 . Next we expandL 1 in a series with the complete set of eigenfunctions of (20) with n > 0. The correction to the ground state function is then and the correction to the ground state energy takes the form We use these expressions below in the discussion of our results.

A. Bloch waves
The spatially periodic function f SkX (32) is the classical static background for fluctuations, which are found in terms of normal modes of (20). The translational symmetry of f SkX results in the same symmetry of A, U and V in equation (20). We look for the normal modes as quasi-periodic solutions where k lies in the first Brillouin zone, index Q runs over the reciprocal lattice, i.e. Q = nb 1 + mb 2 with n, m integer. The sum over Q in (38) corresponds to the Bloch function, V k (r), periodic in r-space. We find the Fourier transform for all potentials in Eq. (20) according to where v = √ 3a 2 /2 is the unit cell area. It is convenient to define a new potential, U = |A| 2 + U . Substituting Eq. (38) into Eq. (20) we obtain after some algebra After regularization, described in Appendix B, the co-efficients U Q , A Q and V Q decrease sufficiently with in-creasing |Q|. It allows us to consider a finite basis in reciprocal space for (40). We choose a simple way to restrict the basis set, Q = nb 1 + mb 2 with n, m = −N, . . . , N , thus considering a large rhombus in reciprocal space. Although the symmetry of such rhombus is lower than hexagonal, the difference becomes inessential for the large N as checked by examining the obtained Bloch functions. The precision of calculation depends on the basis size N b = 2 × (2N + 1) 2 and regularization parameter σ, and exact results would correspond to N → ∞ and σ → 0. The data reported below are obtained for N = 20 and σ = 0.05. In view of relative smallness of coefficients VQ in (40) the upper and lower components of "spinor" The results of numerical calculations for the band structure with different magnetic field b is shown in Fig.3. We describe these findings in the next subsections.

B. Tight binding model
The low-energy band structure depicted in Fig. 3 and also in Fig.5b below is remarkable, as there is no hint of parabolic dispersion expected for uniform ferromagnetic ground state, f 0 ≡ 0. At the same time one sees the emergence of almost flat bands. We tried to fit each band with the tight-binding form of dispersion, including nearest neighbor (t 1 ) and next-to-nearest neighbor (t 2 ) hopping on triangular SkX latticê where t 1 (k) = 2t 1 cos ak y + 2 cos 1 2 ak y cos We fitted the data for several low-lying modes, each with its own set of t j , and found rather satisfactory agreement. The results of this fit are shown in Fig. 4 and can be summarized as follows.
(i) All curves can largely be fitted by only two parameters, t 0 , t 1 . Setting t 2 = 0 worsens the fit insignificantly.
(ii) The curves of "bonding" (t 1 < 0) or "antibonding" (t 1 > 0) character evolve slowly with increasing b. They do not cross each other, as there is anti-crossing property.
(iii) The flat bands of "nonbonding" character (t 1 0) evolve rapidly with t 0 increasing with b. Upon this they pass through other bands, as can be seen in Fig. 4a (blue, orange, cyan curves).
(iv) The finite size of our basis apparently affects two aspects of our data in Fig. 3. First, we see the anticrossing feature, i.e. repulsion of non-bonding and bonding bands, which should presumably be absent. Second, the energy of the lowest band at Γ point is finite instead of expected zero, as discussed below.

C. Classification of the branches
In view of the evolution of the flat bands with b, and their possible crossing, it is helpful to identify the bands according to the behavior of their wave-functions near the center of skyrmions. This identification is most easily done in terms of angular momentum, which provides the connection of SkX excitations to single skyrmion distortions.
Near the skyrmion center, r = 0 the function can be presented as f (z,z) e iφ F (r). Time-dependent fluctuations, f → f + (1 + ff )ψ, describe deformations of skyrmionic shape characterized by angular dependence of ψ. Expanding near the origin, ψ(r) = m exp(imφ)a m (r), we have amplitudes a m (r) corresponding to different partial contributions to the wave function. These coefficients can be calculated at Γ point (k = 0) with the obtained coefficients for the Bloch function (38) as a m (r) = Q exp(−imφ Q )J m (Qr)C (1) Q , with Q = Q(cos φ Q , sin φ Q ), J m (Qr) the Bessel function and we take r a . It turns out that for each branch of the spectrum one of the coefficients a m (r) for some m is significantly larger than the others. Thus we can classify the branches by the corresponding numbers m. In Fig. 3 and Fig. 5 we mark the branches by different colors, with numbers m = 0, . . . 6 corresponding to green, yellow, black, red, blue, orange and cyan, respectively. The modes with m < 0 do not occur at low energies and are not shown in these Figures. Different m are related to different deformations of skyrmions. Consider the equation |f (z,z)| = 1, defining the curve on the plane, where the local magnetization lies in-plane, S 3 = 0. For solitary skyrmion such curve is a circle, whose radius is associated with the skyrmion radius. For SkX the equation |f (z,z)| = 1 yields a periodic structure of disconnected contours whose shape depends on ψ.
The shape of the deformed skyrmion is thus given by |f + (1 + ff )ψ| = 1. Due to the form of f , this equation can roughly be presented as A(r) + B(r) cos(m − 1)φ = 1 with m the above character of ψ. The mode with m = 1 does not change the shape of skyrmion but its radius only, it is hence the breathing mode. For given m the number |m − 1| defines the symmetry of skyrmion deformation: so that for m = −1, 3 it is a second order axis,  The lowest branch shown in Fig. 3 has the angular number m = 2 and corresponds to the zero mode discussed above in Eq. (29). We discuss it in more detail in the next subsection.
Since the dynamics of ψ k is simply given by exp(iω k t), we see that the shape of the skyrmion depends on the combination (m − 1)φ + ω k t with ω k > 0. It is clear then that the mode with m = 2 can be associated with clockwise rotation and the mode with m = 0 with counterclockwise rotation, in accordance with previous findings, see e.g. 22 .
D. The lowest lying mode and correction to the ground state The mode Ψ 1 (k = 0) with the minimum energy, black line (m = 2) in Fig. 3, attains nearly zero energy at Γ point. It can be compared with the analytic expression (29). Since f 0 ∼ 1/z, the behavior of Ψ ø in (29) is consistent with m = 2. Normalizing Ψ ø and Ψ 1 , we find almost perfect overlap between these two modes, e.g. Ψ ø |Ψ 1 (k = 0) = 0.993 for b = 0.52. The imperfect matching is due to the fact, that the main component of the true zero mode (29) is not continuous at the skyrmion positions, namely, Ψ ø ∼ z/z 0 .
Concerning the correction to the ground state due to imprecise character of our trial function, discussed in Sec. III B, we expanded the first variation (34) in our basis, according to (35). We found the small overlap ofL 1 only with the states with m = 1, particularly with the yellow branch in Fig. 3. Using Eq. (37) we evaluate the energy density correction as δE 0 = 7.1 · 10 −4 per unit cell, i.e.  negligibly small value.

A. Link-variable method
The nontrivial topological ground state leads to the appearance of the gauge potential A in Eq. (15). As a result, the emerging band structure has unusual properties, characterized by Berry curvature and Chern numbers of the bands. Nontrivial topology of the bands can leads to interesting physical phenomena, such as heat Hall conductance 24 . The topological properties of exci-tations in SkX were discussed previously 21 and the advantage of the stereographic projection method used here is to provide a better accuracy and details in description of these properties.
One can calculate Berry curvature and Chern numbers, having analytical expressions for wave functions of exact band as a function of wave vector k. In our case we numerically calculate eigenvectors of hamiltonian for any exact value of k and use a special method for discretized Brillouin zone 27 , as described below.
We consider Bloch state Ψ nk = e ikr V nk (r) referring to nth band and assuming that it is a smooth function of k. We use the expression for Berry connection with ∂ µ = ∂/∂k µ and The Berry curvature is Since in 2D system the only non-trivial component Ω n,µν is Ω n,12 , we write below simply Ω n ≡ Ω n,12 . The Chern number for nth band is obtained after the integration over BZ : In numerics, V nk (r) is not analytically known and even the smoothness of V nk (r) is not guaranteed. There are different approaches to calculation of C n on some mesh over BZ and we use the method by Fukui et al. 27 .
To shorten the notation, we write V nk (r) ≡ |n, k and consider a discrete mesh over BZ: k = i δk 1 + j δk 2 , with δk µ = b µ /N and indices i, j running from 0 to N .
For small enough mesh we have an expansion and for normalized functions, n, k|n, k = 1, we can write ln n, k|n, k + δk µ ≈ n, k| ∂ µ |n, k δk µ , so that the Berry connection is A n,µ (k) ≈ 1 δk µ Im ln n, k|n, k + δk µ .
Since A n,µ (k) is real-valued, and the absolute value of n, k|n, k + δk µ can be ignored, we define a new link variable as : U n,µ (k) = n, k|n, k + δk µ | n, k|n, k + δk µ | , Within the same accuracy, we have for the Berry connection : so that the Berry curvature can be written as plaquette combination and when integrated over BZ, becomes the Chern number C n ≈ i 2π BZ ln U n,1 (k + δk 2 )U n,2 (k) U n,1 (k)U n,2 (k + δk 1 ) .
We describe the results of application of this method below.

B. Results
In order to avoid the ambiguity with identification of the intersecting bands, we chose the value of the field b = 0.52, with the dispersion of eight lowest-lying bands shown in Fig. 5b. The topological properties of these bands are shown in Fig. 5c, characterized by the Berry curvature and the Chern number. It can be seen that the flat bands, n = 4 and n = 7 are topologically trivial, with Ω n (k) 0 for any k. This is not so for the lowest bands with pronounced dispersion, n = 1 and n = 2, with sizable Ω n (k), which is however zero on the average, with the Chern numbers C 1 = C 2 = 0. Here and below the integer values of C j are confirmed to the accuracy 5 · 10 −3 . The behavior of the Berry curvature of the remaining bands, n = 3, 5, 6, 8, is rather interesting. It is characterized by peaks at the symmetry points, Γ , K, K and nearly constant background within the whole BZ.
The peak in Ω 3 (k) at Γ with positive amplitude provides roughly half of the integral weight, C 3 = 1, the rest coming from the smooth background. The positive peak in Ω 3 (k) finds its negative counterpart of the same amplitude in Ω 5 (k) at Γ point, but the curvature of the 5th band, Ω 5 (k), is peaked with positive amplitude also at the points K, K so the overall weight is again positive, C 5 = 1.
Further, the 5th band is close in energy to the 6th band at points K and K , which is manifested by the negative amplitudes Ω 6 (k) there. This is however accompanied by large positive contribution Ω 6 (k) at Γ point and C 6 = 1 again. Finally, the large negative contribution Ω 8 (k) at Γ point is compensated by rather smooth positive background, with C 8 = 1.
The analysis of the bands with higher energies is more complicated, because it is hard to find the values of b when the intersection of the bands is avoided. This intersection may produce ambiguities in the definition (49).
Finally, we do not observe the closing and reopening of the gap between the bands, reported elsewhere for Néeltype skyrmions 23 for b 0.9. If this phenomenon would take place in our situation, this value of the field would correspond to the instability region, when the SkX phase is less favorable in energy than the uniform ferromagnetic state. 28

VI. CONCLUSION
We demonstrated the application of the stereographic projection method to the study of the dynamics of SkX in 2D system with DM interaction and in perpendicular magnetic field. The appearing band structure reveals a tight-binding character of the spectrum, with two nearly non-hybridized sets of bands, behaving differently with the increase of the magnetic field. One set of bands consists of flat bands, is topologically trivial, and rapidly evolves with the field. Another set is characterized by pronounced dispersion and non-trivial topological prop-erties. The Berry curvature may be sizable and signreversal in the Brillouin zone, with maximum values near symmetry points, whereas the integral Chern number is always non-negative for several lowest bands. These properties should result in existence of an edge states and unusual heat transport phenomena. Particularly, possible consequences of the sign-reversal character of Berry curvature are yet to be explored.
The developed theory describes the dynamics of SkX in terms of Schrödinger type equation for Bogoliubov spinors. The form of the potentials in this equation is fully defined by the complex-valued function f used for description of the classical ground state. In our particular study we let f be the ansatz for hexagonal ordering of Bloch-type skyrmions. From a more general viewpoint, the magnetic interactions of other types will lead to different topological structures, e.g. Néel-type skyrmions, antiskyrmions, quadratic ordering or skyrmion liquids. In these cases the form of the stereographic function f will also be different, however, the equation for dynamics should be mainly intact, with obvious changes taking into account additional interactions.