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

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,⁸ and narrow ribbons⁹ 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.¹⁹ 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²⁰ and those propagating along Bloch domain walls.²¹ 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,²² 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.

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 μ₀ 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.²⁵ 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.²⁶ 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 = γμ₀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 · ∇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, composed of the internal and the external fields. This (unitless) effective field is obtained from Gibb’s free energy as the variational derivative,²⁷ 

The second term T_G, in torque equation (4), represents viscous (Gilbert) damping and relaxes the magnetization to equilibrium after a certain time.²⁸ 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.²⁹ 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₀(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 |m| = 1, the equilibrium direction and the dynamical component must be orthogonal, m₀ ⊥ δm. 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

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ν*$⁠. Here, the asterisk denotes complex conjugation.

The Hamiltonian operator $Ω̂$ is a tensor operator given by

where h₀ = m₀ · h_eff(m₀) is the projection of the unitless static effective field (including any static external field) onto the equilibrium direction m₀ and $Î$ is the identity operator. The operator $Ω̂$ is self-adjoint and, as long as m₀ is a local minimum of Gibb’s free energy $G$⁠, it is also positive definite for vectors m_ν perpendicular to m₀. 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 $Ω̂$ is made dimensionless.

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₀(ρ) 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₀(ρ) 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 $Ω̂$ is positive definite, it follows that $Ω̂k$ is too, i.e., for some arbitrary lateral profile v = e^−ikzw that satisfies v ⊥ m₀ (with w being the corresponding full volumetric profile)

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 − ω_νkt)], 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).

Recall that we want to solve the plane-wave eigenvalue problem [Eq. (11)] for eigenvectors orthogonal to the equilibrium magnetization, m₀ ⊥ η_ν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₀ or by rotating the eigenvalue problem into the local reference frame of m₀ and reducing the dimension of the eigenvectors.¹⁴ The rotation is performed using a suitable unitary transformation $R̂$ into a local orthonormal system {e₁, e₂, e₃}, where m₀ = e₃ 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 dynamic-matrix operator

and the dagger denoting Hermitian conjugation. Due to the rotation of the mode profiles into the local reference frame of m₀, 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., $η̃ν⋅e3≡0$⁠. 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₀. In the local reference frame, the cross product m₀ × ⋯ takes the simple form of a multiplication with the matrix (written in {e₁, e₂} basis)

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̂k$ of the magnetic tensor for the different interactions considered. For each interaction, the contribution to the equilibrium-field projection h₀ can of course be obtained by applying the operator $N̂k$ to the equilibrium configuration m₀ at k = 0. For example, if only an external field and exchange interaction would be considered, one would have

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 $Ω̂k(uni)$ does not change its form when transformed to the waveguide cross section, i.e., $N̂(uni)=N̂k(uni)$ is not dependent on the wave vector.

The exchange interaction within the continuum limit is written as

with $λex=2Aex/μ0Ms2$ being the exchange length, A_ex the exchange stiffness constant of the material, and ∇² the Laplacian acting on three-dimensional vector fields in three-dimensional space. Projecting the operator into the waveguide cross section yields

Here, $∇ρ2$ denotes the part of the Laplacian operator that operates on a three-dimensional vector fields in the xy plane.

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

For a numerical solution of the eigenvalue problem in Eq. (16), the equilibrium magnetization m₀(ρ), mode profiles η_νk(ρ), potentials ψ_νk(ρ), and 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 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.¹⁸ 

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 (28) with the appropriate Green’s function

using the modified Bessel function of second kind and zeroth order K₀. On a discrete mesh with n nodes, this approach requires the computation of $O(n2)$ matrix elements of Green’s function.²² 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,²² 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 ψ = ψ₁ + ψ₂. The first potential ψ₁ 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 ψ₂ solves a homogeneous Dirichlet problem with boundary conditions obtained from ψ₁ at the boundary ∂A such that it compensates the jump of ψ₁ 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 ψ₂. After having calculated ψ₁ 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 ds′ 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²²] 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 $B̂k̲$ is a dense matrix. As mentioned above, we retained the index k at Green’s function G_k(ρ, ρ′) and the dense matrix $B̂k̲$ 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 ψ₁, 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 ψ₂ 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̂k̲$⁠. 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⁶ + 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.

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.

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₀ 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¹ with effective dipolar pinning conditions introduced by Guslienko et al.⁴ 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 λ_ν = 2W_eff/(ν + 1). To compare our results with the experiments of Roussigné et al.,³⁶ we apply an external field of μ₀H_z = 55 mT and set a saturation magnetization of M_s = 621 kA m⁻¹ and a reduced gyromagnetic ratio of γ/2π = 29.76 GHz T⁻¹. The exchange stiffness is set to a typical value for Ni₈₀Fe₂₀ (permalloy) of A_ex = 13 pJ m⁻¹. 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(η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 · ∇m = 0, which favors completely unpinned waves.

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₀ 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⁻¹, exchange stiffness to A_ex = 13 pJ m⁻¹, and reduced gyromagnetic ratio to γ/2π = 28 GHz T⁻¹. An external field of μ₀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₀(ρ) 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,¹² 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.

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³. 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(η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.¹¹ 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.

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.

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 ψ₁ is canceled out by the second one ψ₂. Finally, the dynamic micromagnetic simulations used for comparison in Sec. V are described.

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.

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