Finite-element dynamic-matrix approach for spin-wave dispersions in magnonic waveguides with arbitrary cross section

,


I. INTRODUCTION
Over the last decades, a number of powerful analytical and numerical tools have been developed to describe the linear propagation characteristics of spin waves, the fundamental small-amplitude excitations in magnetically ordered substances.For example, micromagnetic approaches, which treat the magnetization of a solid as a continuous vector field, have been used to derive approximate analytic expressions for the dispersion relation and spatial mode profiles of spin waves with wavelengths in the nano-and micrometer range.These models allow to describe the spin-wave propagation in simple geometries, such as bulk systems, infinitely extended monoor bilayers [1][2][3] or waveguides with rectangular cross section 4,5 .Apart from flat geometries, in the emerging field of curvilinear magnetism, which includes studying curvature-induced effects on spin-wave propagation, the dispersion in cylindrical waveguides, magnetic nanotubes with thin mantle 6,7 , curved nanowires 8 and narrow ribbons 9 , have been described.
Although the developed analytic approaches are versatile and applicable to many different cases, they tend to give only approximate dispersions and mode profiles.Moreover they are often not available for more complex systems.Therefore, numerical approaches are needed to complement and extend analytical or experimental studies.Dynamic micromagnetic simulations, which rely on a rigorous time integration of the equation of motion of the magnetization on a discrete mesh, have been established as one of the standard tools to study spin-wave propagation numerically [10][11][12] .In such dynamic simulations, spin waves are typically obtained by exciting the magnetic system with a microwave field pulse, letting the sys-* ) Electronic mail: l.koerber@hzdr.detem evolve according to the appropriate equation of motion, later on extracting the spatial mode profiles and computing the wavelengths and frequencies by means of Fourier analysis.This, however, requires an approximate a priori knowledge of the mode profiles and does not always allow to separate degenerate modes.
Alternatively, a numeric micromagnetic method known as dynamic-matrix approach can be used to obtain spin-wave mode profiles and frequencies directly by numerically solving a linearized version of the equation of motion of the magnetization using an appropriate eigensolver 13,14 .One advantage of such a method is that it directly yields frequencies and mode profiles without the need for additional post processing.Moreover, degenerate modes as well as modes with non-trivial spatial profile can be resolved.The dynamic-matrix approach has already been implemented successfully to study standing waves in confined magnetic elements, using both finiteelement and finite-difference methods of discretization [15][16][17][18] .Recently, in an excellent work, Henry et al. 19 succeeded to extend this dynamic-matrix approach for propagating spin waves in systems with a translational invariant magnetic equilibrium, using plane-wave demagnetization tensors.They employed a finite-difference method to efficiently obtain the spin-wave dispersion by numerically solving the linearized equation of motion only in the cross-section of the magnetic medium perpendicular to the propagation direction, modelling infinite magnetic slabs or thin flat waveguides using rectangular elements.This approach was already successfully applied e.g. to study spin waves in multilayer systems 20 and those propagating along Bloch domain walls 21 .Next to the aforementioned benefits, this propagating-wave dynamic-matrix approach allows for almost arbitrary wave-number precision in the propagation direction, something which can become extremely costly when modelling a full three-dimensional magnetic specimen.
Although the finite-difference-based propagating-wave dynamic-matrix approach is very efficient and allows to study e.g.complex multilayer systems, it is not suitable for waveguides with arbitrary (bounded) cross section, for example tubular systems, or, in general, systems with surface curvature.The aim of this paper is to provide a complimentary approach to the one presented in Ref. 19 which allows to study exactly such systems.For this purpose, we employ a different discretization type, namely the finite-element method (FEM), by modeling the cross section of an arbitrary waveguide using triangular elements.As a key difference, instead of using plane-wave demagnetization tensors, we obtain the dipolar potential of the individual spin-wave modes using the Fredkin-Koehler method 22 , which is a hybrid finiteelement/boundary-element method to solve the Poisson equation on a finite volume.Here, for the first time, we present an extension of this method for the screened Poisson equation, which governs the lateral dipolar field of propagating waves.
To this end, in Sec.II, we review the theoretical basis of the eigenvalue problem to calculate spin-wave dispersions within the context of micromagnetism.Following, in Sec.III, the different magnetic interactions considered in the present work are introduced for the case of propagating waves.The numerical implementation, including the plane-wave Fredkin-Koehler method, is described in Sec.IV.Finally, we apply our approach to different waveguide cross sections in Sec.V. We believe that this work is of particular interest for the emerging field of curvilinear magnonics, but also allows to calculate spin-wave dispersions in various standard magnonic problems.

A. Basics of micromagnetism
Within the theory of micromagnetism 23,24 , the magnetization M (r, t) of a magnetic body under isothermal conditions (far below the Curie temperature of the material at hand) is treated as a continuous vector field, subject to the constraint |M | = M s , with M s being the saturation magnetization of the material.At a constant temperature and exposed to some external field H ext (r, t), a stable magnetic equilibrium within a given specimen is reached once its Gibb's free energy is at a local minimum.In terms of the normalized magnetization field m = M /M s , the Gibb's free energy G(m) of a magnetic specimen is given by with µ 0 being the vacuum permeability, V the volume of the magnetic body, h ext = H ext /M s the unitless external field and h int the unitless internal magnetic field produced solely by the magnetization itself 25 .We have chosen the unitless representation of the magnetization and field in anticipation of the numerical implementation.Moreover, we assume to have a homogeneous saturation M s in the material.Extending this method to inhomogeneous materials does not require any difficult steps.
For most common magnetic interactions, the internal field can be expressed as being linear in the magnetization, using the self-adjoint operator N which describes the magnetic self interactions 26 .This includes the symmetric exchange interaction, the dipolar interaction, and uniaxial crystal anisotropy, The details of each interaction will be described in Sec.III.A prominent exception to the linear representation Eq. ( 2) is for example the cubic magnetic anisotropy, which can, however, be included approximately using a linearization.Next to cubic anisotropy, other interactions such as the Dzyaloshinskii-Moriya interaction can also be included in the same way.However, in the present work, we limit ourselves to the ones mentioned above, with the main focus being on the dipolar interaction.
Once the magnetization field m is not in a minimum of the free energy G, its temporal evolution is given by the Landau-Lifshitz-Gilbert equation of motion with ω M = γµ 0 M s being the characteristic magnetic frequency and γ the modulus of the gyromagnetic ratio in the material at hand.The equation of motion must be equipped with the boundary condition n • ∇m = 0 which is automatically satisfied by most choices of basis functions in a FEM approach.The first term in Eq. ( 4) describes the precessional motion around the so-called effective field h eff , comprised of the internal and the external fields.This (unitless) effective field is obtained from the Gibb's free energy as the variational derivative, 27 The second term T G , in torque equation Eq. ( 4), represents viscous (Gilbert) damping and relaxes the magnetization to equilibrium after a certain time 28 .When searching for the spin-wave dispersion in a certain magnetic body, we are interested only in the small-amplitude conservative dynamics and can therefore neglect this term.The linear damping rates of individual spin-wave modes can also be recovered from their spatial profiles 29 .To map the spin-wave spectrum, Eq. ( 4) can be solved using a numerical time integration, as it is done in many micromagnetic codes [10][11][12] .However, as discussed above, it is desirable to obtain stationary small-amplitude solutions (the spin-wave modes) directly from an eigenvalue problem.Following e.g.Refs.14,16,30,31, the nonlinear torque equation Eq. ( 4) can be linearized in the usual way, by separating the time-dependent magnetization into a static and a dynamic part according to with m 0 (r) being the normalized direction of the equilibrium magnetization and δm(r, t) being the dimensionless smallamplitude variation from the equilibrium.Since the norm of the total magnetization vector is constrained by |m| = 1, the equilibrium direction and the dynamical component must be orthogonal, m 0 ⊥ δm.Inserting the separation Eq. ( 6) into the torque equation Eq. ( 4) and neglecting all nonlinear terms, and finally expanding the variable magnetization into linear waves according to leads to the linearized Landau-Lifshitz equation which is a standard eigenvalue problem that yields the spinwave mode frequencies ω ν as well as their spatial profiles m ν .
It is clear from the demand that δm is a real vector, that if ω ν is an eigenvalue of Eq. ( 8) with corresponding eigenvector m ν , then so is −ω ν with eigenvector m * ν .Here, the asterisk denotes complex conjugation.
The Hamiltonian operator Ω is a tensor operator given by here h 0 = m 0 •h eff (m 0 ) is the projection of the unitless static effective field (including any static external field) onto the equilibrium direction m 0 and Î is the identity operator.The operator Ω is self-adjoint and, as long as m 0 is a local minimum of the Gibb's free energy G, it is also positive definite for vectors m ν perpendicular to m 0 .As a consequence, one can show that all eigenfrequencies ω ν are real valued and, moreover, eigenvectors for different eigenvalues satisfy orthogonality relations 14,16 .Let us note, that we write the linearized equation Eq. ( 8) in a similar form as e.g.Refs.16,30,31 with the only difference that our operator Ω is made dimensionless.

B. Reduction to waveguide cross section
In order to solve numerically the eigenvalue problem defined in Eq. ( 8), typically, the whole magnetic volume needs to be discretized.Depending on the implementation of the dipolar interaction one even needs to model a large "air box" outside of the specimen.The problem can be significantly simplified when studying the spin-wave spectrum in an extended waveguide with translationally invariant equilibrium along the length of the waveguide.In other words, the equilibrium m 0 (ρ) does not change along the propagation direction of the spin waves which we choose to be e z in the lab system, with ρ being a spatial vector in the xy plane.Let us stress the point, that the equilibrium magnetization can be inhomogeneous within the waveguide cross section, i.e. depend on the coordinates ρ.We consider a long waveguide with constant cross-section area A and length L → ∞ (see Fig. 1(a)).To avoid inflating notation, we denote a set as well its Hausdorff measure (e.g. the 2D measure area for A) with the same symbol.In every instance, the difference will be clear from the context.When letting the length of the waveguide L go to infinity, one can transform the linearized equation into a planewave problem and solve it only in a single cross section of the waveguide.Of course, for this case, the Gibb's free energy G of the whole magnetic element diverges.Thus, the equilibrium configuration m 0 (ρ) can be obtained by minimizing the energy density over length g(m 0 ) = G(m 0 )/L which remains finite, as L → ∞.
The spin-wave mode profiles about such a translationally invariant equilibrium can be expressed as with k being the wave number in z-direction, η νk the lateral mode profiles which are complex vector fields in the xy plane.
Here, ν is the lateral mode index.For simplicity, we will drop the tilde again and, from now on, use ν only as the lateral mode index.Using this definition, Eq. ( 8) becomes with the plane-wave Hamiltonian operator which is the original operator transformed into the waveguide cross section.Clearly, this operator is still self-adjoint.Moreover, if Ω is positive definite it follows that Ωk is too, i.e. for some arbitrary lateral profile v = e −ikz w which satisfies v ⊥ m 0 (with w being the corresponding full volumetric profile) with respect to the L 2 (S) scalar product Therefore, the spectral properties of the initial eigenvalue problems are recovered.By solving Eq. ( 10) numerically with the wave vector k as a parameter, one obtains the dispersion ω ν (k) ≡ ω νk of different mode branches ν within the waveguide at hand.Due to the complex eigenmode ansatz proportional to exp[i(kz − ω νk t)], again, all lateral mode profiles come in complex pairs.If η νk is the corresponding eigenvector to the physical eigenvalue ω νk , then η * νk is the eigenvector to the conjugate eigenvalue −ω ν,−k .In other words, the graphs (k, ω νk ) of the individual dispersion branches ν are, in general, point symmetric with respect to (0, 0).

C. Rotated eigenvalue problem
Recall that we want to solve the plane-wave eigenvalue problem Eq. ( 11) for eigenvectors orthogonal to the equilibrium magnetization, m 0 ⊥ η νk .The implementation of this constraint can be achieved either by introducing an operator which projects vectors into the plane locally transverse to the equilibrium direction m 0 or by rotating the eigenvalue problem into the local reference frame of m 0 and reducing the dimension of the eigenvectors 14 .The rotation is performed using a suitable unitary transformation R into a local orthonormal system {e 1 , e 2 , e 3 } where m 0 = e 3 everywhere (see Fig. 1(b)).Such a basis is for example given by By left-multiplying Eq. ( 11) with this operator R and using its unitarity, i.e.R † R = Î, one obtains the rotated eigenvalue problem with the rotated eigenvectors ηνk = Rη νk , the dynamicmatrix operator and the dagger denoting Hermitian conjugation.Due to the rotation of the mode profiles into the local reference frame of m 0 , one realizes that, under the constraint that the mode profiles must be locally perpendicular to the static magnetization, their third component must vanish everywhere, i.e. ην • e 3 ≡ 0. This allows to reduce the dimensionality of the profiles from three-dimensional to two-dimensional vectors.
This transformation can be performed using the matrix operator which is unitary for vectors orthogonal to m 0 .In the local reference frame, the cross product m 0 × ... takes the simple form of a multiplication with the matrix (written in {e 1 , e 2 } basis)

III. COMPONENTS OF THE PLANE-WAVE HAMILTONIAN OPERATOR
After the plane-wave eigenvalue problem has been formulated, we can introduce the components of the plane-wave operator by obtaining the plane-wave versions Nk of the magnetic tensor for the different interactions considered.For each interaction, the contribution to the equilibrium-field projection h 0 can of course be obtained by applying the operator Nk to the equilibrium configuration m 0 at k = 0.For example, if only an external field and exchange interaction would be considered, one would have

A. Uniaxial crystal anisotropy
The uniaxial magnetocrystalline anisotropy arising from spin-orbit coupling of the spin system with the crystal lattice can be expressed as with K u being the first-order uniaxial-anisotropy constant, e u being the normalized direction of anisotropy and ⊗ denoting the tensor product.Since this interaction does not involve any spatial derivatives the resulting Hamiltonian operator Ω(uni) k does not change its form when transformed to the waveguide cross section, i.e.N(uni) = N(uni) k is not dependent on the wave vector.

B. Symmetric exchange interaction
The exchange interaction within the continuum limit is written as with λ ex = 2A ex /µ 0 M 2 s being the exchange length, A ex the exchange stiffness constant of the material and ∇ 2 the Laplacian acting on three-dimensional vector fields in threedimensional space.Projecting the operator into the waveguide cross section yields Here, ∇ 2 ρ denotes the part of the Laplacian operator which operates on a three-dimensional vector fields in the xy plane.

C. Dipolar interaction
For an arbitrary magnetic specimen, the dipolar interaction can be expressed using the so-called magnetostatic potential φ(r) as The magnetostatic potential is obtained by solving the Poisson equation with appropriate continuity and jump conditions at the boundary ∂V (see for example Eqs.(4-5) in Ref. 22).For the potential of an individual spin wave with mode profile according to Eq. ( 10), we employ the ansatz with the complex lateral potential ψ νk (ρ).Inserting this ansatz into the initial Poisson equation, we obtain Equation ( 28) is a so-called screened Poisson equation (or Yukawa equation) which, analytically, is equivalent to the inhomogeneous Helmholtz equation by making the substitution k → ik.We discuss the numerical solution of this equation in more detail in Sec.IV.Once the lateral potential ψ νk of a given mode is calculated, the action of the plane-wave dipolar operator on a lateral mode profile is given by

IV. NUMERICAL IMPLEMENTATION AND THE FREDKIN-KOEHLER METHOD FOR PLANE WAVES
For a numerical solution of the eigenvalue problem in Eq. ( 16), the equilibrium magnetization m 0 (ρ), mode profiles η νk (ρ), the potentials ψ νk (ρ) as well as all involved operators in the dynamic matrix D can be discretized within the cross section A of the waveguide by projecting them onto triangular elements with linear shape functions.Here, we use a mass-lumping technique to obtain the discretized differential operators.After discretization, all operators including ∇ or the transformation R (except the dipolar operator) take the form of sparse matrices and their actions are simply matrix-vector multiplications.Therefore, the numerical implementation can be kept close to the mathematical formulation.Once the dynamic matrix is constructed, it can be diagonalized using standard iterative eigensolvers such as the Arnoldi-Lánczos method 32,33 .
The evaluation of the dipolar interaction requires special care.The lateral potential ψ νk in the dipolar operator N(dip) could be readily obtained by the convolution of the right-hand side of the screened Poisson equation Eq. ( 28) with the appropriate Green's function using the modified Bessel function of second kind and zeroth order K 0 .On a discrete mesh with n nodes, this approach requires the computation of O(n 2 ) matrix elements of the Green's function 22 .Moreover, due to the logarithmic singularity at |ρ − ρ | = 0 and the rapid decay for |ρ − ρ | > 0 of the Green's function, it is hard to evaluate the convolution integral without adding up a considerable amount of noise.An alternative of lower complexity and much higher precision is to solve for the potential directly using a hybrid boundary-element/finite-element method, commonly referred to as the Fredkin-Koehler method 22 , which only requires computing O(n 4/3 ) matrix elements.In the following, we will expand this method to plane waves.Let us note that, in the finite-difference case of rectangular elements, the plane-wave dipolar tensor can be obtained directly as a dense matrix, as was done for the first time by Henry et al. in Ref. 19.
Typically, if one is to solve for the potential using FEM, the boundary conditions for the potential at infinity need to be specified.To avoid having to model a large air box, these boundary conditions can be mapped onto the cross-section boundary ∂A using the Fredkin-Koehler method for plane waves.This allows to solve for the potential by integrating over the sample cross-section area, only.
Within this section, we drop the indices ν and k for the sake of visual clarity.The potential subject to the screened Poisson equation Eq. ( 28) has to be equipped with the following continuity and jump conditions at the boundary of the magnetic element and where ψ in and ψ out are the lateral potential inside and outside of the magnetic specimen, respectively.These two conditions follow directly from the properties of the full volumetric potential φ(r, t) when the surface-normal vector field n(ρ) ⊥ e z .Furthermore, we require that the potential vanishes as |ρ| → ∞.
The motivation behind the Fredkin-Koehler method is to map the Dirichlet boundary conditions at infinity to the boundary of the magnetic sample by splitting the potential into two parts ψ = ψ 1 + 2 .The first potential ψ 1 satisfies an inhomogeneous Neumann problem inside the magnetic volume and is zero outside of the specimen such that it produces the jump in the normal derivative according to Eq. ( 32).The second potential ψ 2 solves a homogeneous Dirichlet problem with boundary conditions obtained from ψ 1 at the boundary ∂A such that it compensates the jump of ψ 1 and therefore guarantees the continuity of the full lateral potential ψ according to Eq. ( 31) (see Fig. 2(a)).In summary, we have the system of equations Here, u(ρ) is the Dirichlet boundary condition for ψ 2 .After having calculated ψ 1 numerically using FEM, this boundary condition can be obtained by Here Φ(ρ) is the angle subtended by the boundary point ρ within an arbitrary cross section (see Fig. 2(b)), which is π for a smooth boundary and ds is the line element on ∂A.The Dirichlet boundary condition in Eq. ( 36) is the plane-wave version of the initial relation for the regular three-dimensional Poisson problem (equation 14 in the seminal paper by Frekdin and Koehler 22 ) and is rigorously derived in the supplementary material).In a numerical implementation, it can be expressed using a dense-matrix multiplication in the sense where ψ 1,2 are the mesh vectors of the two potentials at the boundary and Bk is a dense matrix.Above, we retained the index k at the Green's function G k (ρ, ρ ) and the dense matrix Bk to highlight the fact that both are wave-vector dependent.Let us note, that this method only requires to calculate O(n 2 B ) matrix elements per wave vector (with n B being the number of boundary nodes) instead of O(n 2 ) when storing the whole Green's function.As a result, the presented method is most effective for cross-section meshes with a large surfaceto-boundary-node ratio.
A note has to be made on the numerical evaluation of the boundary integral in Eq. ( 36).As a specific property of the inhomogeneous Neumann problem for ψ 1 , the first potential in the plane-wave Fredkin-Koehler method diverges as k → 0 (see supplementary material).This divergence is counteracted by the second potential ψ 2 such that the total potential ψ remains finite.Thus, the correct treatment of the dipolar interaction in our method is very sensitive on the proper calculation of the boundary values, i.e. on the dense matrix Bk .For this reason, we recommend at least a sixth-order segmentation of the boundary elements to evaluate the boundary integral, weighted linearly on the 2 6 + 1 resulting nodes.Due to the reasonably small number of boundary nodes and elements the dense matrices are in general small in storage requirements and can be computed on the fly for the different k values.

V. APPLICATIONS
To showcase and validate our method, we calculate the spin-wave dispersion for three different examples and compare the results for different methods or theoretical predictions typically used to investigate spin-wave dynamics.We do not intend to discuss the physical implications of the different examples in detail, but rather highlight the capabilities of our numerical approach.

A. Longitudinally magnetized rectangular waveguide
As a first example, we calculate the spin-wave dispersion for a standard magnonic problem, namely a rectangular waveguide with fixed thickness T = 29 nm and different widths W, where the equilibrium magnetization m 0 is aligned along the propagation (z) direction (see Fig. 3(a)).In such a case, the spin waves exhibit a negative slope of the dispersion ω ν (k) for small wave vectors, and are typically referred to as backward-volume modes.In the case of a waveguide with finite width, the dispersion is split into multiple branches ν, which belong to standing waves along the width of the waveguide (see Fig. 3(b)).For small aspect ratios T /W 1, their dispersion is analytically described using the theory of Kalinokos and Slavin 1 with effective dipolar pinning conditions introduced by Guslienko et al. 4 .These pinning conditions lead to an effective width of the waveguide W eff which is typically larger than the physical width W , resulting in increased wavelengths along the width λ ν = 2W eff /(ν + 1).To compare our results with the experiments of Roussigné et al. 34 , we apply an external field of µ 0 H z = 55 mT, set a saturation magnetization of M s = 621 kA m −1 and a reduced gyromagnetic ratio of γ/2π = 29.76GHz T −1 .The exchange stiffness is set to a typical value for Ni 80 Fe 20 (permalloy) of A ex = 13 pJ m −1 .The triangular mesh used has an average edge length of 5 nm.
In Figs.3(c-e), we show the dispersion calculated using our eigensolver for different widths of the waveguide W for the lowest four branches ν = 0, 1, 2, 3, overlayed with theoretical and experimental data.Overall, we achieve a very good agreement with the experimental data from Ref. 34.As expected, the agreement with the theoretical calculations (performed according to Refs.1,4) increases as the aspect ratio of the waveguide decreases.For the smallest aspect ratio, i.e. the largest width W = 1.5 µm, we show the real part of the numerically calculated out-of-plane component Re(η z ) of the spatial mode profiles at k = 0 as line scans along the width (x) direction, together with the theoretical prediction.

B. Edge shape in transversally magnetized rectangular waveguide
After having showcased a regular rectangular cross section, we now present the transition to curved geometries.For this purpose we calculate the spin-wave dispersion for a transversally magnetized waveguide of 50 nm thickness and 256 nm width, once with sharp corners and once with round corners, as shown in Fig. 4(a).In the case of a transversally magnetized stripe, the equilibrium magnetization m 0 is perpendicular to the propagation direction and the spin waves, also referred to as magnetostatic surface waves, exhibit a positive dispersion at small wave numbers.We set the saturation to M s = 796 kA m −1 , exchange stiffness to A ex = 13 pJ m −1 and reduced gyromagnetic ratio to γ/2π = 28 GHz T −1 .An external field of µ 0 H x = 600 mT is applied along the width (x) direction.At this field, the waveguide is almost completely saturated except for small regions at the sides (known as flower state).The equilibrium distribution m 0 (ρ) is found by minimizing the energy length density g(m 0 ) = G(m 0 )/L using a conjugategradient method.
We compare our results for the waveguide with sharp edges to a numerically calculated dispersion assuming the same material parameters, obtained using the GPU accelerated micromagnetic solver MuMax 312 which employs a finite-difference method to solve the Landau-Lifshitz-Gilbert equation Eq. ( 4) on a rectangular grid.Spin waves were excited using an out-of-plane oscillating field which is homogeneous along the width (for details, see supplementary material).As a consequence, only modes with an even (symmetric) spatial profile along the width can be excited.Modelling a waveguide with round edges is computationally much more expensive in a finite-difference code such as MuMax 3 since it requires a much finer discretization along the thickness of the sample.In our finite-element dynamic-matrix approach, both cross sections require almost the same computational effort.
In Fig. 4(b), we show the Fourier magnitude P (k, ω) obtained with MuMax 3 as a heatmap, overlayed with the the corresponding dispersions (sharp and round edges) obtained with our eigensolver.Considering the different techniques of discretization, we achieve a very good agreement between both methods.Notably, multiply avoided level crossings between the various dispersion branches are recovered (see also Figs. 4(c)).We see that round edges at the waveguide sides mainly result in a shift of the overall frequencies which can be a blue-or a red shift depending on the branch.However, as seen in 4(d), round edges can have a dramatic impact on the avoided crossings between the branches.
Let us repeat here an another benefit of the propagatingwave dynamic-matrix approach over the full threedimensional approach used by MuMax 3 .Since, in the latter case, wave vectors are obtained by means of spatial Fourier transform, the length of the waveguide has to be increased in order to achieve more wave-vector precision.Depending on the problem at hand, this may result in considerably more computational time.In the plane-wave approach used here, the wave vector simply appears only as a parameter in the resulting eigenvalue problem and, thus, can be varied continuously.

C. Magnetic nanotubes with easy-plane anisotropy
As a last example, we go fully into the field of curvilinear magnonics and calculate the spin-wave dispersion of a magnetic nanotube with thin mantle magnetized in the vortex state (see the inset in Fig. 5(a)).We consider a tube with 60 nm outer and 40 nm inner diameter.At such a size and in the absence of any external field, a vortex state can be stabilized using an easy-plane anisotropy which is, in our case, a uniaxial anisotropy along the z axis with negative constant K u = −50 kJ/m 3 .Other than this, we use the same permalloy material parameters as in the previous example.
The spin-wave dispersion in such a system exhibits a curvature-induced asymmetry (ω ν (k) = ω ν (−k)) of purely dipolar origin (see Otálora et al. 6,7 ).For the case of a thin mantle, the lateral profiles can be described approximately as homogeneous along the thickness ρ and being proportional to η ν ∝ exp(iνϕ) with ϕ being the azimuthal angle with respect to the z axis.Here, the lateral mode index ν counts the periods along the azimuthal direction and can take positive and negative values.In the vortex state, modes with opposite ν are degenerate, i.e their frequencies are equal ω ν = ω −ν .In Fig. 5(a), we compare the asymmetric dispersion calculated using our eigensolver with the theoretical predictions by Otálora, showing an almost perfect agreement.Deviations in the higher-order azimuthal modes ν = ±3 are of exchange origin and can be decreased by further decreasing the discretization of the mesh, which in our case was 3 nm.Finally, Fig. 5(b) shows the real part of z component Re(η z ) of the numerically obtained spatial mode profiles which correspond, as expected, to harmonic waves running in azimuthal direction.Whereas, here for the case of tubes with thin mantle, our numerical approach serves to confirm a theoretical model, it may become useful e.g. for the study of magnetic nanotubes with thick mantle for which theoretical models are not available at the time and full three-dimensional dynamic simulations become computationally demanding.As a comment on the performance of our propagatingwave finite-element approach: We obtained the same dispersion as in Fig. 5(a) -without being able to separate degenerate branches -using the highly-optimized GPU-accelerated three-dimensional finite-element code TetraMag 11 .This was done by using a mesh with similar number of elements within the cross section, but extruded to a certain length in propagation direction to give an adequate wave-vector resolution.Even with the efficient field-pulse method described in the supplementary material and exciting several azimuthal branches up to |ν| = 3 at the same time the computation (including post-processing) took several days on a high-performance TITAN Xp GPU.Using the also GPUaccelerated finite-difference code MuMax 3 requires a similar amount of time due to the large number of rectangular cells needed in the cross section to adequately model the curved surface of the tube.For comparison, obtaining the dispersion and the corresponding mode profiles with a higher wave-vector resolution using an unoptimized implementation of our finite-element propagating-wave dynamic-matrix approach on a standard laptop CPU took less than 15 minutes.

VI. CONCLUSION
In summary, we have presented a finite-element propagating-wave dynamic-matrix approach to efficiently calculate the spin-wave dispersion in waveguides with arbitrary (bounded) cross section and translationally invariant magnetic equilibrium along the propagation direction.This was achieved by numerically solving a plane-wave version of the linearized equation of motion of the magnetization.In contrast to dynamic micromagnetic simulations, spin-wave frequencies and mode profiles are obtained without any post-processing.Important characteristics of the spin waves such as linear damping rate or dynamic susceptibility can be calculated from the spatial mode profiles (see e.g.Ref. 29  and Ref. 30, respectively).Our approach differs from the finite-difference approach presented in Ref. 19 mainly in the discretization type and in the way the dipolar potential of the spin-wave modes is calculated at a given wave vector.For this purpose we presented an extension to the Fredkin-Koehler method to the screened Poisson equation for propagating waves, which is suitable to model waveguides with arbitrarily shaped bounded cross section.For a number of known systems, our results were validated using theoretical predictions, dynamic micromagnetic simulations and experimental results.Although, in the magnetic interactions considered, we have restricted ourselves to uniaxial anisotropy as well as dipolar and exchange interaction, other contributions like the Dzyaloshinskii-Moriya interaction can be included in the same way.Due to the possibility to model waveguides with surface curvature, we believe that our approach is of particular relevance for the emerging field of curvilinear magnetism but also useful for standard problems in magnetization dynamics.

SUPPLEMENTARY MATERIAL
See supplementary material which contains the derivation of the Dirichlet boundary condition Eq. (36) within the Frekdin-Koehler method for plane waves.Moreover, we show how the divergence of the first lateral potential ψ 1 is canceled out by the second one ψ 2 .Finally, the dynamic micromagnetic simulations used for comparison in Sec.V are described.

Figure 1 .
Figure 1.(a) Schematics of an infinitely extended waveguide with arbitrary cross section A and translationally invariant equilibrium magnetization m0.(b) Illustration of the transformation R from lab system {x, y, z} into local coordinate system {e1, e2, e3} attached to the equilibrium magnetization.

Figure 2 .
Figure 2. (a) Schematics of the different lateral potentials in the Fredkin-Koehler method for plane waves, shown inside and outside of the magnetic sample when crossing the boundary ∂A.For simplicity, we only show the real parts of the potentials here.(b) Definition of the angle subtended at a certain boundary node point.

Figure 3 .
Figure 3. (a) Schematics of the longitudinally magnetized rectangular waveguide overlayed with the equilibrium magnetization.(b) Numerically and theoretically calculated spatial mode profiles of the first four dispersion branches ν at k = 0 for a width of the waveguide W = 1.5 µm, shown as line scans across the width of the waveguide.(c-e) Dispersion of the same branches for different widths of the waveguide, obtained using our eigensolver (solid lines) and compared with theoretical prediction (dashed) calculated according to Refs.1,4 and experimental data (rhombus) from Ref. 34.

Figure 4 .
Figure 4. (a) Schematics of the two different waveguide cross sections overlayed with the equilibrium magnetizations.(b) Dispersion map (wave-vector-and frequency-dependent Fourier magnitude) for the waveguide with sharp edges obtained with MuMax 3 and shown as a colormap.Overlayed are the dispersion curves for sharp (solid) and round (dotted) corners, obtained with our eigensolver.(c,d) show zoom-ins of the same data, highlighting the avoided level crossings of the different dispersion branches as well as the influence of the edge roundness on the nature of the crossings.

Figure 5 .
Figure 5. (a) Numerically (solid) and theoretically (dashed) obtained asymmetric spin-wave dispersion in a vortex-state magnetic nanotube with easy-plane anisotropy.The branches correspond to the first five azimuthal mode indices around ν = 0, while modes with opposite index ±ν are degenerate.The inset shows the cross section of the nanotube overlayed with the equilibrium magnetization.(b) Numerically obtained spatial mode profiles (here only the real part of the z component ηz is shown).