Chiral standing spin waves in skyrmion lattice

The shape and size of magnetic samples define many of their static and dynamic properties. The most representative example is the so-called bubble materials^1,2—thin ferromagnetic films with strong perpendicular anisotropy. The size of the magnetic domains in these systems depends on the film thickness. The variation of the thickness changes the energy balance between demagnetizing fields effects, Heisenberg exchange, and magnetic anisotropy. The magnetic phases observed in the bulk crystals and thin films of bubble materials are very different. The ground state of very thin films with a thickness comparable to the domain wall width is the homogeneous ferromagnetic state, whereas, for thicker films, the multidomain states—stripes or bubbles—become more energetically favorable.

Another example of the systems where the shape of the samples plays a crucial role is the so-called chiral magnets, where the main properties are defined by the competition between the Heisenberg exchange interaction and chiral Dzyaloshinskii–Moriya interaction (DMI).^3,4 At zero external magnetic field, in bulk crystals of B20-type FeGe,^5–8 MnSi,^9–12 Fe_1−xCo_xSi,^13,14 and other compounds from that class,¹⁵ the ground state of the system is the helical spiral. With increasing external magnetic fields, the spiral continuously transforms into the cone phase and then into the field saturated state. These three phases are the main phases observed in the bulk crystals of isotropic chiral magnets in a whole range of temperatures. There are, however, a few exceptions; for instance, a tiny pocket near the critical temperature on the phase diagram of most of the bulk chiral magnets known as the anomalous phase,^9–11 and additional phases that are observed at low temperature in Cu₂OSeO₃.^16,17 The experimental observations in thin films of chiral magnets show the emergence of another phase—the skyrmion lattice, which remains stable in a wide range of magnetic fields.^{5–8,12–15} A theoretical explanation of this phenomenon is based on the effect known as the chiral surface twist,^18,19 which provides a significant skyrmion energy gain when the film thickness equals a few periods of helical modulations.^20,21 The energy gain provided by the chiral twist near the surface reduces with the film thickness, and above a critical thickness, the skyrmion phase becomes energetically unfavorable.²¹ 

It is natural to expect that besides the static properties, the system’s dimensionality and the edge effects also significantly affect the dynamic properties of the chiral magnets. For instance, it is known that in pure ferromagnetic films at particular boundary conditions,²² the spin-wave resonance^23,24 can take place. In chiral magnets, the presence of DMI causes an effective “pinning” of the spins at the boundaries.^18,25–28 The latter can lead to the formation of the standing spin-waves (SSW), where all the spins precess with identical frequency, but the relative phase and the amplitude of the spin precessions vary across the thickness. At the same time, the positions of the minima and maxima of the amplitude, nodes and antinodes, respectively, remain fixed in time. Analytical models based on the linear approximation suggest that SSW can arise in chiral magnets in both the field saturated state and the helical spin spiral state.^25,29–31 On the other hand, in some works, it was argued that SSW with fixed in time positions of the antinodes does not exist in chiral magnets³² or at least cannot be considered as conventional SSW.³³ One of the aims of this work is to confirm the presence of SSW in chiral magnets and to highlight some essential aspects that have not been discussed in the previous studies.

Here, we present the results of the theoretical study of the spin-wave resonance in the three-dimensional (3D) skyrmions lattice (SkL) representing a statically stable configuration in the plate of isotropic chiral magnets. We show that the spectra of resonance frequencies, in this case, differ from the 2D model, and the additional resonance modes can be explained in terms of chiral SSW.

We consider the classical spin model of the simple cubic lattice described by the following Hamiltonian:

where n_i = M_i/μ_s is a unit vector of the magnetic moment at site i, μ_s is the value of the magnetic moment of the lattice site, $ij$ means the summation of overall nearest–neighbor pairs, J is the exchange coupling constant and D_ij is the Dzyaloshinskii–Moriya vector defined as D_ij = Dr_ij, D is the DMI scalar constant, and r_ij is the unit vector between sites i and j. The external magnetic field, B(t) = B_DC + B_AC(t), contains the static (DC) and time-dependent (AC) components. To ensure the generality of the results, we use reduced units of the external magnetic fields and the distances with respect to the saturation field B_D = D²/(μ_sJ) and equilibrium period of helical modulations in the continuum limit, L_D = 2πaJ/D, respectively. While the values of B_D and L_D are unique for a particular chiral magnet, the results provided in reduced units are general for the whole class of magnets described by the model (1). Using the reverse conversion for a particular compound with B_D and L_D known from the experiment, one can obtain the absolute values of the external magnetic fields and the sizes of the sample required for the observation of the theory’s predicted effects.

We excluded the dipole–dipole interactions since they do not affect the discussed phenomena significantly (see  Appendix B).

We consider an extended film of finite thickness with periodic boundary conditions (PBC) in the xy − plane and free boundaries along the z-axis [Fig. 1(a)]. The shape and size of the simulated domain in the xy − plane is chosen to fit the equilibrium size of the unit cell of the hexagonal skyrmion lattice, $Lx/Ly=78/45≈3$⁠. The thickness of the film L_z varies by changing the number of atomic layers, d, from 1 to 140, which for a chosen coupling parameter corresponds to L_z = 0 (in the single atomic layer limit—2D case) and L_z = 4L_D. For definiteness, we set J = 1 and D = 0.18, which results in L_D ≈ 35a and $BD=0.0324μs−1$⁠. The static magnetic field is always perpendicular to the film plane, B_DC‖e_z. In the following, we study both cases of in-plane and out-of-plane alternate magnetic fields: B_AC ⊥B_DC and B_AC‖B_DC.

The magnetization dynamics is described by the Landau–Lifshitz–Gilbert (LLG) equation, which can be written in the dimensionless form

where t is a dimensionless time scaled by γE₀/μ_s, with γ being the gyromagnetic ratio and E₀ being the reference energy, which we set equal to exchange constant J since the Heisenberg exchange is the leading energy term in (1), α is the Gilbert damping, and $Hi=−E0−1∂E/∂ni$ is a dimensionless effective field on the ith lattice site. The solution of the LLG equation provides time dependent magnetization vector field n(t) that, in the case of small oscillations, can always be decomposed into static, n₀, and dynamic, δn(t), components: n(t) = n₀ + δn(t) [see Fig. 1(b)]. For visualization of the fields n₀ and δn, we use the color scheme explained in Fig. 1(c).

The calculation of the absorption spectra in Figs. 2(a) and 2(b) is based on Fourier analysis of the dynamic component of magnetization induced by an external magnetic field pulse (see  Appendix A for details). The absorption spectra in Figs. 2(a) and 2(b) can be thought of as the system response to the external magnetic field oscillating with the frequency ω. The picks of th imaginary component of complex dynamic susceptibility, Imχ(ω), correspond to those eigenfrequencies of the system, which can be excited in the skyrmion lattice by the uniform AC field.

First, we consider the case of B_AC ⊥B_DC. The dynamic susceptibility, Imχ_xx(ω), as a function of the film thickness, is shown in Fig. 2(a). In contrast to the 2D case³⁴ (see left panel), the spectrum for the 3D SkL has a larger number of resonance modes marked by index p. Each mode has high- and low-frequency branches that merge with decreasing thickness. The exception is the pair of modes with p = 0, which in the limit of the very thin film converges to the clockwise (CW) and counter-clockwise (CCW) precession modes of the 2D SkL.³⁴ Another remarkable feature of the pair of modes p = 0 is the presence of the crossover point at the thickness of about L_z = 0.62L_D (d = 21), at which the frequency of the CW mode becomes lower than that for the CCW mode. This effect is absent for the modes with p ≥ 1.

Modes belonging to the high-frequency branches have a lower intensity that quickly decays with the thickness. For p > 2, the high-frequency and low-frequency branches overlap and hybridize, which makes identifying a particular mode challenging. Because of that, in the following, we focus on the modes with p ≤ 2 and film thicknesses below 2L_D only.

The functional dependence of the resonance frequency on the thickness is ω_p(d) ∼ p/d, which is most prominent for low-frequency modes. Such behavior of the resonance frequencies resembles the behavior of SSW. We show below that observed phenomena indeed can be explained in terms of SSW. Furthermore, we demonstrate that the presence of the DMI plays an essential role in this phenomenon and gives rise to several interesting effects that essentially distinguish SSW in chiral magnets from SSW observed in non-chiral magnets.

The dependence of Imχ_zz(ω) on the film thickness for B_AC‖B_DC is provided in Fig. 2(b). In the limiting case of one layer, d = 1 (see left panel for the 2D SkL), there is only one resonance mode—the so-called breathing mode.³⁴ With increasing thickness, the number of resonance modes increases [Fig. 2(b)]. Similar to the case of B_AC ⊥B_DC, the resonance frequency and the intensities of those modes decrease with increasing thickness. On the other hand, the functional dependencies for ω_p(d) for the cases of in-plane and out-of-plane AC fields are different.

To illuminate the physical origin of the spectra, shown in Figs. 2(a) and 2(b), we analyze the profile of the dynamic component of the magnetization δn for different resonance modes. For both cases of in-plane and out-of-plane AC fields, we have chosen three representative thicknesses and the frequencies corresponding to the first three modes, p = 0, 1, and 2 [see red hollow circles in Figs. 2(a) and 2(b)]. On the left panel in Figs. 2(c)–2(k), we show representative snapshots of δn(x, y, z), for each of those modes. Since the δn is not a unit vector field, we use a special color scheme that reflects both the direction and the length of the vectors [see Fig. 1(c)]. When |δn| → 0, the color converges to gray, and for |δn| → δn_max, the saturation of the color reaches the maximum value. Here, δn_max is the maximal value of the |δn| over the whole domain. The details of δn(x, y, z) calculations are provided in  Appendix A.

In the case of B_AC ⊥B_DC, for low-frequency modes [see Figs. 2(c)–2(e)], the |δn| is maximum in the inter-skyrmion area where magnetization is parallel to the B_DC. In the case of high-frequency modes [Figs. 2(f)–2(h)], the maximum of the |δn| is mainly located in the core of the skyrmions, where magnetization is antiparallel to the B_DC. Since the high-frequency modes are localized in a much smaller volume, their intensities in the absorption spectrum are much lower than that of low-frequency modes.

For B_AC‖B_DC, |δn| takes its maximal values in the skyrmion shell where magnetization is mainly lying in the plane, n ⊥e_z. Nevertheless, the distribution of δn across the thickness is similar to the case of B_AC ⊥B_DC and is characterized by the presence of clearly visible minima and maxima. In the right panel of Figs. 2(c)–2(k), we provide the profiles of the dynamic magnetization averaged in a certain area Ω in the xy-plane, ⟨δn⟩≡⟨δn⟩(t, z) = ∑_(x,y)∈Ωδn(t, r), which aim to make the presence of the minima and maxima more evident. The details of the ⟨δn⟩ profiles calculation and choosing the area Ω are provided in  Appendix A. The chains of red dots correspond to the snapshots of the ⟨δn⟩(z) at a fixed moment in time. The yellow–blue surfaces of revolution are created by a set of such snapshots taken at equal time intervals (see  Appendix A).

The dynamic magnetization surfaces have nodes and antinodes. The number of nodes equals the index of the corresponding mode, p. A major difference between the low-frequency and the high-frequency modes is the opposite sense of ⟨δn_xy⟩ rotation. The vectors ⟨δn_xy⟩ rotate counterclockwise (CCW) in the case of low-frequency modes (see supplementary material Movies 1–3) and clockwise (CW) in the case of high-frequency modes (see supplementary material Movies 4–6). For breathing modes, the spins rotate CCW about the chosen precession axis (see supplementary material Movies 7–9). The sense of ⟨δn⟩ rotation is defined by the LLG equation and orientation of the projection plane normal with respect to the effective field H.

Since the position of the nodes and antinodes remains stable in time, for the case of B_AC ⊥B_DC, we identify the corresponding excitation modes as chiral SSW. We use the term “chiral” because the profile of the dynamic magnetization ⟨δn⟩(z) has a helical profile with fixed chirality. On the contrary, in the case of SSW in ferromagnet without DMI, the vectors ⟨δn⟩ always lie in the same plane.^22,35 Moreover, for the case of B_AC ⊥B_DC in both cases of the low-frequency and high-frequency modes depicted in Figs. 2(c)–2(h), the chirality of the ⟨δn⟩(z) profile across the thickness is defined by the chirality of the DMI. The period of such chiral modulation of the ⟨δn⟩(z) profile equals L_D.

For the case of B_AC‖e_z, the profile of ⟨δn(z)⟩ is also characterized by the presence of nodes and antinodes [Figs. 2(i)–2(k)], and, thus, it can be thought of as SSW localized in the shells of skyrmion tubes. The ⟨δn(z)⟩ profile, in this case, does not have well-defined chirality as in the case of the CCW and CW modes excited by the in-plane field. On the other hand, the ⟨δn(z)⟩ profile is not flat as in the case of a pure ferromagnet. The latter can be explained by the presence of the chiral surface twist^20,21 and the curvature of the skyrmion shell.

An interesting aspect of the problem is the oscillatory fading and increasing of the dynamic susceptibility seen in the low-frequency branches of the absorption spectrum for the in-plane AC field [see Fig. 2(a)]. To explain this phenomenon, we use the standard approach for analysis of the spin waves.^22,36,37 Since the low-frequency modes of the skyrmion lattice are localized in the regions where n₀‖e_z, we will consider the limiting case of the uniformly magnetized film and resonance of the spin-waves propagating in opposite directions along the z axis. Following the approach of Refs. 25 and 32, we derived the analytical equations for the dynamic magnetization of the eigenmodes,

where the term cos(πpz/d) defines the position of the nodes and antinodes, while Dz/J is an additional phase shift across the film thickness (for details, see  Appendix C). The resonance frequency of such standing spin-wave is a function of the thickness d and the mode’s index p,

Figure 3(a) shows good agreement between the results obtained with spin-dynamics simulations and the analytical expression (4) for B_DC = 1.1B_D. The only exception is the limit of very thin films where the continuum model and discreet spin lattice model diverge.²¹ 

In Figs. 3(c)–3(k), we show a few representative examples of the dynamic magnetization profiles for different modes and thicknesses, which were calculated with (3) and (4). The red lines represent the instant dynamic magnetization profiles at the fixed moment in time, and the surfaces are formed by a complete set of such snapshots over one period. It is easy to see that analytically derived dynamical magnetization profiles in Figs. 3(f)–3(h) are consistent with that obtained in the numerical experiment for the skyrmion lattice under in-plane excitation [see Figs. 2(c)–2(e) and 2(f)–2(h)]. For comparison, in Figs. 3(c)–3(e), we show the δn profiles corresponding to the case of SSW in pure ferromagnets, where, at any fixed time, all the vectors of dynamic magnetization lay in the same plane. From the projections of the δn profiles into the plane orthogonal to the z axis [Figs. 3(c)–3(e) on the right], it is seen that the total dynamic magnetization, $∫0dδndz$⁠, is nonzero only for p = 0. Because of that, in pure ferromagnets, the modes with p > 0, strictly speaking, cannot be excited by a uniform AC field.²² In chiral magnets, it is not the case because the δn profile is not flat. On the other hand, the surface of the revolution created by the instant snapshots of the dynamical magnetization profiles is identical to chiral and flat standing spin-waves. Thereby, in the most general case, one should take into account a complete 3D profile of the dynamical magnetization to identify the position of the nodes and antinodes. Otherwise, when only one component of the dynamic magnetization is taken into account, one can come to a contradicting statement that the position of the antinodes is not fixed in time, and thus, the observed spin-waves cannot be identified as standing.³² The variation of the projected on one plane dynamic magnetization of the chiral SSW is illustrated in supplementary material Movies 10.

To understand the reasons for the fading of the intensity of the chiral standing spin-waves, one should take into account the total dynamic magnetization. It should be noted that the projection of the δn into the plane orthogonal to the z axis has a curved shape [Figs. 3(f)–3(h) on the right]. Because of that, the total dynamic magnetization is not zero, and chiral standing spin-waves can be easily excited even by the uniform AC field. On the other hand, because of the additional phase shift in (3), the total dynamic magnetization is a function of the thickness. In Figs. 3(i)–3(k), we provide three representative examples of the δn profiles at particular thicknesses, where projected δn represents closed loops, and the total dynamic magnetization integrated over the thickness tends to zero. The latter explains the fading zones in the absorption spectrum, which are marked in Fig. 3(a) by (i)–(k), respectively. When the total dynamic magnetization of a particular mode tends to zero, the response of the system to the uniform AC field reduces, and we observe the fading of the intensity of that eigenmode in the absorption spectrum.

To quantify this effect, we estimated the intensity I_p of the resonance modes using the following well-known relation:^22,36,37

where δn ≡ δn(p, z, t, d) for the particular mode p and thickness d are defined by (3) and (4). In Fig. 3(b), we show the dependencies I_p for the modes with index p from 0 to 6. The position of the minima and maxima in I_p(d) dependencies is fully consistent with the fading and increase of the Imχ(ω) in the simulated absorption spectrum in Fig. 3(a). One can conclude that the effect of the fading of the intensities of particular chiral SSWs originates from the chiral profile of the dynamic magnetization. Noticeably, this effect also appears in the spectrum of the conical phase (see  Appendix D). The chiral profile of the dynamic magnetization, in this case, appears not because of the additional phase shift across the thickness but due to the chiral modulations in the static equilibrium state.

The above analytical approach can also be applied to the case of n₀ ⊥e_z, which can be compared to the modes excited in the skyrmion shell. The dispersion relation ω(k) of spin-waves propagating along the z axis, in this case, is symmetric as in the pure ferromagnet.³⁸ In this case, the solutions corresponding to SSWs are identical to those for pure ferromagnet, and the dynamic magnetization profile is flat as in Figs. 3(c)–3(e). Thereby, SSWs are not chiral in this case. For more details, see  Appendix C. On the other hand, in the skyrmion shell, the static equilibrium magnetization across the film thickness is not ideally collinear. Instead, it exhibits an additional chiral twist near the free surface of the sample.^20,21 Because of the modulation of n₀ ≡ n₀(z), similar to the case of the cone phase, the total dynamic magnetization becomes not zero, which allows exciting these modes with the uniform AC field. Moreover, the magnetization in the skyrmion shell has a small out-of-plane component. The spin-wave spectrum ω(k) in this case becomes asymmetric,³⁸ and the dynamic magnetization ⟨δn_zy⟩ gets a weak twisting across the film thickness. The twist induced by the magnetization with the positive and negative out-of-plane components has an opposite sign. This effect is compensated in the ideal case of an isolated flat domain wall. However, since the shell of the skyrmion has a cylindrical shape, the volumes with the positive and negative out-of-plane magnetization are not identical, and the volume that corresponds to the outer region dominates. The latter gives rise to a weak twist of dynamical magnetization in the case of SSWs in the skyrmion shell.

In summary, we studied the eigenmodes of the 3D skyrmion lattice in the film of the isotropic chiral magnet. To illustrate the features of the spin-wave resonance in a 3D skyrmion lattice, we calculated the absorption spectra for the cases of in-plane and out-of-plane excitations. In the case of in-plane excitation, the absorption spectrum calculated as a function of the film thickness is characterized by a set of low-frequency and high-frequency modes that appear in pairs. In this case, the maxima of the spin-wave excitations are localized in the inter-skyrmion area or in the skyrmion core, where magnetization is normal to the film plane. We identify the corresponding resonance modes as chiral standing spin waves. In contrast to standard standing spin waves in ferromagnets without DMI, the profiles of dynamic magnetization of chiral standing spin waves are characterized by helical modulation with the pitch determined by the competition between DMI and Heisenberg exchange. An important consequence of such modulation of dynamic magnetization is the periodic fading of the absorption spectra intensity with the variation of the plate thickness. A detailed analysis revealed that this effect takes place not only in the 3D skyrmion lattice but also in the case of the filed saturated state and the conical phase excited by the uniform field.

Under out-of-plane excitation, the absorption spectrum also demonstrates the appearance of standing spin waves, which are localized in the skyrmion shell where spins lie in the plane of the plate. The helical modulations of the dynamic magnetization profile are not present. Thus, standing spin waves, in this case, have more in common with standard standing spin waves in pure ferromagnets. We show that the effect of the chiral surface twist plays an essential role in this case.

See supplementary material Movies 1–9 [URL:] for more information about the rotation of dynamic magnetization vectors from Figs. 2(c) to 2(k). supplementary material Movie 10 [URL:] demonstrates the variation of dynamic magnetization in 3D space compared to its 2D projection for the case of Fig. 3(g).

The authors thank S.V. Tarasenko for fruitful discussion. The authors acknowledge financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant No. 856538, project “3D MAGiC”). F.N.R. acknowledges support from the Swedish Research Council.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The LLG equation [Eq. (2)] for underlying Hamiltonian (1) was numerically solved by the fourth-order Runge–Kutta method. The spectra were calculated with Excalibur code³⁹ and double-checked with the publicly available code Mumax.⁴⁰ The modes visualization was performed with the code Magnoom.⁴¹ 

The following parameters were used: B_DC = 0.5B_D, D/J = 0.18 (L_D = 34.91a), α = 0.01, and the time step dt = 0.05. To excite the oscillations in the given ground state, we follow the standard approach.^34,42 First, we excite the system with the time-step function⁴² with the amplitude B_AC = 0.03B_DC. After the excitation field is switched off at t = 0, we record the total magnetization of the system at each iteration. To ensure a high resolution of the spectra, we set the total number of iterations to 2¹⁹ ≈ 5 × 10⁵. Using the discrete Fourier transform of n(t), we calculate the dynamic magnetic susceptibility χ(ω) as a function of oscillation frequency ω. The absorption spectra presented in the main text show the imaginary part of dynamic susceptibility Imχ(ω) for different film thicknesses. We also tested the magnetization excitation by the sinc function^43,44 and found that, in this case, it does not qualitatively influence the Imχ(ω) dependencies. As follows from (2), the unitless frequency ω is related to frequency f in Hz as ω = f(2πμ_s)/(γJ). For instance, for J = 0.4 meV, ω = 0.01 corresponds f ≈ 1 GHz.

To identify the distribution of dynamic magnetization δn depicted on the left panel of Figs. 2(c)–2(k) and the corresponding profiles of averaged dynamic magnetization ⟨δn⟩(z), the relaxed 3D SkL was excited by an AC magnetic field B_AC = A sin ω_pt at the corresponding resonance frequency ω_p. Since we consider small oscillations, the amplitude of exciting field A = 0.01B_DC for B_AC‖e_x and A = 0.00125B_DC for B_AC‖e_z. We perform simulations in the long time range to reach the dynamic equilibrium regime. After that, we took the snapshots of the entire magnetization vector field with the time interval Δt = 2π/(ω_pN), where N = 32 is the number of the snapshots. To ensure precision in these calculations, we fix the time step dt in the LLG simulation to be multiple to Δt, meaning Δt/dt is always an integer. In addition, the snapshots of the magnetization are averaged over ten periods. Then, for each snapshot, we calculate the normalized dynamic component of magnetization at each ith site, δn_i(t) = (n_i(t) − n₀)/δn_max, where δn_max = max |n_i(t) − n₀| is the maximal dynamic magnetization over the whole sites and whole period of oscillation, i.e. |δn_i|⩽1 for any t and i. The normalized that way vector field δn_i(t) can be visualized using the color scheme explained in Fig. 1(c). In the left panel of Figs. 2(c)–2(k), we show a representative example of the δn distribution over the simulated domain at a fixed moment in time.

Because the magnetization precession axis and localization of the modes are different [Fig. 2], the area Ω for calculation of the average dynamic magnetization profiles has been chosen differently, as illustrated in Fig. 4. For the low-frequency modes, Figs. 2(c)–2(e), this area occupies the whole xy-plane in the simulation domain, Ω = L_x × L_y. Although δn, in this case, is mainly localized in the inter-skyrmion region, the contribution from other parts of the domain is negligibly small, and for simplicity, the averaging area Ω can be extended to the whole domain.

For the high-frequency modes [Figs. 2(f)–2(h)], we bound the averaging area by a disk, Ω = 2πR², where R is defined as a distance from the skyrmion center at which the dynamic magnetization reaches minima, δn(|r| = R) → min. Such a choice of Ω allows us to exclude small excitations that appeared in the inter-skyrmion region because of the coupling to more intense low-frequency modes.

In the case of the in-plane AC field, the dynamics magnetization mainly lies in the xy-plane (δn_z component is negligibly small). We assume that the normal to the xy-plane is parallel to e_z in both cases of the low-frequency and high-frequency modes. As seen from Figs. 4(a) and 4(b), in this case, the projected into xy-plane dynamics magnetization precesses CCW and CW, respectively.

The contrary to above, for modes excited by the out-of-plane AC field [Figs. 2(i)–2(k)], there is no common precession axis for the spins lying in the skyrmion shell. In this case, we choose the area Ω as a line segment connecting the skyrmion center and an arbitrary point on the circle of radius L_y/2 around it. A particular choice of the point on the circle does not play a role. For definiteness, we select that line segment such that it is parallel to the y axis, as shown in Fig. 4(c). In this case, the dynamics magnetization projected to the yz-plane exhibits CCW rotation because we select the normal to the yz-plane to be along −e_x. Otherwise, the profiles of the average dynamics magnetization would precess CW.

To estimate the role of the dipole–dipole interaction (DDI), we used the micromagnetic simulations in Mumax.⁴⁰ We used the fourth-order Runge–Kutta method implemented in Mumax to solve the micromagnetic LLG equation. The relationships between the micromagnetic constants and the constants of the model (1) have the form⁴⁵ 

where $A$ and $D$ are micromagnetic constants for exchange and DMI, respectively, M_s is the magnetization of the material, and a is a lattice constant. The dimensionless time t in our spin model simulations and dimensional time τ in micromagnetic simulations are related as

For definiteness, we used the following micromagnetic parameters:⁴⁶ $A=4.75$ pJ/m, $D=0.853m$J/m², and M_s = 384 kA/m and performed the simulations on the mesh with 78 × 45 × 50 cuboids in the x, y, and z directions, respectively. Following the procedure described in  Appendix A, we calculate the absorption spectra for in-plane excitation (B_AC ⊥B_DC‖e_z, B_DC = 0.5B_D) for both cases with and without DDI. The results are shown in Fig. 5(a). It is shown that DDI results in a slight shift of some eigenmodes frequencies, but the number of the resonant peaks in a given frequency region is the same. For comparison, in Fig. 5(b), we provide the spectrum calculated for the spin model at the same parameters.

We identified the in-plain dynamic magnetization profiles of eigenmodes for the case with DDI, following the method described in the main text and  Appendix A. We show these profiles in the insets of Fig. 5(a) for the first and second low-frequency eigenmodes. Similar to the cases without DDI in Fig. 2, the dynamic magnetization surfaces also have nodes and antinodes and thus correspond to SSWs.

A slight shift of the resonance frequencies occurs because the dynamic magnetization profiles, ⟨δn⟩(z) in the cases with and without DDI, are slightly different. In particular, in the case of DDI, the dynamic magnetization near the free surfaces twists a bit faster than in the case without DDI.

For comparison, in Fig. 5(b), we show the spectra calculated in the spin lattice model without DDI. The discrepancy between the spectra calculated without DDI in the micromagnetic model and the spin lattice model, Figs. 5(a) and 5(b), respectively, can be explained by a relatively high ratio of D/J used in our simulations. When D/J gets smaller, the spectra calculated with two models converge. Alternatively, one can show the consistency of the simulations performed in different software by correcting the finite difference scheme in calculating effective fields. The solid blue line in Fig. 5(b) corresponds to the spectra obtained in Mumax with such a correction. One can see the perfect agreement. For more details, see the Mumax script provided in the supplementary materials.

Here, we recall the results of Refs. 25 and 32 and show that the analytical solutions for the profile of dynamic magnetization are fully consistent with the results of our spin-dynamic simulations. We consider the case when the magnetization n₀ in static equilibrium is parallel to e_z in a whole sample. The spin wave propagates along the z axis. Since the amplitudes of dynamic magnetization δn are small, one can take into account only orthogonal components δn_x and δn_y of dynamic magnetization and neglect the δn_z component. For unconfined crystals, this leads to the following dispersion relation:²⁵ 

where $k̄=ka$ and ω are the absolute values of the dimensionless wave vector and the frequency of a spin-wave, respectively. H in (C1) is the absolute value of the effective (internal) magnetic field, which exists in static equilibrium H_i‖n_0i at any site. The sign “±” in (C1) defines the spin-wave propagation direction. In the limiting case of D = 0, the dispersion relation (C1) converges to the well-known solution for a bulk ferromagnet,

where $k̄FM$ is the wave vector in pure ferromagnet. Using circular variables, δn = δn_x + iδn_y, the superposition of two spin waves propagating in the opposite directions along the z axis in the case of a chiral magnet can be written as

where $k̄r$ and $k̄l$ are the wave vectors for spin waves propagating parallel and anti-parallel z axis. In the case of a pure ferromagnet, the dynamic magnetization corresponding to the superposition of two waves takes the form

Using (C1)–(C4), the vectors $k̄r$ and $k̄l$ can be written as $k̄r=k̄′+k̄FM$ and $k̄l=k̄′−k̄FM$⁠, respectively.

As follows from (C1), for small D/J, the wave vector shift, up to the first leading term, is $k̄′=D/J$⁠. Let us assume that the sample has a shape of a film in the xy-plane with the boundaries located at z = 0 and z = d. The boundary conditions are²⁵ 

By substituting (C3) in (C5), we obtain the set of solutions with nonzero amplitudes, $k̄FM=πp/d$⁠, where p = 0, 1, 2, …. Then, the in-plane components of the dynamic magnetization for the SSW take the form of (3) in the main text.

Now let us consider the case when the magnetization lies in the film plane, n₀ ⊥e_z. For definiteness, we will assume that n₀‖e_x [see Fig. 4(c)]. For spin waves propagating along the z axis, the dispersion relation ω(k) is symmetric.³⁸ In the case of an unconfined crystal of a chiral magnet, the dispersion relation has the form (C2). For the plate, using (C4) and boundary conditions (C5), we get the solution for the SSW with the “flat” profile of ⟨δn_zy⟩,

Thereby, the SSW is not chiral in this case. The same is true for SSW localized in the shell of a skyrmion.

Let us consider the case when n₀ in static equilibrium represents a cone phase with the wave vector Q = 2π/L_D parallel to the z axis, $n0=sinΘ⁡cosΦ,sinΘ⁡sinΦ,cosΘ$⁠, where Θ and Φ = Qz are polar and azimuthal angles, respectively. For definiteness, in the following, we assume that the static magnetic field B_DC‖e_z and B_DC = 0.85B_D are present, and the conical phase is the ground state of the system in a wide range of film thickness.²¹ When spin wave propagates along e_z, the dispersion relations for bulk chiral magnet have the form²⁵ 

where $Q̄=Qa$⁠. Using (C5), one can show that the eigenfrequency of the SSW is²⁵ 

which remains valid only for B_DC ⩽ B_D.

Figure 6(a) shows the absorption spectrum Imχ_xx for cone phase obtained in numerical simulations as a function of the film thickness. The eigenfrequencies of the SSWs according to Eq. (D2) show good agreement with spin-dynamics simulations. Similar to the case of the spectrum for the field saturated state in Fig. 3, the only exception is the limit of very thin films where the continuum model and discreet spin lattice model diverge.²¹ According to Eq. (D2), the eigenfrequency ω₀ = 0 at any thickness, while in numerical simulations at small thicknesses, it has a nonzero value.

Similar to the case of the saturated state, the absorption spectra for the cone phase are also characterized by periodical fading of the mode intensities. To make the periodical decay of Imχ_xx more visible, we show the analytical dependencies (dashed red lines) only for p ≥ 4. In contrast to the saturated state (C1), the ω(k) relation for the cone phase (D1) is symmetric. The SSW, in this case, is formed by two opposite propagating spin waves with identical k-vector. To explain the presence of the fading zones following the period of cone modulations, L_D in Fig. 6(a), one should take into account the spiraling of the magnetization in a static equilibrium state [Fig. 6(b)]. To take this into account, we apply the rotation matrix associated with the cone phase to the two counterpropagating spin-waves in a ferromagnetic state (C4). Then, by substituting the obtained δn(z, t) into (5), we calculate the dependencies of the eigenmode intensities as functions of the film thickness, which are depicted in Fig. 6(c). The position of the minima and maxima in Figs. 6(a) and 6(c) show excellent agreement. Interestingly, in the coordinate system related to the static equilibrium state of the cone phase illustrated in Fig. 6(b), the δn profiles of SSW are flat [Fig. 6(d)], as in the case of eigenmodes of the pure ferromagnet. The profiles of δn in Fig. 6(d) are obtained in numerical simulations. The coordinate transformation depicted in Fig. 6(b) effectively unwinds the cone phase into the state with n₀‖e_z. Thus, the nonzero intensities of the SSW exited by the uniform AC field in our numerical experiment are solely explained by the twist of magnetization n₀ in the static equilibrium state. The fading of the intensities following the period of cone modulations, L_D, appears at the thickness, which satisfies the criteria I_p(d) → 0 in (5).