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 materials1,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 Fe1−xCoxSi,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 Cu2OSeO3.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 resonance23,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 magnets32 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 ni = Mi/μs is a unit vector of the magnetic moment at site i, μs is the value of the magnetic moment of the lattice site, means the summation of overall nearest–neighbor pairs, J is the exchange coupling constant and Dij is the Dzyaloshinskii–Moriya vector defined as Dij = Drij, D is the DMI scalar constant, and rij is the unit vector between sites i and j. The external magnetic field, B(t) = BDC + BAC(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 BD = D2/(μsJ) and equilibrium period of helical modulations in the continuum limit, LD = 2πaJ/D, respectively. While the values of BD and LD 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 BD and LD 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, . The thickness of the film Lz varies by changing the number of atomic layers, d, from 1 to 140, which for a chosen coupling parameter corresponds to Lz = 0 (in the single atomic layer limit—2D case) and Lz = 4LD. For definiteness, we set J = 1 and D = 0.18, which results in LD ≈ 35a and . The static magnetic field is always perpendicular to the film plane, BDC‖ez. In the following, we study both cases of in-plane and out-of-plane alternate magnetic fields: BAC ⊥BDC and BAC‖BDC.
(a) The simulated domain with the size fits the equilibrium period of the 3D skyrmion lattice with periodic boundary conditions in the xy plane. The color code represents the magnetization field n0(x, y, z). The static magnetic field is perpendicular to the film plane, BDC‖ez. (b) In the excited state, the spins oscillate about the direction of the effective field and can be decomposed into static component n0 and dynamic component δn at resonance. (c) Color code used for a unit vector field of magnetization n and its dynamic component of magnetization δn.
(a) The simulated domain with the size fits the equilibrium period of the 3D skyrmion lattice with periodic boundary conditions in the xy plane. The color code represents the magnetization field n0(x, y, z). The static magnetic field is perpendicular to the film plane, BDC‖ez. (b) In the excited state, the spins oscillate about the direction of the effective field and can be decomposed into static component n0 and dynamic component δn at resonance. (c) Color code used for a unit vector field of magnetization n and its dynamic component of magnetization δn.
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 γE0/μs, with γ being the gyromagnetic ratio and E0 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 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, n0, and dynamic, δn(t), components: n(t) = n0 + δn(t) [see Fig. 1(b)]. For visualization of the fields n0 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.
(a) Absorption spectra Imχxx(ω) of the 3D skyrmion lattice for different thicknesses d = Lz/a of chiral magnet film under an in-plane exciting field BAC‖ex. (b) The thickness dependent spectra of Imχzz(ω) in the case of an out-of-plane exciting magnetic field BAC‖ez. The white digits in (a) and (b) correspond to the indices of modes. The left panels in (a) and (b) correspond to the limiting case of the 2D system. The red–white dots in (a) and (b) indicate the values of the thickness and frequency used for the visualization of dynamic magnetization distributions δn(x, y, z) and spin-wave profiles depicted below in (c)–(k). (c)–(e) Snapshots of δn(x, y, z) in the simulated domain at a fixed time and profiles of the average value of ⟨δnxy⟩ along the film thickness for low-frequency modes with indices p = 0, 1, 2 correspondingly for BAC‖ex. (f)–(h) Snapshots of δn(x, y, z) and ⟨δnxy⟩ for high-frequency modes with indices p = 0, 1, 2 for BAC‖ex; white dashed line indicates the average region for ⟨δnxy⟩; (i)–(k) Snapshots of δn(x, y, z) and ⟨δnzy⟩ of the modes with indices p = 0, 1, 2 for BAC‖ez.
(a) Absorption spectra Imχxx(ω) of the 3D skyrmion lattice for different thicknesses d = Lz/a of chiral magnet film under an in-plane exciting field BAC‖ex. (b) The thickness dependent spectra of Imχzz(ω) in the case of an out-of-plane exciting magnetic field BAC‖ez. The white digits in (a) and (b) correspond to the indices of modes. The left panels in (a) and (b) correspond to the limiting case of the 2D system. The red–white dots in (a) and (b) indicate the values of the thickness and frequency used for the visualization of dynamic magnetization distributions δn(x, y, z) and spin-wave profiles depicted below in (c)–(k). (c)–(e) Snapshots of δn(x, y, z) in the simulated domain at a fixed time and profiles of the average value of ⟨δnxy⟩ along the film thickness for low-frequency modes with indices p = 0, 1, 2 correspondingly for BAC‖ex. (f)–(h) Snapshots of δn(x, y, z) and ⟨δnxy⟩ for high-frequency modes with indices p = 0, 1, 2 for BAC‖ex; white dashed line indicates the average region for ⟨δnxy⟩; (i)–(k) Snapshots of δn(x, y, z) and ⟨δnzy⟩ of the modes with indices p = 0, 1, 2 for BAC‖ez.
III. RESULTS
A. The resonant spectra
First, we consider the case of BAC ⊥BDC. The dynamic susceptibility, Imχxx(ω), as a function of the film thickness, is shown in Fig. 2(a). In contrast to the 2D case34 (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 Lz = 0.62LD (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 2LD 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 BAC‖BDC 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 BAC ⊥BDC, 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| → δnmax, the saturation of the color reaches the maximum value. Here, δnmax 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 BAC ⊥BDC, 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 BDC. 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 BDC. 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 BAC‖BDC, |δn| takes its maximal values in the skyrmion shell where magnetization is mainly lying in the plane, n ⊥ez. Nevertheless, the distribution of δn across the thickness is similar to the case of BAC ⊥BDC 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 ⟨δnxy⟩ rotation. The vectors ⟨δnxy⟩ 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 BAC ⊥BDC, 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 BAC ⊥BDC 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 LD.
For the case of BAC‖ez, 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 twist20,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 n0‖ez, 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 BDC = 1.1BD. The only exception is the limit of very thin films where the continuum model and discreet spin lattice model diverge.21
(a) Absorption spectra Imχ(ω) for a film of chiral magnet in the field polarized phase (BDC‖ez, BDC = 1.1BD) under in-plane excitation (BAC ⊥ez). The red dashed lines are the resonance frequencies (4). The symbols indicate the corresponding figures with visualizations of dynamic magnetization profiles across the film thickness. (b) The normalized intensity of calculated with (5) as a function of the film thickness for the modes with p = 0–6. (c)–(e) Dynamical magnetization profiles (C4) for the case of a ferromagnetic film without DMI. (f)–(h) Profiles of eigenmodes excited by in-plane uniform field BAC in the magnetic film with DMI in field polarized phase. (i)–(k) Profiles of eigenmodes unexcited by in-plane uniform field BAC in the magnetic film with DMI in field polarized phase. The solid black circle and hollow red circle in (c)–(k) correspond to the ends of the vectors ⟨δnxy⟩ at z = 0 and z = d, respectively.
(a) Absorption spectra Imχ(ω) for a film of chiral magnet in the field polarized phase (BDC‖ez, BDC = 1.1BD) under in-plane excitation (BAC ⊥ez). The red dashed lines are the resonance frequencies (4). The symbols indicate the corresponding figures with visualizations of dynamic magnetization profiles across the film thickness. (b) The normalized intensity of calculated with (5) as a function of the film thickness for the modes with p = 0–6. (c)–(e) Dynamical magnetization profiles (C4) for the case of a ferromagnetic film without DMI. (f)–(h) Profiles of eigenmodes excited by in-plane uniform field BAC in the magnetic film with DMI in field polarized phase. (i)–(k) Profiles of eigenmodes unexcited by in-plane uniform field BAC in the magnetic film with DMI in field polarized phase. The solid black circle and hollow red circle in (c)–(k) correspond to the ends of the vectors ⟨δnxy⟩ at z = 0 and z = d, respectively.
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, , 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 Ip 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 Ip for the modes with index p from 0 to 6. The position of the minima and maxima in Ip(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 n0 ⊥ez, 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 n0 ≡ n0(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 ⟨δnzy⟩ 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 code39 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: BDC = 0.5BD, D/J = 0.18 (LD = 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 function42 with the amplitude BAC = 0.03BDC. 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 219 ≈ 5 × 105. 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 function43,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 BAC = A sin ωpt at the corresponding resonance frequency ωp. Since we consider small oscillations, the amplitude of exciting field A = 0.01BDC for BAC‖ex and A = 0.00125BDC for BAC‖ez. 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, δni(t) = (ni(t) − n0)/δnmax, where δnmax = max |ni(t) − n0| is the maximal dynamic magnetization over the whole sites and whole period of oscillation, i.e. |δni|⩽1 for any t and i. The normalized that way vector field δni(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, Ω = Lx × Ly. 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.
The area Ω in the xy-plane for the averaged dynamic magnetization, ⟨δn⟩(t, z) = ∑(x,y)∈Ωδn(t, r), for different modes. (a) and (b) correspond to the case of BAC ⊥BDC for the low-frequency and high-frequency modes, respectively. (c) corresponds to the case of BAC‖BDC. The top row of images shows the simulated domain, the arbitrarily chosen xy-plane, and the area Ω bounded by the blue line. The bottom row of the images shows the snapshots of the δn distribution provided for the reference [see also Figs. 2(c), 2(f), and 2(i)]. In (a), Ω corresponds to the whole xy-plane, in (b), Ω is bounded by the circular domain of the radius corresponding to the lowest value |δn| around the skyrmion core, and in (c), Ω corresponds to the line segment drawn from the center of the simulated domain and parallel to the y axis.
The area Ω in the xy-plane for the averaged dynamic magnetization, ⟨δn⟩(t, z) = ∑(x,y)∈Ωδn(t, r), for different modes. (a) and (b) correspond to the case of BAC ⊥BDC for the low-frequency and high-frequency modes, respectively. (c) corresponds to the case of BAC‖BDC. The top row of images shows the simulated domain, the arbitrarily chosen xy-plane, and the area Ω bounded by the blue line. The bottom row of the images shows the snapshots of the δn distribution provided for the reference [see also Figs. 2(c), 2(f), and 2(i)]. In (a), Ω corresponds to the whole xy-plane, in (b), Ω is bounded by the circular domain of the radius corresponding to the lowest value |δn| around the skyrmion core, and in (c), Ω corresponds to the line segment drawn from the center of the simulated domain and parallel to the y axis.
For the high-frequency modes [Figs. 2(f)–2(h)], we bound the averaging area by a disk, Ω = 2πR2, 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 (δnz component is negligibly small). We assume that the normal to the xy-plane is parallel to ez 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 Ly/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 −ex. 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 form45
where and are micromagnetic constants for exchange and DMI, respectively, Ms 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 pJ/m, J/m2, and Ms = 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 (BAC ⊥BDC‖ez, BDC = 0.5BD) 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.
(a) The absorption spectra Imχxx(ω) of SkL in the plate of isotropic chiral magnetic with the thickness Lz = 1.42LD calculated under in-plane excitation (BAC ⊥BDC‖ez) in Mumax3 with and without dipole–dipole interactions. The abbreviations “LF” and “HF” denote the low-frequency and high-frequency modes with the corresponding index according to Fig. 2(a). The bottom axis corresponds to the dimensionless frequency ω [see (2) and Appendix A], while the frequency in the top axis is given in GHz and corresponds to FeGe parameters.46 The insets show the profiles of the dynamic magnetization across the plate thickness. (b) The comparison of the absorption spectra calculated in the spin-lattice model and Mumax with the correction to the finite difference scheme implemented through a built-in Mumax function that allows user-defined effective fields.
(a) The absorption spectra Imχxx(ω) of SkL in the plate of isotropic chiral magnetic with the thickness Lz = 1.42LD calculated under in-plane excitation (BAC ⊥BDC‖ez) in Mumax3 with and without dipole–dipole interactions. The abbreviations “LF” and “HF” denote the low-frequency and high-frequency modes with the corresponding index according to Fig. 2(a). The bottom axis corresponds to the dimensionless frequency ω [see (2) and Appendix A], while the frequency in the top axis is given in GHz and corresponds to FeGe parameters.46 The insets show the profiles of the dynamic magnetization across the plate thickness. (b) The comparison of the absorption spectra calculated in the spin-lattice model and Mumax with the correction to the finite difference scheme implemented through a built-in Mumax function that allows user-defined effective fields.
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 n0 in static equilibrium is parallel to ez 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 δnx and δny of dynamic magnetization and neglect the δnz component. For unconfined crystals, this leads to the following dispersion relation:25
where 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 Hi‖n0i 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 is the wave vector in pure ferromagnet. Using circular variables, δn = δnx + iδny, 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 and 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
As follows from (C1), for small D/J, the wave vector shift, up to the first leading term, is . 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 are25
By substituting (C3) in (C5), we obtain the set of solutions with nonzero amplitudes, , 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, n0 ⊥ez. For definiteness, we will assume that n0‖ex [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 ⟨δnzy⟩,
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 n0 in static equilibrium represents a cone phase with the wave vector Q = 2π/LD parallel to the z axis, , where Θ and Φ = Qz are polar and azimuthal angles, respectively. For definiteness, in the following, we assume that the static magnetic field BDC‖ez and BDC = 0.85BD 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 ez, the dispersion relations for bulk chiral magnet have the form25
which remains valid only for BDC ⩽ BD.
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.
(a) Absorption spectra Imχxx(ω) of chiral magnetic film in the conical spin spiral state for different thicknesses d = Lz/a (dark–blue–white color map) and the analytical spectrum (D2) (red dashed lines) for BDC = 0.85BD; white digits correspond to indices of modes. (b) The equilibrium magnetic ordering n0 in the conical spin spiral state along film thickness and the local coordinate system x′y′z′ associated with n0i in each ith layer (z′‖n0i). (c) The values of (5) as a function of the film thickness for modes with indices p = 1, 2, 3, 4, 5 in the conical spin spiral state. (d) Snapshots of dynamic magnetization ⟨δnx′y′⟩ in the local system x′y′z′ at fixed time moment for modes with indices p = 1, 2, 3, 4, 5.
(a) Absorption spectra Imχxx(ω) of chiral magnetic film in the conical spin spiral state for different thicknesses d = Lz/a (dark–blue–white color map) and the analytical spectrum (D2) (red dashed lines) for BDC = 0.85BD; white digits correspond to indices of modes. (b) The equilibrium magnetic ordering n0 in the conical spin spiral state along film thickness and the local coordinate system x′y′z′ associated with n0i in each ith layer (z′‖n0i). (c) The values of (5) as a function of the film thickness for modes with indices p = 1, 2, 3, 4, 5 in the conical spin spiral state. (d) Snapshots of dynamic magnetization ⟨δnx′y′⟩ in the local system x′y′z′ at fixed time moment for modes with indices p = 1, 2, 3, 4, 5.
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, LD 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 n0‖ez. 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 n0 in the static equilibrium state. The fading of the intensities following the period of cone modulations, LD, appears at the thickness, which satisfies the criteria Ip(d) → 0 in (5).