Controlling the three dimensional propagation of spin waves in continuous ferromagnetic films with an increasing out of plane undulation

The role of three-dimensionality in a ferromagnetic medium in ruling the propagation properties of spin-waves (SW) has been one of the main focuses of the research activity in recent years. In this context, we investigate the evolution of the SW dispersion (frequency vs wave vector) induced by a progressive vertical undulation of a ferromagnetic film. The geometric undulation is taken along a single direction and is periodic with constant period, while the amplitude (differential maximum height with respect to the film thickness) is gradually increased from 0 to 60 nm. We study the characteristic modification of the internal effective field and link it to the resulting SW dispersions and spatial profile. These systems display at once features both of a planar film and a discretized medium, and the dispersion curves change not only when SWs propagate along the undulation direction, but also perpendicular to it. We discuss the geometric and magnetic conditions for having either the invariance of the SW group velocity with respect to even major changes in the undulation, or a large group velocity for some edge modes. We address a potential dual-band activity, namely the simultaneous propagation of two independent SW-signals, with separated frequency bands and disjoint oscillation regions.


Controlling the three dimensional propagation of spin waves in continuous ferromagnetic films with an increasing out of plane undulation Federico Montoncello 1* & Gianluca Gubbiotti 2
The role of three-dimensionality in a ferromagnetic medium in ruling the propagation properties of spin-waves (SW) has been one of the main focuses of the research activity in recent years. In this context, we investigate the evolution of the SW dispersion (frequency vs wave vector) induced by a progressive vertical undulation of a ferromagnetic film. The geometric undulation is taken along a single direction and is periodic with constant period, while the amplitude (differential maximum height with respect to the film thickness) is gradually increased from 0 to 60 nm. We study the characteristic modification of the internal effective field and link it to the resulting SW dispersions and spatial profile. These systems display at once features both of a planar film and a discretized medium, and the dispersion curves change not only when SWs propagate along the undulation direction, but also perpendicular to it. We discuss the geometric and magnetic conditions for having either the invariance of the SW group velocity with respect to even major changes in the undulation, or a large group velocity for some edge modes. We address a potential dual-band activity, namely the simultaneous propagation of two independent SW-signals, with separated frequency bands and disjoint oscillation regions.
The propagation of waves in periodically structured materials is of fundamental interest in modern physics and technology, and when spin waves (SWs) in ferromagnetic media are concerned, the corresponding research field is Magnonics 1-3 . In Magnonics, SWs play the role of information carriers, which are exploited to process and store information and operate in ultra-low-power conditions. When propagating along waveguides consisting of artificial magnetic materials with properties periodically varied in space (magnonic crystals), SWs undergo Bragg diffraction and acquire new properties, absent in bulk materials, due to the formation of frequency bands and band gaps. Since these properties depend on geometry and applied magnetic field, they can be easily designed and reconfigured. For these reasons, Magnonics is considered particularly promising for a possible beyond-CMOS technology 4 . In the last decades, the conceptual problem was just studying SW propagation, which is inherently a two-dimensional problem, since a signal travelling along a waveguide can be embedded in one-(1D) or twodimensional (2D) space, depending on the curvature [5][6][7][8][9] . Hence, the research exploration has mainly involved structures displaying 1D/2D magnetization textures, either periodic [10][11][12][13] or aperiodic [14][15][16][17] , while structures with in-plane periodicity, but inhomogeneous along the perpendicular direction, have been poorly explored so far. However, exploring the third dimension is crucial indeed, because a 3D distribution of the magnetization offers a new degree of freedom, which in general allows, from one side, to fit more functionality into a smaller space, and hence to considerably increase the density of elements in magnonic devices 18,19 , and from the other, further possibilities for the SW dynamics and transport (e.g., vertical magnon transport, nonreciprocal coupling [20][21][22] ). More recently, the opportunity of designing 3D ferromagnetic systems where ferromagnetic layers are linked by a large number of vertical connections has been proposed in the perspective of the so-called neuromorphic networks [23][24][25] .

Results
System modeling. As discussed in section "Methods", we use a micromagnetic simulator to compute the equilibrium magnetic configuration, and the dynamical matrix method (DMM) to get the SW frequency and profile at any given wavevector value and direction: both calculation tools are based on the discretization of the sample into a tridimensional mesh of prisms. Given such a framework, we conceived the undulation as a sequence of infinitely long elemental stripes with thickness L = 30 nm but varying width w and displacement d along the z coordinate, depending on the geometric, out-of-plane, undulation amplitude ( Fig. 1). Each stripe is made of three 10-nm thick layers, and the displacement along z between adjacent stripes is h = 10 nm, so that L = 3 h and the generic displacement of any stripe, with respect to the ground (z = 0), is d = Mh, with integer M. We choose to fix the undulation period to a = 720 nm. Hence, if N is the number of stripes forming the primitive cell, the overall extent of the system along the z direction is d MAX N = N 2 h + L , while the undulation amplitude is defined as A = d MAX N − L as shown in Table 1. In practice, the planar film has zero undulation amplitude, while the 12-stripe system attains the largest amplitude, i.e., A = 60 nm.
The magnetization distribution has peculiar features, depending on the direction of the applied field (with constant magnitude B 0 ≡ µ 0 H 0 = 0.1T ) with respect to the undulation direction, which is kept along x-axis. This is because the geometric arrangement acts as a constraint for the magnetic moments, and determines a specific internal effective field µ 0 H eff , which is the net sum of the external field (in our case, 0.1 T), and the internal dipolar and exchange fields. µ 0 H eff is a reference quantity, because it rules both the static distribution of the magnetization and its dynamics.
If the magnetization M is along y, i.e., perpendicular to the undulation direction, then no geometric disturbance of the magnetic moment alignment occurs, and hence the distribution is uniform ( ∇ · M = 0 , i.e., the magnetization flux lines are parallel) and so is the internal effective field ( ∇ · H eff = 0 , see Fig. 2a). The only change from the planar arrangement Fig. 1a to the undulated ones (b-e) is the height of the adjacent stripes along x-direction, which implies a small increase of the exchange field due to the lack of any adjacent moments at the edges of each stripe, but a considerable decrease of the total dipolar field since, on the average, the adjacent stripes (and the magnetic moments inside) are farther when misaligned than when aligned to form the planar film. Altogether, this implies a progressive increase of µ 0 H eff with respect to the undulation amplitude, resulting in an overall increase of the SW frequency.
Conversely, if the magnetization M is along x, i.e., parallel to the undulation direction, then the distribution is no more uniform, and magnetic charges form, with density ρ M = −∇ · µ 0 M ; in turn, also the internal effective field varies along x : in Fig. 2b-e, in each primitive cell H eff (x) displays plateaux in the stripe bulk and peaks at the stripe edges, repeated for the number of stripes forming the primitive cell. In Fig. 2f we illustrate the behavior of the magnetization inside each stripe, and plot the magnetization components M z and M x along the x direction. Both plots exhibit peaks at the same heights, showing unambiguously that the actual lattice constant, for a unitary primitive cell, is different from the geometric lattice constant, and it corresponds to the width of the single stripe w This observation supports the finding, discussed in the following, that a gap opening (mode repulsion) would be possible only if the effective primitive cell was w, and not a, and hence zone boundary was π w and not π a . In an ideally continuous film w → 0 , and no gap is ever possible, despite the nonuniformity.
Dispersion relations and mode profiles. Through the DMM, we calculated the dispersion curves, in the five different systems (Fig. 1a-e), and the different magnetization configurations (i.e., M y and M x ), for SWs propagating either perpendicular or parallel to the magnetization direction, (i.e., k ⊥ M or k M ). where δm k (r) is the cell function, of which we'll be showing the real z-component profiles in the following, r = x, y is the in-plane position vector and k = k x , k y is the wavevector. In our calculations, we investigated only on the two main limit directions, i.e., k = 0, k y and k = (k x , 0).
(1) δm k (r) = δm k (r)e ik·r Figure 1. Geometric structures used in the study: (a-e), in perspective, the six primitive periodic cells (darker blue) and a repetition of each (lighter blue), as an illustration of the infinite repetitions along both in-plane coordinates; the periodic cell is 720 nm wide, and the geometric undulation is attained through a number of elemental stripes of width w. (a) is the planar film, with w = 720 nm equal to the (virtual) periodic cell width; (b) 4-stripe system, with w = 180 nm; (c) 6-stripe system, with w = 120 nm; (d) 8-stripe system, with w = 90 nm; (e) 12-stripe system, with w = 60 nm. (f) The reference system and the geometric displacement h = L/3 of adjacent stripes, with constant thickness L = 30 nm and varying width w. The undulation amplitude A is specified below each panel. www.nature.com/scientificreports/ As a first point, intending to provide a reference, we plot in Fig. 3 the analytical dispersion relation of the surface Damon-Eshbach (DE) mode (k ⊥ M) and the volume backward (BA) mode (k M) in a wide range of wavevector values. For the latter (BA), the range is enough to display both the negative slope at low wavevector (due to dipolar effects) and the parabolic rise above 2.5 × 10 7 rad/m (due to exchange effects). The equations for the dipolar-exchange SWs are 38,39 : and where ex is the exchange length, which defines the range below which is not likely to have a nonuniform spin distribution, namely: is uniform independently of the undulation amplitude (from the planar film to the 12-stripe system); despite the many geometric discontinuities, due to the different heights of the adjacent stripes, H y eff (x) is uniform (computational confidence within 0.05 kA/m) independently of the applied field magnitude. In (b-e), when M x , instead, H x eff (x) is highly nonuniform: (b) 4-stripe system; (c) 6-stripe system; (d) 8-stripe system; (e) 12-stripe system. In (f) we plot both magnetization components M z (x)(orange curve, left axis) and M x (x)(blue curve, right axis) of the top layer in the (representative) 6-stripe system, when M x ; a sketch of the magnetization along the primitive cell is juxtaposed, illustrating the behavior of its out-of-plane component across the vertical sections of the 6 stripes. Since our purpose is to compare the dispersions to the cases with periodic undulations, it is useful to fold the dispersions to the reduced Brillouin zone [0, π a ] (in our case a = 720 nm). Besides, in calculations, we do use periodic boundary conditions even when simulating the flat, infinite film, so we can plot in Fig. 4 both the analytical curves and the results of the calculations with the DMM: the agreement is very good, despite the completely different approaches, and the unavoidable computational errors.
As a general comment about Fig. 3, with a particular emphasis on the BA curve, we notice that the exchange contribution becomes appreciable above 2.5 × 10 7 rad/m, where the curve slope apparently changes sign. Such a wavevector value corresponds to a wavelength of 250 nm, hence we cannot expect exchange effects for larger wavelengths.
Now we comment the results of the DMM applied to the systems of Fig. 1, concerning the effects of the undulation on the dispersions, when the magnetization is either perpendicular or parallel to the undulation direction.
Configuration M y. With reference to the first Brillouin zone, we observe that the increasing undulation amplitude produces an almost linear rigid shift in frequency of the dispersions, upward for DE mode, downward for the BA mode, with almost no alteration of the curvature, particularly for the BA mode (Fig. 5). Both  www.nature.com/scientificreports/ modes have an ideal uniform cell function (Eq. 1), even though DMM output shows an artefact undulation, characteristic of the numerical calculation ( Fig. 6), which is not a real effect. For either modes, at k = 0 the average incremental frequency shift, on increasing undulation amplitude A , is almost constant and equal to about �ν/A = 0.014GHz/nm (which means 0.14GHz for every 10-nm-increase, following Table 1). This evidence implies that the group velocity v g is almost not sensitive (i.e., is "robust") 40 to the exact amplitude of the undulation transverse to the propagation direction, and might be considered a topological invariant of the system (see Table 2). This means that imprecision in the geometry does not impact the dynamics (propagation speed and phase relationship), which is a rather attractive feature, when SWs are operating as information carriers in 3D networks [41][42][43] . The absence of any frequency gap at zone boundary is not surprising, since the magnetization is uniformly distributed (full saturation along y-direction), and the geometric discontinuity along x-direction is not creating any magnetic poles.
Configuration M x. In this case, when undulation has non-zero amplitude, both magnetization M and H eff (x) are no more uniform, since the primitive cell is made of elemental stripes each magnetized along the width (hard axis). As it appears in Fig. 2c-e, H eff (x) displays pronounced discontinuities and sharp peaks, which separate regions where SW modes can localize; as a consequence, we found the occurrence of two main families of modes, the bulk modes (BM), localized in the center of each stripe, and the edge modes (EM) localized at the boundaries between adjacent stripes (Fig. 7). Now, for M x , the general behavior of the dispersions, independently of the SW wavevector direction and mode character BA (Fig. 8) or DE ( Fig. 9) is: from the planar film to the 4-stripe system, when the undulation amplitude A varies from zero to 20 nm, the dispersion experiences a significant frequency drop due to the enhanced dipolar fields, having moved from a uniform to a non-uniform magnetization. Then, for increasing A (and hence, the number of stripes N , Table 1) the frequency increases due to the decreasing dipolar energy (discussed above), since adjacent stripes shift far apart. As far as the BA mode is concerned (Fig. 8), this effect is accompanied by a progressive decrease of the mode bandwidth (hence, the group velocity). In fact, the increasing undulation amplitude apparently decreases the coupling among the magnetic moments, since the same amount of magnetic moments is distributed over an increasingly larger arc (the magnetic moment linear density decreases in moving from a straight line to the arc), and this impacts the dynamic fields too, which are responsible for the mode bandwidth. An additional effect arises here: no bulk uniform mode can exist below a given critical stripe width. Indeed, since each stripe is magnetized along the narrow width, the strong dipolar (demagnetizing) fields prevent any coherent spin oscillation involving the bulk of each stripe. For this reason, no bulk solution is found for the 12-stripe system. In this case, only edge solutions exist (see either Fig. 7e or Fig. 10d, which are the same), as found in the literature for similar or analogue systems, with a strongly varying H eff (x) [44][45][46][47][48] . Frequency increases about 0.14 GHz every A = 10 nm. The curves corresponding to the 2nd zone are less accurate, because of the finite discretization of the sample. Despite propagating along the direction where magnetization is uniform (y), the BA mode is nonetheless sensitive to the undulation amplitude changes. www.nature.com/scientificreports/ Comments about the occurrence of a frequency gap. It is worth to recall how the opening of a gap at zone boundary comes as a result of two conceptually independent effects. The first one is the linear superposition of two propagating plane waves (with = 2a, i.e., k = ± π/a), incident and reflected by the medium discontinuities, which interfere because of periodicity and combine into two stationary waves with a quarter of wavelength path difference. The second one relies on the medium physical properties, which must be differently experienced by the two mentioned stationary waves: otherwise, the two waves are indistinguishable with respect to the medium and degenerate in frequency: in order this could occur, the periodic cell must be such that its center has different properties from the edges. While it is clear that a uniform magnetization (configuration M y ), despite distributed in a curved geometry, cannot offer such symmetry breaking, it is less trivial when M x and k x (BA mode), and must be based on the specific magnetization profile. In a recent paper, this effect was already seen but not explained 49 . Details will be given in the "Discussion" section below.
As far as the DE mode is concerned (Fig. 9), we found a similar frequency drop, and a similar SW bandwidth reduction, moving from the planar film to the 4-stripe system, i.e., varying A from zero to 20 nm. In addition, an apparent, remarkable effect is that the dispersion for the 8-stripe system ( A = 40 nm) seems intermediate between the 4 and the 6-stripe systems, depending on the wavevector. This effect is here just spotted, but would require wider studies to be fully understood. Nevertheless, we observe how the frequency variations among the three systems (4, 6, and 8-stripe systems) are in any case very small in this configuration, around or below  Table 2. Values of the average group velocity (m/s) of the spin wave propagating along either coordinates, for the two magnetic configurations. As explained in the article body, the largest value is got by an edge mode, surprisingly. In configuration M x no bulk cell function is found. The negative sign in the values indicates the negative slope of the dispersion (BA character). www.nature.com/scientificreports/ 0.5 GHz, so that the frequency variations due to different undulation amplitudes, at a specific k, is just a matter of computational error, or, in possible measurements, experimental resolution. Finally, we discuss the dispersion curves for the uniform edge mode (profile shown in Fig. 10) either when k y (Fig. 11, showing no periodicity) or k x (Fig. 12, showing periodicity since oscillation is along the undulation direction). The edge mode is a Bloch spin wave with the cell function showing significant amplitude only In (e), the stripe width is so small (60 nm) that no bulk-like solution is possible even in the single elemental stripe, and the most uniform mode shows nonzero amplitude only at the edges of each elemental stripe: this profile Is qualitatively different from the others. Note that each map is a 3 × 3 repetition of the primitive cell, as an indicative extract of the ideally infinite system. www.nature.com/scientificreports/ at the elemental stripe edges. Edge modes can exist only because of a non-uniform H eff (x) , hence only when the film is not a straight plane. We see how the increase in the undulation amplitude causes a frequency increase for all 4, 6, and 8-stripe system, but not for the 12-stripe system, which shows a different behavior in both cases. This diversity occurs because the width of the elemental stripe is too small 45 (i.e., w = 60 nm) and the corresponding dipolar fields are very stronger. For this reason, in both Figures the curve relevant to the 12-stripe system departs from the others. Remarkably, the 12-stripe system shows a significant bandwidth (independently of the propagation direction), which is a consequence of the large dynamic stray fields caused by the dynamic magnetization.

Discussion
As a consequence of the specific representation we chose for the film undulation, we have the following results. Firstly, the spin wave periodicity is observed only along the undulation direction, expressly on the DE mode when the film is magnetized perpendicular to the undulation direction (y-axis), and the BA mode when the  www.nature.com/scientificreports/ film is magnetized parallel to the undulation direction (x-axis). Moreover, the variation in frequency is experienced by spin waves independently of the propagation direction or the magnetization configurations, because the geometric changes influence the whole system energetics. In other words, the undulation changes vary the whole system total energy, which in turn can make the system either more stable (stiffer), determining a general frequency increase, or less stable (softer), determining a general frequency decrease 50,51 . Interestingly, due to the specific symmetry of the primitive cell that we chose, we do not have any frequency gap across the first and second  www.nature.com/scientificreports/ Brillouin zone for any propagating spin wave. The two wavefunctions, forming at zone boundary ( = 2a ) as the consequence of diffraction, happen to have both node and maxima at the edges across two adjacent stripes, and independently of its direction (x or y) the magnetization shows the same flux line orientations in each stripe along the primitive cell (see Fig. 13), so that a shift of λ/4 would not change how a solution experiences the energetics of the system: hence, they happen to be degenerate in frequency in principle. Since the magnetization distribution along x direction is the same independently of the position of the elemental stripe in the primitive cell, in general, Bragg diffraction occurs when and in our case N (number of elemental stripes in a primitive cell) is always even, so we would not see any mode repulsion at any zone boundary, i.e., at any N 2 π a . In this class of problems, i.e., undulated films, the absence of any mode repulsion at zone boundary is independent of the elemental stripe size w though which it is constructed, and, in the limit w → 0 (smooth continuous film), mode repulsion must be still absent and hence the forbidden gap cannot be found at any k, despite periodicity, as a very general feature. Periodicity is in fact a necessary but not sufficient condition for a forbidden gap formation, the other necessary condition is either discontinuity in any physical parameter, say f , across the primitive cell, typically its left limit differing from the right limit f a − � = f a + (as in 1-D lattice of stripes of different materials), or non-uniformity in case of continuity of f , e.g. f a 2 = f (0) , i.e., the value of f at the center of the primitive cell must be different from that at the cell boundaries (as either in nanomagnet or antidot lattices). In this way, since the energetics is ruled by f , the two diffracted waves with an optical path length difference δ = 4 ≡ a 2 would both have = 2a , which corresponds to zone boundary ( k = π a ), but would experience differently the physical parameter f , getting a different frequency. In our problems, f is the magnetization function M(x) , and, importantly, due to the chosen geometry, not only it is always continuous (i.e., M a − = M a + ), but furthermore, the rule f a 2 = f (0) , i.e., in our case, M a 2 = M(0) , holds even when the undulated film is magnetized along the undulation direction, with a non-uniform magnetization. Besides, as already remarked, both diffracted waves experience the same magnetization flux lines, even if at different points in the primitive cell. Hence, we can't expect any forbidden gap in any case in principle.
Edge modes. Going a little further, we discuss a few other remarkable features. Despite localized in narrow regions, the edge modes in configuration M x show appreciable frequency bandwidth, in particular it increases on increasing undulation amplitude. Usually localization decreases the dynamic coupling fields, but in our geometry, the number of edges increases with the undulation amplitude (i.e., the number of stripes) and hence, the more the edges between adjacent stripes, the larger the dynamic coupling fields and hence the mode bandwidth. Their appreciable average group velocity allows to look to these modes as information carriers too, as done for their bulk counterparts. We computed the average group velocities of all the SW modes in either configurations M y and M x , collected in Table 2. While, as expected, the fastest bulk wave is a DE mode in the planar film, surprisingly the absolute fastest spin wave is an edge mode for the maximum undulation amplitude k = π w = Nπ a Figure 13. Sketch of the correspondence between magnetization undulation and spin oscillation wavelength at zone boundary (k = π/a) for the system made of 6-stripes. Black dashed lines set the reference system, the blue and red lines represent the two diffracted waves at zone boundary, having a phase difference of π/2. Apparently, even when M is not uniform, the two solutions experience point by point exactly the same magnetization flux lines (illustrated by the blue arrows), resulting indistinguishable and hence in principle degenerate in frequency. This happens any time that λ/4 is an integer of the stripe width w. www.nature.com/scientificreports/ in configuration M x . In ordinary conditions, spin wave localization fits with low frequencies as a consequence of high dipolar fields 52 : but in our geometry, due to the large number of stripes forming the primitive cell (in the 12-stripe system), and magnetized along the hard axis, there is a summation effect which enhances the dynamic stray fields. In a sense, limiting to our results, we might say in a catch phrase that undulation accelerates a special class of spin waves, the backward edge mode.
The reader might observe that any real sample, in configuration M x , would not have the strong discontinuities we found that are associated to the serrated shape of our stripe-based representation: the real magnetization should in fact vary smoothly. We reply noticing that nevertheless a divergence would display peaks at each change of curvature of the magnetic undulation (specifically, the component M x ), signaling the presence of magnetic poles and effective walls for spin waves, dividing the primitive cell into bulk and edge regions (Fig. 14). Hence, even in a smoothly varying geometry, undulated films can host edge modes as well, similarly to what found in our calculations. The difference is that only two edge regions per primitive cell can appear, while in our discretized case we had two edge regions per stripe. Besides, a device actually built of stripes might be of interest in the context of a dual-band signal delivery: as apparent by our calculations, edge mode frequencies are below 8 GHz, while bulk BA/DE modes are beyond 9 GHz, so there is a gap of more than 1 GHz, and the excitation areas are complementary (bulk/edges) and the related SW intensity is considerable: hence the perspective of simultaneous delivery of two signals is conceivable in appropriate devices inspired by our composite structure, as suggested in the introductory paper ref. 53 .
Finally, to broaden the generality of our results, we suggest that the continuous magnetization undulation could be attained not only by a geometric deformation of the planar film but also through any appropriate nanomagnet lattice layers on a planar film, or even through a multiferroic multilayer that, via magnetoelastic coupling, can provide a periodic effective undulation amplitude of the magnetic anisotropy, and hence magnetization, which can be tuned on demand by any external electric fields (electrically driven magnonics).
In conclusion, we investigated how a geometric undulation in a continuous film modifies the phase profile and dispersion relation for the SW not only propagating along the undulation direction, but also perpendicular to it. We studied the cases when the film is magnetized parallel or perpendicular to the undulation direction, i.e., when the magnetization is uniform or non-uniform, and provide rules for the occurrence (or absence) of a forbidden gap at zone boundary. We discussed the coexistence of periodic waves together with non-periodic waves, and, on the other side, of bulk waves together with edge waves when the magnetization is not uniform. We addressed the conditions at which the SW group velocity becomes topologically invariant against undulation alterations. We found how special edge modes, forming in such a 3D structure, show a large frequency bandwidth comparable to or even larger than the bulk wave ones, hence allowing an additional channel for simultaneous information delivery.

Methods
We used the object oriented micromagnetic framework (OOMMF) simulator 54 , with implemented periodic boundary conditions, to compute the distribution of the equilibrium magnetization for each system, as well as the internal effective field profiles. The primitive cell, a square with side 720 nm and a number of 10-nm-thick layers depending on the overall height d MAX N (Table 1), was discretized through micromagnetic cells with dimensions 5 × 5x10 nm 3 and represented by a series of adjacent, 30-nm thick stripes following the geometry shown in Fig. 1. The magnetic parameters were typical for permalloy material: exchange stiffness constant A = 10.0pJ/m , saturation magnetization M S = 800kA/m and the gyromagnetic ratio γ = 185radGHz/T . The magnetization maps were the input of the DMM software 55,56,59 , which computes the frequency and phase profiles of the dynamic magnetization both at any given applied field and wavevector magnitude and direction. In our case the applied field was µ 0 H 0 = 0.1T , while the wavevector magnitude k was taken at discrete points, from 0 to 0.5 with steps of 0.1 in units of the Brillouin zone width, and then converted into absolute values rad/m. www.nature.com/scientificreports/ Outline of the DMM. As discussed at the beginning of the Results Section, we divide the primitive cell into a mesh of micromagnetic cells (prisms), inside which magnetization is uniform: this assumption is justified by the small size of these micromagnetic cells, comparable to the exchange length. Hence, our problem consists both of many micromagnetic cells interacting with each other inside the primitive cell, and also many primitive cells (at different lattice points) interacting with each other. Each micromagnetic cell, across the mesh, has a fixed magnetic moment, which precesses around its local equilibrium direction, identified by the equilibrium azimuthal and polar angles ϕ, θ (found through OOMMF). The equilibrium direction in a specific cell is determined by the effective field generated by all the other magnetic moments (through dipolar and exchange interactions) and the applied external field (if present). In the precession motion of a single spin S , the generalized coordinates q n are the small variations both of azimuthal (q 1 = δϕ) and polar q 2 = δθ angles around the corresponding rotation axis, and the conjugate momenta (p 1 and p 2 ) are the projection of the related spin variation δS to the corresponding rotation axis ( e rotation ) 57,58 : with n = 1, 2 . The precessing spin S is linked to the magnetization M = M s m through the following rule: where µ is the magnetic moment of a single elemental prism (micromagnetic cell), V is the volume of the prism, γ is the gyromagnetic ratio, M s is the saturation magnetization of the material, and the unit vector m is expressed in spherical coordinates. With this relation, the generalized coordinates and conjugate momenta can be written as (Fig. 15): Then, the Hamilton equations, relevant for a single micromagnetic cell (i.e., a single magnetic moment) read: where H is the Hamiltonian of the system, i.e., the energy function including Zeeman, exchange and dipolar interactions (for the micromagnetic representation of these interactions see refs. 55,59 ). After substitution we get:  www.nature.com/scientificreports/ Now, we have to consider that the Hamiltonian includes interactions between a given micromagnetic cell and all the (for dipolar) or some (for exchange) others. Hence it is necessary to index the cells: we address a given reference micromagnetic cell with index k, while all the others, interacting with that, will be indexed with j. The total number of cells in the micromagnetic mesh is N (maximum value of the indices, and, hence, system size).
Furthermore, we exploit the fact that, for harmonic precession, we also have: where δϕ 0 k and δθ 0 k are the oscillation amplitudes and ω/2π the frequency. Then, we introduce the energy density U = H/V , and, since a harmonic behavior involves small oscillations, we expand it in Taylor series up to second order. Considering that we deal with equilibrium systems (found by OOMMF) we retain only the quadratic term, and then, with just a simple algebraic manipulation, we get a linear system, describing the precession of the magnetization in the kth-cell, which is interacting with all the other jth-cells of the system: which can be used to get the SW frequencies ω k /2π , and the eigenvectors: which represent the precession angles of the magnetization in each micromagnetic cell and can be transformed into the full SW mode spatial profiles considering that k goes from 1 to N. Note that the extension of the ferromagnetic body along the three cartesian coordinates x, y, z is just a matter of indexing, so that position r of any micromagnetic cell inside the three-dimensional mesh can be seen as r x, y, z ≡ r k . Clearly, the system size N determines also the number N of possible normal modes (v k ).
When dealing with periodic systems, we have a set of lattice points at position R (any linear combination of the primitive vectors). Note that our systems involve all the three dimensions, but periodicity is only along x, y . At each lattice point R we have exactly the same magnetization distribution that we have in the primitive cell at R = 0 (and the same index k to address all the micromagnetic cells at position r k for any R = 0 ). The linear system above [Eq. (4)], in this new case, corresponds to the dynamics within a single primitive periodic cell only, and hence must be extended to include the interactions of this primitive cell (usually taken, as a reference, at R = 0 ) with all the other periodic cells at other lattice points R ′ (introducing a summation over R ′ ). Since the magnetization has the same distribution M x, y, z in any periodic cell, exploiting the Bloch theorem and δϕ k = δϕ 0 k e iωt δθ k = δθ 0 k e iωt www.nature.com/scientificreports/ introducing the wavevector K, this can be easily shown to give an extra factor (i.e., e iK·R ′ ) in Eq. (4), and the dynamic equation for periodic systems is transformed into: Here, index k runs for the micromagnetic cells inside the reference primitive cell ( R = 0 ), while index j runs for the micromagnetic cells in any primitive cell at lattice point R ′ (including R ′ = 0 ). Even in this case, after an appropriate change of rows/columns, the linear system can be cast into an eigenvalue/eigenvector problem, which has exactly the same system size N of the previous one Eq. (4), thanks to the translational symmetry and the Bloch theorem. This time, we have an additional parameter, which is the wavevector K at which the calculation is performed, which must be provided as an input, and in general can be any arbitrary point in the reciprocal space. Hence, the solutions (eigenvectors) must be intended as Bloch waves following Eq. (1), which (with the present notation) reads: being δm K (r) = sinθ k δϕ k δθ k the cell function of the dynamic magnetization, and r the generic position vector, which, in the discretized periodic sample, depends on the micromagnetic cell index k (going from 1 to N in the primitive cell), and possibly includes a lattice vector R ′ to span across the full periodic system, i.e., r x, y, z = r k + R ′ . The computation and diagonalization of the dynamical matrix, out of any equilibrium magnetization distribution, is performed by an appropriate software, developed in parallel with the theoretical model at the University of Ferrara.