Magnon spectrum in two- and three- dimensional skyrmion crystals

TL;DR

This paper analyzes the low-energy magnon spectrum of skyrmion crystals in 2D and 3D, revealing stable localized states, their hybridization into a tight-binding spectrum, and the effects of magnetic field and temperature on skyrmion stability.

Contribution

It introduces a disk approximation for the SkX lattice, computes the magnon spectrum, and discusses quantum corrections and melting phenomena in 2D and 3D skyrmion crystals.

Findings

Localized magnon states are stable at certain fields.

Hybridization leads to a tight-binding magnon spectrum.

Quantum corrections can favor skyrmion ground state.

Abstract

We study the low-energy magnon spectrum of the skyrmion crystal (SkX) ground state, appearing in two-dimensional ferromagnet with Dzyaloshinskii-Moriya interaction and magnetic field. We approximate SkX hexagonal superlattice by a set of overlapping disks, and find the lattice period by minimizing the classical energy density. The determined spectrum of magnons on the disc of optimal radius is stable and only two lowest energy levels can be considered as localized. The subsequent hybridization of these levels in the SkX lattice leads to tight-binding spectrum. The localized character of the lowest magnon states is lost at small and at high fields, which is interpreted as melting of SkX. The classical energy of SkX is slightly above the energy of a single conical spiral, and a consideration of quantum corrections can favor the skyrmion ground state. Extending our analysis to three-dimensional case, we argue that these quantum corrections become more important at finite temperatures, when the average spin value is decreased.

Figures

Equations

$$H=\sum _{⟨ij⟩}\left[J_{ij}{\hat{\mathbf{S}}}_i{\hat{\mathbf{S}}}_j+{\mathbf{D}}_{ij}\cdot {\hat{\mathbf{S}}}_i\times {\hat{\mathbf{S}}}_j\right]-\mathbf{B}\sum _i{\hat{\mathbf{S}}}_i$$

$$E_{cl}=s^2\int d\mathbf{r}\left(\frac{C}{2}{\left(\nabla \mathbf{m}\right)}^2+D\mathbf{m}\left[\nabla \times \mathbf{m}\right]-s^{-1}\mathbf{B}\mathbf{m}\right),$$

$${\tilde{S}}_j^z$$

$${\tilde{S}}_j^{-}$$

$$\mathbf{m}=s\left(\begin{array}{ccc}\sin \beta \cos \alpha ,& \sin \beta \sin \alpha ,& \cos \beta \end{array}\right)$$

$$Q$$

$$\displaystyle=\frac{1}{4\pi}\alpha\left(\varphi\right)\big{|}_{\varphi=0}^{\varphi=2\pi}\cdot\cos\beta\left(r\right)\big{|}_{r=0}^{R}$$

$$\alpha \left(\phi \right)=Q\phi +{\alpha }_0,\beta (0)=\pi ,\beta (R)=0.$$

$${\left(\nabla \mathbf{m}\right)}^2$$

$$\mathbf{m}\left[\nabla \times \mathbf{m}\right]$$

$$\mathbf{B}\mathbf{m}$$

$$D\mathbf{m}\left[\nabla \times \mathbf{m}\right]=D\sin \left({\alpha }_0+(Q-1)\phi \right)\left(\frac{d\beta }{dr}+Q\frac{\sin 2\beta }{2r}\right)$$

$${\rho }_c$$

$$+\left\{r\frac{d\beta }{dr}+\frac{\sin 2\beta }{2}\right\}-br\left(\cos \beta -1\right))$$

$$\frac{d^2\beta }{d{r}^2}+\frac{1}{r}\frac{d\beta }{dr}-\frac{\sin \beta \cos \beta }{r^2}+\frac{2{\sin }^2\beta }{r}-b\sin \beta =0,$$

$$H_2=\int d\mathbf{r}\left(a_{\mathbf{r}}^{†}\hat{F}({\hat{L}}_z)a_{\mathbf{r}}+\frac{1}{2}\left(G{a}_{\mathbf{r}}^{†}{a}_{\mathbf{r}}^{†}+H.c.\right)\right)$$

$$F({\hat{L}}_z)$$

$$+\left\{-\frac{3\sin 2\beta }{2r}+\frac{2\sin \beta }{r}{\hat{L}}_z-\frac{d\beta }{dr}\right\}+b\cos \beta$$

$$G$$

$$a_{\mathbf{r}}=\frac{1}{\sqrt{2\pi }}\sum _m\sum _{n\ge 0}{e}^{im\phi }{f}_n^{(m)}(r)a_{m,n},$$

$$\left[a_{m,n},a_{m^{\prime },n^{\prime }}^{†}\right]={\delta }_{m{m}^{\prime }}{\delta }_{n{n}^{\prime }}$$

$$a_{m,n}$$

$$\displaystyle\sum_{k}\big{(}u^{(m)}_{nk}u^{(m)}_{n^{\prime}k}-v^{(m)}_{nk}v^{(m)}_{n^{\prime}k}\big{)}$$

$$\displaystyle\sum_{k}\big{(}u^{(m)}_{nk}v^{(-m)}_{n^{\prime}k}-v^{(m)}_{nk}u^{(-m)}_{n^{\prime}k}\big{)}$$

$$\sum _{k\ge 0}^{\infty }\left(\begin{array}{cc}{F}_{nk}(m)& G_{nk}\\ -G_{kn}& -F_{nk}(-m)\end{array}\right)\left(\begin{array}{c}{u}_k^{(m)}\\ v_k^{(-m)}\end{array}\right)={\epsilon }_{m,n}^0\left(\begin{array}{c}{u}_n^{(m)}\\ v_n^{(-m)}\end{array}\right)$$

$$a_{\mathbf{r}}=\frac{1}{\sqrt{2\pi }}\sum _m\sum _{k\ge 0}{e}^{im\phi }({\varphi }_k^{(m)}(r)c_{m,k}+{\tilde{\varphi }}_k^{(m)}(r)c_{-m,k}^{†}),$$

$${\varphi }_k^{(m)}(r)=\sum _n{f}_n^{(m)}(r)u_{nk}^{(m)},{\tilde{\varphi }}_k^{(m)}(r)=\sum _n{f}_n^{(m)}(r)v_{nk}^{(m)},$$

$$\left(\begin{array}{cc}F(m)-{\epsilon }_{m,n}^0,& G\\ G,& F(-m)+{\epsilon }_{m,n}^0\end{array}\right)\left(\begin{array}{c}{\varphi }_n^{(m)}\\ {\tilde{\varphi }}_n^{(-m)}\end{array}\right)=0.$$

$$H_2=\frac{1}{2}\int d\mathbf{r}(a_{\mathbf{r}}^{†},a_{\mathbf{r}})\left(\begin{array}{cc}\hat{F},& G\\ G,& \hat{F}\end{array}\right)\left(\begin{array}{c}{a}_{\mathbf{r}}\\ a_{\mathbf{r}}^{†}\end{array}\right)$$

$$E_q=\frac{1}{2}\left(\sum _{m,n}{\epsilon }_{m,n}^0-\sum _m\text{mbox}Tr{F}_m\right).$$

$${\rho }_{sp}$$

$$=0,b\ge 1$$

$$⟨a_{\mathbf{r}}^{†}{a}_{\mathbf{r}}⟩=\sum _{mk}{\left|{\tilde{\varphi }}_k^{(m)}(r)\right|}^2$$

$${\epsilon }_m(\mathbf{k})$$

$$\gamma (\mathbf{k})$$

$$(-{\partial }_r^2+\hat{M}-{\tau }_3{\epsilon }_{m,0}){\tilde{\psi }}_{m,0}=0$$

Full text

Introduction

Topological properties of the condensed matter systems is under active investigations over the last decade. Topological objects in magnetically ordered systems include various exotic spin structures, which are interesting both theoretically and experimentally, see J. (2016) and references therein. One example is a so-called skyrmion crystal, which is a regular array of magnetic vortices. The possibility of skyrmion textures has been envisioned in the earlier works.Belavin and Polyakov (1975); Bogdanov and Hubert (1994a); Ivanov (1995); Roessler et al. (2006)

The existence of skyrmion lattices was experimentally proven recently, in particular in the magnetic B20 compounds. Tchoe and Han (2012); Mühlbauer et al. (2009); Yu et al. (2010a, b); Grigoriev et al. (2007) Nowadays materials with confirmed skyrmion lattice include, e.g., metals $\mathrm{Mn\,Si}$ and $\mathrm{Fe\,Ge}$, semiconductors $\mathrm{Fe_{1-x}\,Co_{x}\,Si}$ Münzer et al. (2010); Yokouchi et al. (2014), and insulators $\mathrm{Cu_{2}\,O\,Se\,O_{3}}$.Nagaosa and Tokura (2013)

One of the specific property of one skyrmion (Fig. 1) is its statical stability which is provided by topological protection and by the fact that this configuration of spins corresponds to a local minimum of the classical energy. Belavin and Polyakov (1975)

Topological protection of skyrmion configurations is interesting not only on the fundamental reasons, but also from a technical point of view. It is expected, that one can create a super-dense, long-term storage media Kiselev et al. (2011); Fert et al. (2013) and transistors Zhang et al. (2015) based on skyrmions. In addition, a skyrmion lattice strongly affects spin current. An electron moving through a skyrmion lattice changes its spin orientation multiply, in order to adjust it to a local magnetization vector. The skyrmion lattice causes an effective force, which changes an electron direction of motion that should manifest itself macroscopically as a kind of Hall effect. Ochoa et al. (2016)

The static stability of skyrmion lattice was analyzed in Bogdanov and Hubert (1994a) in the classical limit, which limit is a usual theoretical approximation. Similar approach was applied to various magnetic vortex structures Metlov and Guslienko (2002); Metlov (2010) and to arrays of magnetic dots. Galkin et al. (2006) The semi-classical quantization method can be further used for the analysis of magnon spectrum, which was done for single skyrmion configuration in Schütte and Garst (2014); Aristov et al. (2015); Aristov and Matveeva (2016). Various aspects of dynamics in skyrmionic vortex structures were studied in Ivanov and Zaspel (2002); Butenko et al. (2010); Dai et al. (2013); Ivanov (1995); Ivanov et al. (1999).

In this work we investigate the low-energy spin-wave spectrum in the hexagonal skyrmion crystal. Our analysis is done in three steps. We start by considering the two-dimensional (2D) magnet with exchange interaction, Dzyaloshinskii-Moriya interaction placed in the uniform magnetic field. We study the system in classical limit and minimize its energy thus finding an optimal period of skyrmion structure. A hexagonal unit cell of a skyrmion lattice is approximated by a disc and the semi-classical quantization is used to obtain magnon spectrum on it. The energies and wave functions for magnons are found from the Schrödinger equation solved on the disc numerically. This solution implies that we consider the hexagonal cells as isolated, by effectively applying infinite magnetic field on the borders, so that the magnon from one cell cannot pass to another. We relax this effective field in the further treatment by invoking the idea of small tunneling of low-energy magnons between cells. The amplitude of this tunneling is found from the form of the wave-function on the border, and we come to his the tight binding form of the spectrum Abrikosov (1988). These findings are comparable to other studies of magnon spectrum on a skyrmion lattice. Mochizuki (2012)

In addition to that we also discuss effects connected to quantum nature of spin. First is the zero-point motion of spins arising in the non-parallel spin configuration. This motion contributes to the ground state energy of the system, although this quantum correction is formally small by inverse value of spin, $s^{-1}$. We notice however, that the purely classical energy of skyrmion lattice is slightly higher than the energy of single conical spiral, and the quantum correction is needed to avoid metastability of skyrmionic state. At finite temperatures, $T$, the equilibrium value of spin decreases and the role of these quantum corrections in the stabilization of the skyrmion state should increase.

Another effect is the quantum reduction of equilibrium spin value, which includes the zero-point contribution and the thermal correction. The zero-point contribution turns out to be small and non-uniform. The consideration of finite $T$ is impossible in purely 2D case, as the spectrum is gapless and the fluctuations destroy the long-range magnetic order. The inclusion of third spatial direction into our analysis is not difficult, when we recall that the experimental evidence in B20 compounds $\mathrm{Mn\,Si}$ and $\mathrm{Fe_{1-x}\,Co_{x}\,Si}$ shows the 3D skyrmionic structures as stacks of two-dimensional lattices. Within our model we show that the low-lying modes do not contribute much to the thermal reduction of spin value, due to the small volume of the Brillouin zone in the large-period skyrmion crystal.

The plan of the paper is as follows. We introduce the microscopic model and discuss its continuum limit in Sec. I. The basic ingredients of our analysis are introduced in Sec. I.1 and the determination of the superlattice period of skyrmion crystal is discussed in Sec. I.2. The magnon spectrum on the disc with one skyrmion is discussed in Sec. II. Particularly we analyze the dependence of the spectrum on radius of the disc in Sec. II.1, quantum corrections to the ground state energy in Sec. II.2 and the quantum reduction of spin in Sec. II.3. The band structure of magnons in the SkX is evaluated in Sec. III based on tight-binding model, the consistency of the model is checked here. The extension of our analysis to 3D case is done in Sec. III.1. We present our conclusions in Sec. IV. The technical details of our derivation are given in three Appendices.

I Skyrmion lattice model

I.1 2D chiral magnet and its classical description

We consider two-dimensional (2D) magnetic system without inversion center. In consideration we do not include anisotropy. Wilson et al. (2014) The model Hamiltonian on the square lattice is given by

[TABLE]

with $J_{ij}<0$ the ferromagnetic exchange, external magnetic field $\mathbf{B}$ is directed perpendicular to the 2D plane, and Dzyaloshinskii-Moriya (DM) interaction is characterized by the vector $\mathbf{D}_{ij}=D\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)$ is directed in the 2D plane. Crépieux and Lacroix (1998); Elhajal et al. (2002)

We introduce the local magnetization, $\mathbf{m}(\mathbf{r})=\langle\widehat{\mathbf{S}}_{\mathbf{r}}\rangle$, and assume the semiclassical limit, $s\gg 1$. From this point onwards it is set $\hbar=1$. The exchange $J$ is regarded as the largest interaction, so that adjacent spins are almost parallel, it is convenient to subtract this large energy of uniform ferromagnet, $s^{2}\sum_{ij}J_{ij}$, from the subsequent consideration. Making the gradient expansion of $\mathbf{m}(\mathbf{r})$, as described in Appendix A, we have the usual expression for the classical energy

[TABLE]

with $C\sim|J|$ being spin stiffness constant. When passing to continuum limit, i.e. $\sum_{i}\to a_{0}^{-2}\int d\mathbf{r}$, we introduce the dimension of length into the quantities $D\to D/a_{0}$ and $|\mathbf{B}|=B\to B/a_{0}^{2}$. In what follows we use dimensionless units. The energy is measured in units of $sC$ and the distance in units of $a_{0}C/D$. The classical energy then has the large prefactor $s$ as compared to quantum Hamiltonian below without this factor. The remaining dimensionless parameter $b$ equals $BC/D^{2}s$. In the discussion of the spectra below, we take $b=0.6$ because it produces the results, typical for the general case.

One can find the static configuration of $\mathbf{m}(\mathbf{r})$ producing the energy minimum, and then also determine the equation of motion of fluctuations around it. The spectrum of these fluctuations corresponds to conventional spin waves. However, if we want to retain a possibility to study effects of interaction between spin-waves, Aristov and Matveeva (2016) it is better to adopt another method described below. Upon this the classical configuration and the spectrum of linear spin wave theory remain unchanged.

Assume that the direction of the average magnetization $\mathbf{m}(\mathbf{r})$ changes with $\mathbf{r}$. We define the rotation matrix $\hat{U}\left(\mathbf{r}_{i}\right)$ at each site ${\mathbf{r}_{i}}$ so that $\mathbf{S}_{\mathbf{r}}=\hat{U}\left(\mathbf{r}\right)\tilde{\mathbf{S}}_{\mathbf{r}}$ and the average local spin in the new basis, $\tilde{\mathbf{S}}_{\mathbf{r}}$, is directed along the $\hat{z}$-axis. The position-dependent matrix is parametrized as $\hat{U}(\mathbf{r})=e^{-\alpha\sigma_{3}}e^{-\beta\sigma_{2}}e^{-\gamma\sigma_{3}}$ with generators of $SO(3)$ group $\sigma_{2}$,$\sigma_{3}$ and Euler angles $\alpha$, $\beta$, $\gamma$. We expect that the average spins in new local bases are directed along $\hat{z}$, and we use the Maleyev-Dyson representation for spin operators, preserving the spin commutation relations

[TABLE]

where $s$ is a value of spin, $\tilde{S}^{\pm}=\tilde{S}^{x}\pm i\tilde{S}^{y}$ and $[a_{j},a^{+}_{j}]=1$. In the classical limit, $s\to\infty$, we have $\tilde{\mathbf{S}}=s\hat{z}$ and

[TABLE]

Before we proceed further, we outline our computation strategy. It is known that in a certain range of magnetic fields the Skyrmion crystal (SkX) may be formed, with the hexagonal superlattice characterized by some lattice spacing which is denoted $R$. The average magnetization at the center of each hexagonal cell in perfect SkX is directed antiparallel to the field $\mathbf{B}$, and is parallel to it at the boundary of cells. To find a spectrum of magnetic excitations in such crystal we first determine an optimal $R$, denoted by $R_{0}$, from the minimization of the energy density. Next we find the magnon spectrum for elementary plaquette with the skyrmion at its center. Technically, these plaquettes in the form of discs Bogdanov and Hubert (1994a, b); Bogdanov and Yablonskii (1989) are taken and then the whole plane is paved by these overlapping discs. The overlap between different plaquettes leads us to the tight-binding model on the hexagonal superlattice. Calculating the parameters of the latter model, the band structure for magnons in the whole system is obtained.

I.2 Classical energy of skyrmion and optimal radius

Our first step is to approximate the hexagonal unit cell by a disc of radius $R$. (Fig. 2). The optimal radius of the disc is found by minimization of the average energy density for a single skyrmion $\rho_{c}=E_{sk}/(\pi R^{2})$ on that disc. By doing this we minimize the energy of a whole skyrmion lattice Han et al. (2010), which is approximatly equal to a sum of the energies of separate skyrmions (see below).

For the analysis of a single skyrmion placed at the center of the disc of radius $R$ we use polar coordinates $(x,y)=(r\cos\varphi,r\sin\varphi)$ and parametrize $\mathbf{m}$ by (4). Skyrmions can be characterized by the topological charge Rajaraman (1982):

[TABLE]

where the last line was obtained for centrosymmetric solutions $\alpha=\alpha(\phi)$ and $\beta=\beta(r)$.

Spin configuration on the plane can be viewed as mappings of one spherical surface onto another ($S^{2}\rightarrow S^{\prime 2}$) and can be classified into homotopy sectors. Rajaraman (1982) Mappings within one sector can be continuously deformed onto another and there is infinity number of such homotopy sectors or classes, characterized by integer $Q$. Non-trivial topological configuration of vector field corresponds to the boundary condition:

[TABLE]

Our solution corresponds to $Q=1$ and $\alpha_{0}=\pm\pi/2$ as is shown below.

The terms in Eq.  (2) can be represented in the form

[TABLE]

The exchange part of interaction, $\left(\nabla\mathbf{m}\right)^{2}$, and $\mathbf{B}\,\mathbf{m}$ are independent of the parameter $\alpha_{0}$. The DM part is related with topological charge and $\alpha_{0}$:

[TABLE]

The contribution of this term to classical energy for $D>0$ has a minimum at $Q=1$ and $\alpha_{0}=\pi/2$.

The classical energy (2) is now evaluated with the use of Eqs.  (6), (7), and dividing it by the area of disk $\pi R^{2}$ and by the factor $s$

[TABLE]

with the terms in curly brackets come from the DM part of interaction. For convenience we subtracted here the energy of uniform ferromagnet so that $\rho_{c}=0$ for $\beta\equiv 0$. The resulting Euler-Lagrange equation

[TABLE]

is supplemented by the boundary conditions from (6). The latter equation is not of hypergeometric type thus its solution cannot be expressed in a closed form through hypergeometric functions. However it is readily solved numerically for any particular $R$ by shooting method. Proceeding this way, we get the profile $\beta(r)$ and the dependence of $\rho_{c}$ on the disc radius, $R_{0}$. The results for the energy density are presented in the Fig. 3.

It can be seen in the Fig. 3 that in the range $b\in(0,0.8)$ there is a minimum in the energy density at some optimal radius, $R_{0}$. For magnetic field with larger values, $b\gtrsim 0.8$, the energy of skyrmion for any $R$ is higher than the energy of ferromagnetic ground state, $\beta\equiv 0$. We show below that the value of optimal $R_{0}$ determined from the classical energy can be further refined by consideration of quantum corrections to the ground state.

We also plot the dependence of the determined $R_{0}$ on the magnetic field $b$ in the Fig.  4. We obtain $R_{0}\sim 4$ almost in the whole range of $b$. This indicates that $R_{0}$ is mostly determined by the ratio $C/D$ in accordance with below estimates for the period of single spiral. Notice that $R_{0}$ mildly diverges at $b_{c}\simeq 0.8$, in accordance with Schütte and Garst (2014); Lin et al. (2014), presumably by logarithmic law.

II Magnon spectrum on the disc

The above minimization of the classical energy happens in order $s^{2}$. The obtained dependence $\beta(r)$ is used in the next terms of the expansion of the Hamiltonian in powers of $s^{-1}$. Each boson operator in this expansion is accompanied by the factor $s^{-1/2}$. The linear -in-bosons terms, $\sim s^{3/2}$, are absent as it should be for the extremum configuration in the semiclassical expansion. The linear-in- $s$ terms are quadratic in bosons, and they constitute the linear spin-wave theory (LSWT) Schütte and Garst (2014):

[TABLE]

with

[TABLE]

where $\hat{L}_{z}=-i\frac{\partial}{\partial\varphi}$ and $\nabla^{2}=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial r}\right)-\frac{\hat{L}_{z}^{2}}{r^{2}}$. Here again the curly brackets contain the DM contribution and the integration in (10) is done within the disc $r<R$.

The local magnon operators can be expanded in the basis functions $e^{im\varphi}f_{n}^{(m)}(r)$ forming an orthonormal set on the disk

[TABLE]

At the center of the disc we have $F(\hat{L}_{z}=m)=-r^{-1}\partial_{r}r\partial_{r}+(m+1)^{2}/r^{2}+{\cal O}(1)$ and $G={\cal O}(1)$. It is therefore convenient to choose the Bessel functions with shifted index, $f_{n}^{(m)}(r)\propto J_{m+1}(\kappa_{n}r)$, as discussed in Appendix B. Notice that the normalization of the wave-functions introduces the square of the inverse length scale into the energies. It means that magnon energies are obtained in units of $D^{2}/C$.

Operators $a_{m,n}$ can be further represented in terms of true magnon operators $c_{m,k}$ via Bogoliubov transformation

[TABLE]

where we require

[TABLE]

so that $[c_{m,k},c_{m^{\prime},k^{\prime}}^{\dagger}]=\delta_{mm^{\prime}}\delta_{kk^{\prime}}$.

In terms of these coefficients we come to the eigenvalue problem for each $m$

[TABLE]

where the expressions for $F_{nk}(m)$ and $G_{nk}$ are given in Appendix B.

Alternatively we may define

[TABLE]

which results in the equation

[TABLE]

II.1 Spectrum dependence on the disc radius

Using the method outlined above we computed the magnon spectrum on the disc for various $m$. To show a wider picture, we first relax the condition $R=R_{0}$ and depict the dependence of the energy levels on $R$ in the Fig. 5. For better visibility we show for each $m$ the lowest level for radial quantization.

It was shown in Aristov et al. (2015) that in the absence of DM interaction and of magnetic field the magnon spectrum for Belavin-Polyakov skyrmion possesses three zero modes, corresponding to conformal symmetries of the action. In the present notation these modes have $m=-1,0,1$, for translation, dilatation and special conformal (rotation at infinity) symmetry, respectively. The modes with $m=0$ and $m=1$ acquire the finite energy due to existence of DM interaction and magnetic field. The energy of the mode with $m=0$ is non-zero because the DM interaction defines the characteristic scale, $C/D$. The mode with $m=1$ has positive energy because the direction of magnetic field breaks the rotational symmetry in spin space.

The mode with $m=-1$ has non-zero energy because the translational invariance is lost on the disc. Its energy is comparable to other modes at optimal radius, $R_{0}\simeq 4$, and decays exponentially with $R$ cf. Schütte and Garst (2014). All energies tend to constant values at large $R$. The formal expression existing for zero mode $\psi_{m=-1,n=0}$ at $R\to\infty$ (see Eq. (51) inSchütte and Garst (2014)) is not inapplicable to our case of finite $R$ because $\beta^{\prime}(R)\neq 0$.

The obtained spectrum is always stable, i.e. all energies are positive, for optimal disc radius, $R_{0}$. This was checked for various $b$ and should be compared with the instability, discussed in Schütte and Garst (2014). Namely, negative energies were obtained for $m=-2$, $b\lesssim 0.56$ and for large $R$. The latter instability for a disc with non-optimal $R\neq R_{0}$ indicates that the energetically favorable spin configuration is not a single skyrmion ground state. Schütte and Garst (2014); Ezawa (2011) We show the calculated lowest energy for $m=-2$, $b=0.5$ and various $R$ in the Fig. 6. In this case $\epsilon_{-2,0}>0$ for optimal value $R=R_{0}=3.9$ and becomes negative at larger $R\gtrsim 0.56$.

II.2 Quantum correction to the ground state energy

The LSWT Hamiltonian (10) is normal-ordered in $a_{\mathbf{r}}^{\dagger}$, $a_{\mathbf{r}}$. The diagonalization of it in matrix form is done by first representing

[TABLE]

then expressing it in the basis (15) and subsequently finding the unitary Bogoliubov transformation to come to the form $H_{2}=\sum_{mn}\varepsilon^{0}_{m,n}c_{mn}^{{\dagger}}c_{mn}$, which is diagonal in terms of true magnon operators, $c_{mn}^{{\dagger}}$, $c_{mn}$. It is easily verified that these manipulations involve a loss of normal ordering in terms of $a_{\mathbf{r}}$, and returning to normal ordering in terms of $c_{mn}$. Two appearing commutators are not equal and their difference is the quantum correction, $E_{q}$, to the energy of the ground state. After some calculation we find

[TABLE]

This quantity is calculated on the disc of the radius $R$ using the formulas in Appendix B. The results are shown in Fig.  7 for the energy and its average density, $\rho_{q}=E_{q}/\pi R^{2}$. The dependence of density of quantum correction, $\rho_{q}$ on $R$ is rather pronounced and $\rho_{q}\to 0$ at large $R$. What is interesting is the existence of the minimum at $R\sim 1$, where the density of classical (ground state) energy is high, see Fig. 3.

It means that the above optimal radius, $R_{0}$, calculated purely in classical terms, does not correspond to the minimum of energy density at finite $s$. The consideration of quantum correction shifts the optimal $R_{0}$ to smaller values. This correction can be neglected in case $s\gg 1$, where our theory should work well. It is instructive to evaluate the total ground state energy, $E_{cl}+E_{q}$, for the ultra-quantum case $s=1/2$. As shown in the Fig. 8 the optimal radius in this latter case becomes $R^{uq}_{opt}\simeq 3.5<R_{0}\simeq 4$. Summarizing, we see that the consideration of the quantum correction to the ground state energy lowers the optimal values for the skyrmion crystal lattice spacing, but this shift in $R_{0}$ is not too big even for small $s$. Therefore we can safely use the previously determined spectrum, Fig. 5, which exhibits no big variation around $R=R_{0}\simeq 4$.

It can be argued Maleyev (2006) that the single conical spiral is a good candidate for the true ground state of the system, described by Eq. (2). As we show in Appendix C, the energy density of the conical spiral without the ferromagnetic contribution, $-B$, is given by the expression

[TABLE]

In the Fig. 9 we compare the energy density described by Eq. (19) with the density determined for the skyrmion configuration on the disc, Eq. (8). We see that in the whole range of the fields, $b\in(0,1)$ the energy of the conical spiral is slightly lower than that of the skyrmion configuration. This means that the skyrmion configuration should be considered as metastable. Although this fact does not prevent us from the determination of the magnon spectrum below, it may cast doubts in such procedure. We note here that the discussed quantum corrections to the ground state lowers the energy of the skyrmion configuration and can ultimately make it favorable. The validation of the latter statement requires however a calculation of the quantum correction to the spiral state, which is beyond the scope of the present study.

II.3 Quantum spin reduction

In addition to quantum correction to the ground state energy we can discuss the contribution of the zero point motion to the equilibrium spin value, $\langle\delta S^{z}_{\mathbf{r}}\rangle=s-\langle\tilde{S}^{z}_{\mathbf{r}}\rangle=\langle a^{\dagger}_{\mathbf{r}}a_{\mathbf{r}}\rangle$. It is well known that this correction in case of uniform 2D Heisenberg antiferromagnet on the square lattice is $\langle\delta S^{z}_{\mathbf{r}}\rangle\simeq 0.20$, because the staggered magnetization is not a conserved quantity. In our case this quantum spin reduction appears also because the sum of equilibrium vector spin values, $\sum_{\mathbf{r}}\mathbf{m}(\mathbf{r})\mathbf{S}_{\mathbf{r}}$, does not correspond to the conserved operator.

Using the above formulas (16) we represent the spin reduction in the form

[TABLE]

The obtained quantity does not depend on $s$ and is plotted in Fig. 10. It is seen that even for $s\sim 1$ the contribution of zero-point motion is negligible.

III Band structure of magnons on skyrmion lattice

Until now we discussed the spectrum of magnons confined on the disc with the skyrmion configuration of the classical state. The demand for the spin orientation at the edge be parallel to the field, $\beta(R)=0$, led to the property $\psi_{m,n}(R)=0$. The latter condition may be modeled by the quantum box potential $U_{edge}(r=R)=\infty$. When discussing the magnon spectrum for the skyrmion lattice we first pave the whole plane by the hexagons, which are separated from each other by the infinite edge potentials. The height of these potentials is then gradually reduced to zero, which results i) in the continuum spectrum for delocalized magnon states with $E>b$ and ii) in the tight-binding form of the spectrum for states with energies $E<b$.

Strictly speaking, the estimate for the boundary for the continuum spectrum, $E_{c}=b$, is obtained from (10) by putting there $\beta(R)=0$ and neglecting derivatives, $\beta^{\prime}(R)$, in which case the Hamiltonian is diagonal and corresponds to free motion with the dispersion $E_{k}=Ck^{2}+B$. It turns out that for the disc of optimal radius the derivatives, $\beta^{\prime}(R)$, are sizable and the picture becomes complicated as discussed below. For our estimates it is enough to take $E_{c}\simeq b$.

In the Fig. 5 we see that for optimal radius of the disc only two magnon states with $m=-1$ and $m=-2$ lie below the magnon continuum $E<E_{c}\simeq 0.6$. We show the evolution of these energy modes with field in the Fig. 11. The wave function of these states is mostly present at the center of the disc and is relatively small at its edge. When lowering the confining barrier between discs (hexagons) these states should remain almost localized we expect that the dispersion obtains the tight-binding form :

[TABLE]

with $\ell\simeq 1.94R_{0}$ is the distance between neighboring skyrmions, obtained from equality of the area of disc and hexagon, $\pi R_{0}^{2}=\sqrt{3}\,\ell^{2}/2$. The highest and the lowest band energies are $\varepsilon_{m,n=0}^{0}+6t_{m}$ and $\varepsilon_{m,n=0}^{0}-3t_{m}$ which is achieved at $\mathbf{k}=0$ and $\mathbf{k}=(4\pi/3\ell,0)$, respectively. Our aim now is to determine the on-site energy $\varepsilon_{m,n=0}^{0}$ and the hopping amplitude $t_{m}$. We estimate these quantities by evaluating the wave-functions at the edge.

Consider first a usual example of 1D symmetric two-well potential with $H=-\nabla^{2}+U(x)$, whose low-lying states $\epsilon_{i}$ lie below the barrier potential $U(0)$. In the semiclassical regime the wave-function of the individual well decreases exponentially at the barrier $\psi(x\simeq 0)\sim\exp(-\kappa x)$, with $\kappa^{2}\simeq U(0)-\epsilon_{i}$. The eigenfunctions are obtained as symmetrized combinations $\psi_{s}\sim\cosh(\kappa x)$, $\psi_{a}\sim\sinh(\kappa x)$. The solution $\psi_{a}$ has lower energy and the energy difference is $\delta E=E_{s}-E_{a}\simeq 2\psi_{s}(0)\psi_{a}^{\prime}(0)=2(\psi_{a}^{\prime}(0))^{2}/\kappa$. Landau and Lifshitz (1981)

In comparison with this usual case we have three modifications. First we have to pass from the 2D Hamiltonian, $r^{-1}\partial_{r}r\partial_{r}$, to 1D case, $\partial_{r}^{2}$ by making substitution $\psi(r)=r^{-1/2}\tilde{\psi}(r)$. Second, more importantly, we have a two-component spinor form of wave functions, rather than scalar one. As a result, we have two eigenfunctions $\psi_{1,2}$ at the edge, characterized by two values, $\kappa_{1,2}$. Third modification addresses the parity of the wave-function in the individual well. For odd functions the lower energy combination $\psi_{a}$ is given by the sum of individual functions. In our case it happens for odd $m$, and is accompanied by a change of sign of $t_{m}$ in (21).

We proceed by linearizing the equations (17) at $r=R$

[TABLE]

with $m=-1,-2$ and constant $\hat{M}={\hat{\cal H}}^{(m)}_{r=R}$. The above equation has general solution ${\tilde{\psi}}_{m,0}(r)=\sum_{i=1,2}c_{i}\tilde{\psi}_{i}^{0}\sinh\kappa_{i}(r-R)$, having the property ${\tilde{\psi}}_{m,0}(R)=0$. Knowing the eigenfunctions ${\psi}_{m,0}(r)$ of the full equation (17) we can extract the weights $c_{i}$ in ${\tilde{\psi}}_{m,0}(r)$ and verify that this function is a good approximation for the exact numerical $r^{1/2}{\psi}_{m,0}(r)$ at $r\simeq R$. Finally the energy splitting is found as $\delta E=2\sum_{i=1,2}c_{i}^{2}\kappa_{i}\langle\tilde{\psi}_{i}^{0}|\tilde{\psi}_{i}^{0}\rangle>0$.

Denote the lowest energy level as $E_{1}\equiv\varepsilon_{m=-1,n=0}$. In our approximation by the two-well potential the symmetrical solution has the energy $E_{s}=E_{1}+\delta E_{1}$ and hence the on-site energy is $E_{1}+\tfrac{1}{2}\delta E_{1}$ whereas the hopping is $t_{1}=-\delta E_{1}/2$, with the change of sign discussed above. It follows that the bottom of the lowest band on the hexagonal superlattice has the energy $E_{bottom}=E_{1}-\tfrac{5}{2}\delta E_{1}$, and the top of the band corresponds to $E_{top}=E_{1}+2\delta E_{1}$. Similar consideration holds for $m=-2$, but now $t_{2}=+\delta E_{2}/2$.

If our approximate calculation were exact, we would have $E_{bottom}=0$ for $m=-1$. This property is expected, since the finite energy $\varepsilon_{m=-1,n=0}$ corresponds to a broken translational symmetry on the disc, and consideration of the skyrmion lattice restores the translational symmetry and should result in the Goldstone mode.

We show the energies for the top and the bottom of the bands with $m=-1$, $m=-2$ in the Fig. 11. Both bands become narrow near critical field $b\simeq 0.8$, because the distance between the skyrmions increases and the overlap of the individual wave-functions decreases in this case. In this case the semiclassical treatment is justified and we see that $E_{bottom}$ touches zero for $m=-1$, as expected. However, $E_{bottom}$ becomes negative for smaller values of $b$, where the semiclassical regime fails; it provides the estimate for the accuracy of our calculation.

At lower fields, $b\lesssim 0.5$ the top of the lower bands increases to magnon continuum at $E\sim b$. In this case our initial picture of localized levels in individual wells and their eventual hybridization is hardly justified. This argument can be applied already to Fig. 11 with the lower bound on the field $b\simeq 0.35$, and our calculation of the bandwidths changes this bound to higher values, $b\sim 0.5\div 0.6$. It means that the skyrmion crystal can be considered as such in the range of the fields, $b\in(0.6,0.8)$. Similar estimate is found in Schütte and Garst (2014), however based on the other arguments.

III.1 3D case

The generalization of our analysis to the three-dimensional case is straightforward. We assume that the classical configuration of spins is given by 2D picture extended to third dimension. Instead of points on the plane, representing the centers of skyrmions, we now have the vortex lines in space. The classical treatment of Sec. I is unchanged. The consideration of Sec. II is modified as follows. The Schrödinger equation (10) acquires the term $-\nabla_{z}^{2}$ and is solved by separation of variables. The basis functions in (12) are $e^{ik_{z}z}e^{im\varphi}f_{n}^{(m)}(r)$ with the resulting change in the spectrum

[TABLE]

see Eq.(15). The conclusions of the rest of this section are unchanged, including the discussion around Eq. (19). The analysis on the lattice, Sec. III, now includes the modification of the spectrum (21) according to (22), which does not change the rest of this section.

What is important, is the possibility to discuss in 3D case the case of non-zero temperature without losing the long-range order. In purely 2D case the calculation of the spin reduction gives divergent quantity. We may estimate

[TABLE]

here $N(\omega)=(\exp(\omega/T)-1)^{-1}$, and we ignored the amplitude of wave-functions for simplicity. The singular contribution to averaged fluctuations comes from the Goldstone mode $m=-1$. At small $|\mathbf{k}|$ we have $\varepsilon_{-1}(\mathbf{k})\propto k^{2}$ in (21), which leads to logarithmic divergence at $T\neq 0$ and a loss of magnetic order in skyrmion lattice, in accordance with Mermin-Wagner theorem. The third spatial dimension improves the situation, as we may read from the Fig.  11 for intermediate fields

[TABLE]

Notice here that $\varepsilon_{-1,0}^{0}\sim D^{2}/C$ and $1-\gamma(\mathbf{k})\sim(k\ell)^{2}$ with $\ell\sim C/D$, when we restore dimensional units. The integration over in-plane $\mathbf{k}$ is done within the hexagonal Brillouin zone, $|\mathbf{k}|\lesssim\ell^{-1}$, related to skyrmion lattice with a period $\sim C/D$. Simple estimate gives for this singular contribution

[TABLE]

which is small except for $b\simeq 0.8$. The contribution from the bands with $m\neq-1$ can be estimated by approximate equating it to the case of uniform ferromagnet

[TABLE]

Summarizing here, we showed that the long-range order is stabilized in 3D and is unaffected by the existence of translational Goldstone mode of the skyrmion crystal everywhere except the close vicinity of the critical field.

IV Conclusions

We calculated the spectrum of magnons in the skyrmion crystal by first solving the eigenvalue problem on the disc of finite radius and then considering the hybridization of the discrete energy levels in the hexagonal skyrmion superlattice. Following the earlier works, we determine the superlattice spacing by minimization of classical energy density. The finite superlattice spacing leads to stability of the system at the level of one cell, and only two energy levels lie below the magnon continuum and may be considered as localized. The subsequent hybridization of these levels leads to the bands of finite width in the hexagonal tight-binding scheme. We show that the translational Goldstone mode is approximately restored in our numerical approach. The almost-localized character of the lowest energy levels is lost at smaller values of the field due to significant overlap of these lower bands with the higher-energy continuum. At higher fields $B\simeq 0.8D^{2}/C$ the superlattice period diverges and the quantum fluctuations destroy the long-range order. This may be interpreted as a disappearance of the skyrmion crystal at the fields $B$ outside the interval $(0.5\div 0.8)D^{2}/C$. The classical energy density of the skyrmion lattice lies slightly above the energy of the single conical spiral in the whole range of magnetic fields. The calculated quantum correction of order $1/s$ to the classical energy can favor the skyrmion crystal at higher temperatures, when the average value of spin is reduced. Our findings may indicate the mechanism for creation of skyrmion crystal in the so-called A-phase in B20 compounds observed at high temperatures in the restricted range of magnetic fields. Mühlbauer et al. (2009)

Acknowledgements.

We thank M. Garst, S.V. Maleyev, K.L. Metlov, A.O. Sorokin for useful discussions.

Appendix A From lattice model to continuum model

In this section we derive the continuum version of our model, suitable for possible subsequent analysis of the interaction between magnons. We use below Latin indices for spin space coordinates ($a,b,c,d=1,2,3$), and Greeks indices for physical space ($\alpha,\beta,\gamma=1,2$). We write the Hamiltonian (1) in the form

[TABLE]

with totally antisymmetrical tensor $\varepsilon_{abc}$ and magnetic field $B$ supplied by factor $s$ for further convenience. DM interaction has mixed spin and space indices stemming from its spin-orbital nature.

Let the ground state be characterized by non-collinear spin configuration. We express Eq. (24) in such local basis, where the average local spin, $\tilde{\mathbf{S}}_{\mathbf{r}}$, is directed along the $\hat{z}$-axis, with $\mathbf{S}_{\mathbf{r}}=\hat{U}\left(\mathbf{r}\right)\tilde{\mathbf{S}}_{\mathbf{r}}$. The position-dependent matrix, $\hat{U}(\mathbf{r})=e^{-\alpha\sigma_{3}}e^{-\beta\sigma_{2}}e^{-\gamma\sigma_{3}}$ is defined with generators of $SO(3)$ group $\sigma_{2}$, $\sigma_{3}$ and Euler angels $\alpha$, $\beta$, $\gamma$.

The Euler angle $\gamma$ is not determined by variational equation (9) and we may chose $\gamma=\alpha$ for the continuity of LSWT equation (10) at $r=0$, see Ref. Aristov et al. (2015). In terms of bosons the choice of $\gamma$ is encoded in simple $U(1)$ transformation, $a\to e^{-i\gamma}a$, and we adopt $\gamma=0$ in this paper, in order to make a better comparison with other authors Schütte and Garst (2014).

We assume the long wavelength limit, ${\mathbf{qn}}\ll 1$, which particularly corresponds to smooth variation of $\hat{U}\left(\mathbf{r}\right)$ on the scale of interatomic distances. In this case we can write

[TABLE]

It is convenient to define quantities

[TABLE]

with the explicit expressions for $U(\mathbf{r})$ and $\chi_{1}(\mathbf{r})$, $\chi_{2}(\mathbf{r})$ given in Aristov et al. (2015). We perform the summation over $\mathbf{n}$ according to the rules

[TABLE]

with $J(\mathbf{q})=\sum\limits_{\mathbf{n}}e^{i\mathbf{qn}}J(\mathbf{n})\simeq J(0)+\tfrac{1}{2}Cq^{2}\,.$ As a result we come to the Hamiltonian

[TABLE]

The classical part of the energy, Eq. (2), is obtained by putting $\tilde{\mathbf{S}}_{\mathbf{r}}=(0,0,s)$ in Eq. (26) and using the convention (4). Explicit formulas are given in Eq. (7).

The terms linear in bosons have the factor $s^{3/2}$ and vanish due to the extremum condition on $\beta$ (9) after integration by parts. Quadratic terms in bosons lead to Eq. (10).

Appendix B Numerical methods

The eigenfunctions to (17) are expanded in the basis of the orthonormal functions

[TABLE]

with $k=0,1,2\ldots$, $J_{m}(x)$ the Bessel function and $z_{k}^{(m)}$ its $(k+1)$-th positive zero. The set $\eta_{mk}=e^{i(m+1)\varphi}f_{k}^{(m)}$ satisfies the Laplace equation on a disc of radius $R_{0}$ and with boundary conditions

[TABLE]

with eigenvalues $\lambda_{k}^{(m+1)}=({z_{k}^{(m+1)}}/{R_{0}})^{2}$.

In the above basis the equation (17) acquires the matrix form

[TABLE]

with vectors $\varepsilon^{0}_{m,n}$ and $u_{n}^{(m)}$, $v_{n}^{(-m)}$ of, generally speaking, semi-infinite length. Practically, it suffices to take the basis of 40 lowest $z_{k}^{(m)}$’s to achieve a good accuracy for the low-energy spectrum. The matrix elements are given by

[TABLE]

with

[TABLE]

The spectrum of magnons is numerically obtained by diagonalizing the matrix (15). We calculate the spectrum for $m>0$ only, due to the property $\tau_{1}\tau_{3}{\widehat{H}^{(m)}}\tau_{1}=-{\widehat{H}^{(-m)}}$, with Pauli matrices $\tau_{i}$. The positive eigenvalues then give $\varepsilon^{0}_{m,n}$ for $m>0$ and negative eigenvalues correspond to $-\varepsilon^{0}_{m,n}$ for $m<0$.

Appendix C Energy of single spiral

In this section we evaluate the classical energy of the single conical spiral, whose form is given by the expression (4) with

[TABLE]

Here $\mathbf{qr}=q_{x}x+q_{y}y$, with vector $\mathbf{q}$ determined shortly. The expression (4) is written in orthogonal frame $\hat{e}_{1}$, $\hat{e}_{2}$, $\hat{e}_{3}$, with angle $\Theta$ between $\hat{e}_{3}$ and normal to the plane, $\hat{e}_{z}$, i.e.  $\hat{e}_{3}\dot{\hat{e}}_{z}=\cos\Theta$. Let the field $\mathbf{B}$ be directed at an angle to the plane, $\mathbf{B}\hat{e}_{z}=B\cos\Theta_{1}$, and to the center of the cone, $\mathbf{B}\hat{e}_{3}=B\cos\chi$. Simple calculation with the use of (2) gives for the energy density

[TABLE]

where $\mathbf{d}=D(\sin\Theta,-\cos\Theta,0)$ in the laboratory frame and we chose $\hat{e}_{3}\dot{\hat{e}}_{y}=0$ for definiteness. Clearly, the minimum of $\rho$ happens at $\mathbf{q}=\mathbf{d}/C$ and at $\chi=0$. The latter condition shows that the cone axis is directed along the field, so that $\Theta_{1}=\Theta$ ; we obtain

[TABLE]

Variation over $\beta$ gives $\cos\beta=BC/D^{2}=b$ for $b<1$ and $\cos\beta=1$ otherwise. Subtracting here the energy of uniform ferromagnet $\rho_{FM}=-B$, we obtain the expression (19) for the energy gain of the single spiral. Maleyev (2006)