This work studies the resonance excitations of the three-dimensional skyrmions lattice in the finite thickness plate of an isotropic chiral magnet using spin dynamics simulations. We found that the absorption spectra and resonance modes differ from those predicted by the two-dimensional model and the model of the unconfined bulk crystal. The features observed on the spectra can be explained by the formation of chiral standing spin waves, which, contrary to conventional standing spin waves, are characterized by the helical profile of dynamic magnetization of fixed chirality that is defined by the Dzyaloshinskii–Moriya interaction. In this case, the dynamic susceptibility becomes a function of the plate thickness, which gives rise to an interesting effect that manifests itself in periodical fading of the intensity of corresponding modes and makes excitation of these modes impossible at specific thicknesses.

## I. INTRODUCTION

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−x}Co_{x}Si,^{13,14} and other compounds from that class,^{15} 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_{2}OSeO_{3}.^{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.^{21}

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,^{22} 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^{32} or at least cannot be considered as conventional SSW.^{33} 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.

## II. MODEL

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} = *D***r**_{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*^{2}/(*μ*_{s}*J*) 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\u22483$. 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} = 4*L*_{D}. For definiteness, we set *J* = 1 and *D* = 0.18, which results in *L*_{D} ≈ 35*a* and $BD=0.0324\mu s\u22121$. 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*_{0}/*μ*_{s}, with *γ* being the gyromagnetic ratio and *E*_{0} 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=\u2212E0\u22121\u2202E/\u2202ni$ is a dimensionless effective field on the *i*th 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**_{0}, and dynamic, *δ***n**(*t*), components: **n**(*t*) = **n**_{0} + *δ***n**(*t*) [see Fig. 1(b)]. For visualization of the fields **n**_{0} 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.

## III. RESULTS

### A. The resonant spectra

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^{34} (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.^{34} 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.62*L*_{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 2*L*_{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.^{34} 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.

### B. Visualization of resonance modes

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.

## IV. DISCUSSIONS

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**_{0}‖**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.1*B*_{D}. The only exception is the limit of very thin films where the continuum model and discreet spin lattice model diverge.^{21}

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, $\u222b0d\delta 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.^{22} 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.^{32} 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**_{0} ⊥**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.^{38} 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**_{0} ≡ **n**_{0}(*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,^{38} 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.

## V. CONCLUSIONS

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.

## SUPPLEMENTARY MATERIAL

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).

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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

### APPENDIX A: METHODS

#### 1. Parameters of the model and spectra calculation

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^{39} and double-checked with the publicly available code Mumax.^{40} The modes visualization was performed with the code Magnoom.^{41}

The following parameters were used: *B*_{DC} = 0.5*B*_{D}, *D*/*J* = 0.18 (*L*_{D} = 34.91*a*), *α* = 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^{42} with the amplitude *B*_{AC} = 0.03*B*_{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^{19} ≈ 5 × 10^{5}. 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.

#### 2. Calculation of the distribution and profile of the dynamic magnetization

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 *ω*_{p}*t* at the corresponding resonance frequency *ω*_{p}. Since we consider small oscillations, the amplitude of exciting field *A* = 0.01*B*_{DC} for **B**_{AC}‖**e**_{x} and *A* = 0.00125*B*_{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*π*/(*ω*_{p}*N*), 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 *i*th site, *δ***n**_{i}(*t*) = (**n**_{i}(*t*) − **n**_{0})/*δn*_{max}, where *δn*_{max} = max |**n**_{i}(*t*) − **n**_{0}| 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*^{2}, 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.

### APPENDIX B: EFFECT OF DIPOLE–DIPOLE INTERACTIONS ON RESONANCE EIGENMODES

To estimate the role of the dipole–dipole interaction (DDI), we used the micromagnetic simulations in Mumax.^{40} 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^{45}

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:^{46} $A=4.75$ pJ/m, $D=0.853m$J/m^{2}, 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.5*B*_{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.

### APPENDIX C: RESONANCE EIGENMODES OF SATURATED STATE

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**_{0} 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:^{25}

where $k\u0304=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\u0304FM$ 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\u0304r$ and $k\u0304l$ 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\u0304r$ and $k\u0304l$ can be written as $k\u0304r=k\u0304\u2032+k\u0304FM$ and $k\u0304l=k\u0304\u2032\u2212k\u0304FM$, respectively.

As follows from (C1), for small *D*/*J*, the wave vector shift, up to the first leading term, is $k\u0304\u2032=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^{25}

By substituting (C3) in (C5), we obtain the set of solutions with nonzero amplitudes, $k\u0304FM=\pi 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**_{0} ⊥**e**_{z}. For definiteness, we will assume that **n**_{0}‖**e**_{x} [see Fig. 4(c)]. For spin waves propagating along the *z* axis, the dispersion relation *ω*(*k*) is symmetric.^{38} 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.

### APPENDIX D: ABSORPTION SPECTRA AND RESONANCE EIGENMODES OF CONICAL SPIN SPIRAL STATE

Let us consider the case when **n**_{0} in static equilibrium represents a cone phase with the wave vector *Q* = 2*π*/*L*_{D} parallel to the *z* axis, $n0=sin\Theta \u2061cos\Phi ,sin\Theta \u2061sin\Phi ,cos\Theta $, 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.85*B*_{D} are present, and the conical phase is the ground state of the system in a wide range of film thickness.^{21} When spin wave propagates along **e**_{z}, the dispersion relations for bulk chiral magnet have the form^{25}

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.^{21} According to Eq. (D2), the eigenfrequency *ω*_{0} = 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**_{0}‖**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**_{0} 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).

## REFERENCES

_{c},”

_{2}OSeO

_{3},”

_{2}OSeO

_{3},”

_{2}OSeO

_{3}studied by a novel electron spin resonance technique,”

_{2}OSeO

_{3},”