We present a numerical approach to efficiently calculate spin-wave dispersions and spatial mode profiles in magnetic waveguides of arbitrarily shaped cross section with any non-collinear equilibrium magnetization that is translationally invariant along the waveguide. Our method is based on the propagating-wave dynamic-matrix approach by Henry *et al.* (Ref. 19) and extends it to arbitrary cross sections using a finite-element method. We solve the linearized equation of motion of the magnetization only in a single waveguide cross section, which drastically reduces computational effort compared to common three-dimensional micromagnetic simulations. In order to numerically obtain the dipolar potential of individual spin-wave modes, we present a plane-wave version of the hybrid finite-element/boundary-element method by Fredkin and Koehler which we extend to a modified version of the Poisson equation. Our method is applied to several important examples of magnonic waveguides including systems with surface curvature, such as magnetic nanotubes, where the curvature leads to an asymmetric spin-wave dispersion. In all cases, the validity of our approach is confirmed by other methods. Our method is of particular interest for the study of curvature-induced or magnetochiral effects on spin-wave transport and also serves as an efficient tool to investigate standard magnonic problems.

## I. INTRODUCTION

Over the last few 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 analytical expressions for the dispersion relation and spatial mode profiles of spin waves with wavelengths in the nano- and micrometer range. These models allow us to describe the spin-wave propagation in simple geometries, such as bulk systems, infinitely extended mono- or bilayers,^{1–3} or waveguides with a 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} has been described.

Although the developed analytical 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–12} In such dynamic simulations, spin waves are typically obtained by exciting the magnetic system with a microwave field pulse, letting the system evolve according to the appropriate equation of motion, and 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 us 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 a 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 finite-element and finite-difference methods of discretization.^{15–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, modeling 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 that can become extremely costly when modeling a full three-dimensional magnetic specimen.

Although the finite-difference-based propagating-wave dynamic-matrix approach is very efficient and allows us to study, e.g., complex multilayer systems, it is not suitable for waveguides with arbitrary (bounded) cross sections, 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 that allows us 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 finite-element/boundary-element method to solve the Poisson equation on a finite volume. Here, 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 and also allows us to calculate spin-wave dispersions in various standard magnonic problems.

## II. THEORETICAL MODEL

### 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}, Gibb’s free energy $G(m)$ of a magnetic specimen is given by

with *μ*_{0} being the vacuum permeability, *V* being the volume of the magnetic body, *h*_{ext} = *H*_{ext}/*M*_{s} being the unitless external field, and *h*_{int} being 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\u0302$, 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, *γ* being the modulus of the gyromagnetic ratio in the material at hand, and finally *T*_{G} being a (for this work unimportant) damping term. The equation of motion must be equipped with the boundary condition ** n** · ∇

**= 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**

*m**effective*field

*h*_{eff}, composed of the internal and the external fields. This (unitless) effective field is obtained from Gibb’s free energy as the variational derivative,

^{27}

The second term *T*_{G}, in torque equation (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–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, and 31, the nonlinear torque equation (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 small-amplitude variation from the equilibrium. Since the norm of the total magnetization vector is constrained by |

**| = 1, the equilibrium direction and the dynamical component must be orthogonal,**

*m*

*m*_{0}⊥

*δ*

**. Inserting the separation equation (6) into the torque equation (4), neglecting all nonlinear terms, and finally expanding the variable magnetization into linear waves according to**

*m*lead to the linearized Landau–Lifshitz equation

which is a standard eigenvalue problem that yields the spin-wave mode frequencies *ω*_{ν} as well as their spatial profiles *m*_{ν}. It is clear from the demand that *δ*** m** is a real vector and that if

*ω*

_{ν}is an eigenvalue of Eq. (8) with corresponding eigenvector

*m*_{ν}, then so is −

*ω*

_{ν}with eigenvector $m\nu *$. Here, the asterisk denotes complex conjugation.

The Hamiltonian operator $\Omega \u0302$ is a tensor operator given by

where *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 $I\u0302$ is the identity operator. The operator $\Omega \u0302$ is self-adjoint and, as long as *m*_{0} is a local minimum of 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 (8) in a similar form to, e.g., Refs. 16, 30, and 31 with the only difference that our operator $\Omega \u0302$ 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 coordinate

**. Moreover, it may also contain a non-zero component along the propagation direction,**

*ρ**m*

_{0,z}≠ 0. We consider a

*long*waveguide with constant cross-sectional 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 plane-wave problem and solve it only in a single cross section of the waveguide. Of course, for this case, 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(m0)=G(m0)/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 the *z* direction and *η*_{νk} the lateral mode profiles that are complex vector fields in the *xy* plane. 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 $\Omega \u0302$ is positive definite, it follows that $\Omega \u0302k$ is too, i.e., for some arbitrary lateral profile ** v** =

*e*

^{−ikz}

**that satisfies**

*w***⊥**

*v*

*m*_{0}(with

**being the corresponding full volumetric profile)**

*w*with respect to the $L2(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 $\eta \nu 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 that 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\u0302$ 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\u0302$ and using its unitarity, i.e., $R\u0302\u2020R\u0302=I\u0302$, one obtains the rotated eigenvalue problem

with the rotated eigenvectors $\eta \u0303\nu k=R\u0302\eta \nu k$, the dynamic-matrix 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., $\eta \u0303\nu \u22c5e3\u22610$. This allows us 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 $N\u0302k$ 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 $N\u0302k$ 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 $\Omega \u0302k(uni)$ does not change its form when transformed to the waveguide cross section, i.e., $N\u0302(uni)=N\u0302k(uni)$ is not dependent on the wave vector.

### B. Symmetric exchange interaction

The exchange interaction within the continuum limit is written as

with $\lambda ex=2Aex/\mu 0Ms2$ being the exchange length, *A*_{ex} the exchange stiffness constant of the material, and ∇^{2} the Laplacian acting on three-dimensional vector fields in three-dimensional space. Projecting the operator into the waveguide cross section yields

Here, $\u2207\rho 2$ denotes the part of the Laplacian operator that 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) and (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) that, 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}(

**), potentials**

*ρ**ψ*

_{νk}(

**), and all involved operators in the dynamic matrix $D\u0302$ 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\u0302$ (except the dipolar operator) take the form of sparse matrices and their actions are simply matrix–vector or matrix–matrix 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}Note that the mass-lumping technique mentioned above leads to a small non-Hermitian contribution of the discretized effective field operators in the dynamic matrix.

^{34,35}However, according to Bruckner

*et al.*, this effect is negligible for typical problems.

^{18}

The evaluation of the dipolar interaction requires special care. The lateral potential *ψ*_{νk} in the dipolar operator $N\u0302(dip)$ could be readily obtained by the convolution of the right-hand side of the screened Poisson equation (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(n2)$ matrix elements of Green’s function.^{22} Moreover, due to the logarithmic singularity at |** ρ** −

**′| = 0 and the rapid decay for |**

*ρ***−**

*ρ***′| > 0 of 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)$ matrix elements for two-dimensional meshes. 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 by Henry *et al.* in Ref. 19.

Typically, if one is to solve for the potential using the 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-sectional boundary *∂A* using the Fredkin–Koehler method for plane waves. This allows us to solve for the potential by integrating over the sample cross-sectional 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 (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 the 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 d

*s*′ is the line element on

*∂A*. The Dirichlet boundary condition in Eq. (39) is the plane-wave version of the initial relation for the regular three-dimensional Poisson problem [Eq. (14) in the seminal paper by Fredkin 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 $\psi 1,2\u0332$ are the mesh vectors of the two potentials at the boundary and $B\u0302k\u0332$ is a dense matrix. As mentioned above, we retained the index *k* at Green’s function *G*_{k}(** ρ**,

**′) and the dense matrix $B\u0302k\u0332$ to highlight the fact that both are wave-vector dependent. Let us note that this method only requires us to calculate $O(nB2)$ matrix elements per wave vector (with**

*ρ**n*

_{B}being the number of boundary nodes) instead of $O(n2)$ when storing the whole Green’s function. As a result, the presented method is most effective for cross-sectional meshes with a large surface-to-boundary-node ratio.

A note has to be made on the numerical evaluation of the boundary integral in Eq. (39). 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 the 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 $B\u0302k\u0332$. 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 Kalinikos 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} that is typically larger than the physical width *W*, resulting in increased wavelengths along the width *λ*_{ν} = 2*W*_{eff}/(*ν* + 1). To compare our results with the experiments of Roussigné *et al.*,^{36} we apply an external field of *μ*_{0}*H*_{z} = 55 mT and set a saturation magnetization of *M*_{s} = 621 kA m^{−1} and a reduced gyromagnetic ratio of *γ*/2*π* = 29.76 GHz 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)–3(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, overlaid with theoretical and experimental data. Overall, we achieve very good agreement with the experimental data from Ref. 36. As expected, the agreement with the theoretical calculations (performed according to Refs. 1 and 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(\eta z)$ of the spatial mode profiles at *k* = 0 as line scans along the width (*x*) direction, together with the theoretical prediction. At the edges of the waveguide, there is a visible deviation of the numerical profiles from the analytical prediction, best seen for the highest-order mode ** ν = 3**. This can be explained by the fact that the effective theory of Guslienko

*et al.*does not consider the exchange boundary condition

**· ∇**

*n***= 0, which favors completely unpinned waves.**

*m*### 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(m0)=G(m0)/L$ using a conjugate-gradient method.

We compare our results for the waveguide with sharp edges to a numerically calculated dispersion assuming the same material parameters, obtained using the graphics-processing-unit-(GPU)-accelerated micromagnetic solver MuMax3,^{12} which employs a finite-difference method to solve the Landau–Lifshitz–Gilbert equation (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 the supplementary material). As a consequence, only modes with an even (symmetric) spatial profile along the width can be excited. Modeling a waveguide with round edges is computationally much more expensive in a finite-difference code such as MuMax3 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 MuMax3 as a heatmap, overlaid with the correspondingdispersions (sharp and round edges) obtained with our eigensolver. This magnitude *P* is obtained by first performing a fast Fourier transform (FFT) of the magnetization in time and along the *z* direction and then averaging the resulting spectra in the *xy* plane.

Considering the different techniques of discretization, we achieve very good agreement between both methods. Notably, multiple avoided level crossings between the various dispersion branches are recovered [see also Fig. 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 Fig. 4(d), round edges can have a dramatic impact on the avoided crossings between the branches.

Let us repeat here another benefit of the propagating-wave dynamic-matrix approach over the full three-dimensional approach used by MuMax3. 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 diameters. 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 Refs. 6 and 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(\eta z)$ of the numerically obtained spatial mode profiles that correspond, as expected, to harmonic waves running in the 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 propagating-wave 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 a similar number of elements within the cross section, but extruded to a certain length in the 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 GPU-accelerated finite-difference code MuMax3 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 central processing unit (CPU) took less than 15 min.

## 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., Refs. 29 and 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 and also useful for standard problems in magnetization dynamics.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the derivation of the Dirichlet boundary condition [Eq. (39)] within the Fredkin-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.

## ACKNOWLEDGMENTS

The authors are very thankful to Burkhard Clauß, Steffen Boerm, and Jorge A. Otálora for fruitful discussions. Financial support from the Deutsche Forschungsgemeinschaft within the programs under Grant Nos. KA 5069/1-1 and KA 5069/3-1 is gratefully acknowledged.

## DATA AVAILABILITY

The data that support the findings of this study are openly available in RODARE at https://doi.org/10.14278/rodare.956.^{37}

## REFERENCES

The term *internal* here is meant in a thermodynamic sense (as in *internal energy*) and not in a spatial sense, because obviously the external field is also present within the volume of the sample. In a spatial sense, the internal field would be the effective field defined in Eq. (5).

The negative sign in the internal field is usually chosen by convention in accordance with the demagnetizing/dipolar field.