Finite-element dynamic-matrix approach for propagating spin waves: Extension to mono- and multilayers of arbitrary spacing and thickness

In our recent work [AIP Adv. 11, 095006], we presented an efficient numerical method to compute dispersions and spatial mode profiles of spin waves propagating in waveguides with translationally invariant equilibrium magnetization. Using a finite-element method (FEM) allowed to model two-dimensional waveguide cross sections of arbitrary shape but only finite size. Here, we extend our FEM propagating-wave dynamic-matrix approach from finite waveguides to the important practical cases of infinitely-extended mono- and multilayers of arbitrary spacing and thickness. To obtain the mode profiles and frequencies, the linearized equation of motion of magnetization is solved as an eigenvalue problem only on a one-dimensional line-trace mesh, defined along the normal direction of the layers. Being an important contribution in multilayer systems, we introduce interlayer-exchange interaction into our FEM approach. With the calculation of dynamic dipolar fields being the main focus of this paper, we also extend the previously presented plane-wave Fredkin-Koehler method to calculate the dipolar potential of spin waves in infinite layers. The major benefit of this method is that it avoids the discretization of any non-magnetic material, such as non-magnetic spacers in multilayers. Therefore, the computational effort becomes completely independent on the spacer thicknesses. Furthermore, it keeps the resulting discretized eigenvalue problem sparse, which therefore, inherits a comparably low arithmetic complexity. As a validation of our method (implemented into the open-source finite-element micromagnetic package TetraX), we present results for various systems and compare them with theoretical predictions as well as with established finite-difference numerical methods. We believe this method offers an efficient and versatile tool to calculate spin-wave dispersions in layered magnetic systems.


I. INTRODUCTION
Over the last decades, micromagnetic simulations have become a powerful method to predict new effects in magnetism besides being used to complement experimental results and to validate analytical derivations.The continuous requirement to simulate the static magnetic states and magnetization dynamics in complex systems has lead to the development of highly-optimized micromagnetic codes [2][3][4][5] for classical time-integration-based methods that solve the equation of motion of the magnetization on a discrete mesh.Certainly, these codes have already played a key role in the investigation of the fundamental as well as applied aspects of magnetism, such as magnetization reversal processes, 6,7 domain wall dynamics, [8][9][10][11][12] skyrmion motion driven by spin-polarized currents, 13,14 non-linear magnetization dynamics for reservoir and neuromorphic computing, [15][16][17][18][19] or the hunt for exotic magnetic textures such as hopfions, [20][21][22][23] just to name a few.However, even with such micromagnetic codes, the numerical study of spin-wave propagation in waveguides or extended films remained ineffective due to the required long simulation times and large mesh dimensions for proper frequency and wave-vector resolution.Other pitfalls are the extensive post processing required after each simulation, the necessity * ) Electronic mail: l.koerber@hzdr.de to have a priori knowledge about the mode-profile symmetries, the inability to detect degenerate modes, and other related problems. 1,24This, finally, initiated the development of different approaches.Among those best suited to study linear spin-wave dynamics is the dynamic-matrix method, [25][26][27][28] especially the approach for propagating spin waves in systems with a translationally invariant magnetic equilibrium. 24,29n contrast to standard time-domain micromagnetic simulations, which rely on the time-integration of the nonlinear Landau-Lifshitz-Gilbert (LLG) equation of motion of the magnetization and subsequent post-processing of the simulation data by means of Fourier analysis to obtain spin-wave frequencies and spatial mode profiles, dynamic-matrix approaches (or dynamic-matrix methods) are based on an exact numerical solution of the linearized LLG equation about some (unitary) equilibrium magnetization m 0 (r).The linearized equation of motion for the complex and unitless spatial mode profiles δm(r) of the spin-wave modes reads as with ω M = γµ 0 M s being the characteristic magnetic frequency, γ being the modulus of the gyromagnetic ratio and M s being the saturation magnetization of the magnetic body at hand, the vectors h 0 and δh are the unitless static and dynamic effective fields.Considering a waveguide with at least one of its dimensions extended to infinity, the linearized equation can be trans-arXiv:2207.01519v2[cond-mat.mes-hall]5 Jul 2022 formed into a plane-wave problem.The spin-wave mode profiles for a magnetic waveguide with translationally invariant equilibrium can be written as δm ∝ m k ≡ η k (x, y)e i(kz−ωt) , (2)   with k being the wave number in the direction of propagation (z in this manuscript) and η k being the complex lateral mode profile.Using this definition and after some transformations (see e.g.Ref. 1), the linearized equation Eq. ( 1) becomes a wave-vector-dependent eigenvalue problem that yields the frequencies of the eigenmodes ω(k) and their lateral spatial profiles η k simultaneously for wavelengths far greater than the lattice parameter of the investigated magnetic specimen.Here, Dk is the so-called dynamic matrix which, in the case of propagating waves, depends on the wave vector k.
After spatial discretization, the dynamic matrix Dk can be diagonalized for each k using a suitable numerical eigensolver.Apart from arbitrary frequency-and wave-vector resolution, a major benefit in adapting the numerical modeling to plane waves is that only the lateral spatial directions (x, y) need to be discretized.This means that, for example, in the case of a waveguide, only a single two-dimensional cross section needs to be modeled.In the case of infinite film, even only a linetrace along its normal direction needs to be modeled. 24,29,30ecently, we presented a finite-element dynamic-matrix approach that allows to calculate the spin-wave dispersion and related mode profiles in translationally invariant waveguides with arbitrary cross section. 1 To take into account the (usually) computationally demanding dipolar fields, we presented an extension of the finite-element/boundary-element method by Fredkin and Koehler 31 to calculate the lateral dipolar potentials of propagating spin waves.This method relies on a numerical solution of the Poisson equation, which governs the dipolar potential of each mode, and comes with the major advantage that only the magnetic material needs to be modeled.Moreover, it keeps the resulting eigenvalue problem Eq. (3) sparse, leading to a lower arithmetic complexity than explicitly calculating the dipolar tensor in matrix form, which is usally done in finite-difference codes.Furthermore, discretizing the cross section of a waveguide using triangular finite elements allowed to model the spin-wave dynamics in waveguides of any cross-section shape, such as polygonal tubes, 32,33 thick-shell round nanotubes, 34 rectangular waveguides 1,35 and many more, as long as the area of the cross section remained finite.
Naturally, the important cases of mono-and multilayer systems cannot be captured by this method.Therefore, in this paper, we extend our finite-element dynamic-matrix approach for propagating waves to infinitely extended mono-or multilayers, which, notably, can have arbitrary thicknesses and spacings between each other.We also note that the equilibrium magnetization and spin-wave dynamics can be inhomogeneous along the thickness of each individual layer.For the extension, we derive analytically the most critical part, namely the computation of the dipolar fields, by extending the plane-wave Frekdin-Koehler method to infinite layers.As this method avoids the discretization of the non-magnetic material, the computational effort is completely independent of the spacer thickness, which solely appears as a numerical factor in the computations.To validate our method, we first compare our results for the dipole-exchange spectra of spin waves in monolayers with varying film thickness with well-known analytical formulae from the literature, as well as with established finite-difference propagating-wave dynamicmatrix approaches. 24,36Furthermore, we validate our method for magnetic bilayer systems, when the ferromagnetic layers are separated by a non-magnetic spacer with varying spacing distance, by comparing the dispersion of the symmetric and antisymmetric modes to the analytical derivations from Gallardo et al. 37 As an important contribution for the study of multilayer systems, we also introduce interlayer-exchange interaction into our FEM dynamic-matrix approach.The presented method is readily implemented in the TETRAX opensource micromagnetic modeling package. 38

II. NUMERICAL IMPLEMENTATION OF EFFECTIVE FIELDS
The construction of the dynamic matrix from Eq. ( 3) requires the evaluation of the magnetic interactions and their related stiffness fields as in usual micromagnetic simulations.Therefore one may need to take into account contributions of the various interactions, such as exchange (being symmetric or asymmetric), dipole-dipole, magnetocrystalline, interlayerexchange and Zeeman interaction.Here, we give a detailed consideration of the dipolar field, as being the main focus of the paper, while also introducing expressions for the symmetric (intralayer) exchange as well as the interlayer-exchange interaction.

Screened Poisson equation of propagating spin waves
Taking into account the effect of dynamic dipolar fields on spin-wave dispersions requires calculating the dipolar field (or demagnetizing field) of each mode as with φ k (r, t) being the unitless dipolar potential (or magnetostatic potential) generated by each mode.The linear operator N(dip) k which outputs the dipolar field at wave-vector k is also referred to as the plane-wave dipolar (or demagnetizing) tensor. 24Considering some magnetic waveguide which is translationally invariant along the z direction [see Fig. 1(a)], the potential of spin waves propagating along this direction with wave number k is of the form Here, ψ k (ρ) is the complex lateral potential of the respective mode, which is defined only within the cross section A of the magnetic element, which, in return, is some arbitrarily shaped subset of the (x, y) plane [see Fig. 1(a)].The lateral potential ψ k (ρ) can be obtained by solving the screened Poisson equation with the following continuity and jump conditions at the boundary of the magnetic element: and where ψ in/out is the potential inside/outside of the magnetic material and n is the normal vector.It is also required that ψ(ρ) → 0 as ρ → ∞.

Recap: Fredkin-Koehler method for propagating spin waves
There are several ways to numerically obtain a solution to the boundary-value problem Eqs.(6).In general, it is possible to calculate the potential by convolution with the kernel which is the Green's function of the Yukawa operator ∆ − k 2 in two dimensions, with K 0 being the modified Bessel function of second kind and zeroth order.Performing this convolution numerically is non-trivial due to the divergent and asymptotic behavior of the Green's function.These features can lead to a considerable accumulation of numerical noise and, thus, to an imprecise calculation of the dipolar fields.
When discretizing a magnetic body using a structured mesh, as commonly done in the finite-difference (FD) method, this integration can indeed be made much more stable by averaging within the regular cells of the mesh and calculating the kernel for each cell.This, in return, allows to explicitly calculate the dipolar tensor N(dip) k for each wave vector k in matrix form.This was done, for example, for plane layer stacks or rectangular waveguides in Ref. 24.Note, however, that this makes the eigenvalue problem Eq. (3) dense due to the longrange character of the dipole-dipole interaction.
When working with unstructured meshes, as done in finite element methods (FEM), this explicit calculation of the cell kernels is not possible.Instead, to circumvent the asymptotic characteristics of the Green's function, in FEM, the dipolar field of each propagating mode can be calculated by first calculating the lateral potential by explicitly solving the screened Poisson equation with according boundary-and jump conditions [Eq.(6)].Note that this keeps the eigenvalue problem Eq. ( 3) sparse. 39Recall that to solve the boundary-value problem Eq. ( 6), the boundary conditions at infinity need to be specified.In order to avoid having to model a large "airbox" around the magnetic sample, these boundary conditions can be mapped directly onto the sample surface using a hybrid finite element/boundary element (FEM/BEM) method known as the Fredkin-Koehler method. 31A major benefit of this method is that only the magnetic sample itself needs to be modeled.In our recent work, 1 this method has been extended to a plane-wave Fredkin-Koehler method capable of calculating the lateral potential of propagating waves by considering only a single (finite) cross section of the waveguide.Before extending this method further from finite cross sections (waveguides) to infinite cross sections (extended layer stacks), for the reader's convenience, it is worth recalling the basic ideas of the plane-wave Fredkin-Koehler method.
To improve the readability, the subscript k is suppressed for the next part of the manuscript.The idea of the Fredkin-Koehler method is to divide the lateral potential into two parts, such that the first potential fulfills the jump condition in Eq. (6c) at the boundary, which is given by the magnetic surface charges.The second potential, coupled to the first one through an appropriate Dirichlet boundary condition, ensures the continuity of the whole potential ψ = ψ 1 +ψ 2 .The role of these two potentials is sketched in Fig. 1(b).To obtain them, the following set of equations must be solved for the first potential and for the second potential where u(ρ) denotes the Dirichlet boundary condition for ψ 2 , which is of the form and can be calculated after ψ 1 is computed by solving the sparse linear system corresponding to Eq. (8a).Here Φ(ρ) is the boundary angle subtended by the boundary point ρ within a cross section [see zoom-in in Fig. 1(a)], and ds is the line element on ∂A.Let us note that Eq. ( 9) is the planewave version of the boundary condition, which can be derived (as done in Ref. 1) directly from the original relation for the regular three-dimensional Poisson problem from Fredkin and Koehler. 31When discretizing the magnetic body into finite elements, all differential operators (∇, ∆ and so forth) take the form of sparse matrices.Furthermore, the Dirichlet boundary condition Eq. ( 9) can be expressed using the Dirichlet matrix B such that with ψ 1,2 being the mesh vectors of the two potentials at the boundary.It is clear, that in order to extend the plane-wave Frekdin-Koehler method to infinitely extended layer stacks, we need to derive a tractable expression for the Dirichlet matrix Bk which can be evaluated as the cross section of the waveguide becomes infinite in one direction.Note that, Bk is a dense matrix of size n B × n B (with n B being the number of boundary nodes).This means that the Frekdin-Koehler method re-introduces a dense-matrix multiplication into the solution of the eigenvalue problem Eq. ( 3).We will see soon, however, that the Dirichlet matrix for extended layer stacks will be very small.For the following discussion it is important to remark that this method can also be used if the magnetic region consists of several disjoint regions.

Extension to mono-and multilayers
To extend the plane-wave Fredkin-Koehler method to mono-and multilayers, we consider a stack of rectangular waveguides with finite width W , as exemplified in Fig. 2(a), which we can describe with the aforementioned method.Let A j be the cross section of the j-th rectangular waveguide, such that the full cross section is given by A = ∪ j A j .We now want to calculate the boundary condition for the second potential ψ 2 [Eq.( 9)] analytically while letting the width of all layers go to infinity, W → ∞ along the x direction, as seen in Fig. 2(b).This allows to model the whole layer stack using only a single one-dimensional (1D) line-trace mesh along its thickness (y) direction, as also seen in Fig. 2(b).Note that, since the Fredkin-Koehler method avoids modeling nonmagnetic regions, the mesh does not have to be connected and the resolution can be freely varied for the different layers.
As a first observation, we see that, as W → ∞ the first lateral potential ψ 1 becomes independent of x and can only depend on the thickness coordinate y, thus ψ 1 = ψ 1 (y).This is clear from the fact that the lateral mode profiles can only depend on the thickness direction in a thick film [η = η(y)] since we set the propagation direction in z direction.Therefore, the Poisson-Neumann problem Eqs.(8a-8c) for ψ 1 becomes translationally invariant along the x direction.
Let us proceed with the calculation of the Dirichlet boundary condition u for ψ 2 from ψ 1 (y).We denote the first term of Eq. ( 9) as the Green's function contribution u G (ρ) and the second term as the boundary-angle contribution u Φ (ρ).The latter is quite trivially u Φ (ρ) = −ψ 1 (y)/2, since the boundary angle for a rectangular element is always Φ(ρ) = π, except at the corners of the rectangle which move towards infinity, as W → ∞.Furthermore, the Green's function contribution separates into a sum over the layers j.With this, we have as the boundary condition.What remains is the calculation of the Green's function u G,j (ρ) contribution of each rectangular waveguide.
We can split the integral along the circumference of each cross section into the different edges as shown in Fig. 2(c), Since the Green's function Eq. ( 7) goes to zero as W → ∞ the integrals over the side facets (left and right) vanish and can, therefore, safely be ignored (see supplementary material for proof).With that, we can change the summation of layers and top-and bottom surfaces to a summation over all remaining surfaces S of the layer stack [see Fig. 2 Executing the normal derivative (using K 0 = −K 1 ) and inserting ρ = (x, y, 0), we have as the contribution of each remaining surface Here n = ±1 is the sign (up or down) of the outward-normal direction of each surface [see Fig. 2(c)].As we let W → ∞, the integral over x becomes translationally invariant, and, therefore, independent of x.Therefore, we can safely set x = 0.With the abbreviation ∆y = y − y the Green's-function contribution u G, takes the form with As carried out in the supplementary material, this integral can be solved in a closed form and one arrives at the tractable expression with "sgn" denoting the sign function, using the convention sgn(0) = 0.As a result, the whole Dirichlet boundary condition takes the very simple form with, again, running over the boundary surfaces (boundary nodes of the 1D mesh) and n being the normal direction along the y axis of each surface.It is important to note that for infinite layers (1D samples), the Dirichlet boundary conditions are calculated without having to perform any numerical integration, like is the case for waveguides with finite (2D) cross section or volumetric (3D) samples.Furthermore, we see that, for infinite layers, the boundary matrix B will always be only of size 2N × 2N (with N being the number of layers) and, therefore, its size will be completely independent of the actual thickness of the layers.This fact has considerable implications on the arithmetic complexity of the plane-wave Frekdin-Koehler method for infinite layers.Take, for example, the case of a single monolayer, where B contains only four elements.With increasing number of nodes n along the thickness of the layer, only the discretized differential operators (∆, ∇ etc.), which are all sparse, increase in size.With that, on a 1D mesh, the number of nonzero elements in these sparse matrices only increases as O(n).In contrast to this, explicitly calculating the matrix elements of the dipolar tensor N(dip) for the same layer and number of cells along the thickness leads to a dense matrix with the number of non-zero elements scaling as O(n 2 ).Si and Sj are the surfaces of the bilayer system adjacent to the interlayer, while ci and cj are the "volumes" (being the length in the case of the 1D line-trace mesh) of the Wigner-Seitz cells associated to the boundary nodes.

B. Interlayer-exchange field
To study spin-wave dynamics in multilayer systems it is often also desirable to model a possible interlayer exchange coupling between the different layers.If two ferromagnetic layers are separated by a metallic non-magnetic interlayer, the exchange coupling due to Ruderman-Kittel-Kasuya-Yosida [40][41][42] (RKKY) interaction between them gives rise to the following bilinear energy where J bl (given in J/m 2 ) is the bilinear interlayer-exchange constant, Γ is the surface of the bilayer system adjacent to the interlayer and P (r) maps the point r on the surface of one magnetic layer to the nearest point on the surface of the other magnetic layer. 43Equation ( 20) can be obtained by replacing the nearest-neighbor sum in a Heisenberg Hamiltonian with an integral.Variation of the energy according to Eq. ( 20) leads to the unitless interlayer-exchange field acting only on the surface Γ as On a 1D line-trace mesh and using finite-element discretization, this field can be obtained as a matrix-vector multiplication using the off-diagonal sparse block matrix N(IEC) with the blocks with Î being the identity on R 3 , i and j being the indices of two coupled boundary elements.Here, on a line-trace mesh, c i and c j are the lengths of the Wigner-Seitz cell associated with each boundary node [see Fig. 3].We note that biquadratic interlayer exchange can be considered in a similar way, by including biquadratic terms in the energy functional in Eq. ( 20), calculating the corresponding field and subsequently linearizing it with respect to some magnetic equilibrium state.

C. Symmetric-exchange field and other effective fields
Before being able to calculate the dispersion of spin waves in infinite layers, we also need to consider the (symmetric) internal exchange field of the individual layers.Starting, for example, from the exchange operator for propagating waves in Ref. 1 and considering the fact that the lateral mode profiles η k in infinitely extended layers can only depend on its thickness (y) coordinate, the unitless lateral exchange field can be obtained as with λ ex = 2A ex /µ 0 M 2 s being the exchange length and A ex being the exchange stiffness constant of the material.In order to assemble the dynamic matrix Dk , the differential operator d 2 /d 2 y needs to be discretized on the 1D line-trace mesh of the layer [see again Fig. 2(b)] using finite elements under consideration of the exchange boundary condition We note that, using finite elements, the second y derivative on the considered line-trace mesh will be quite similar to the finite-difference version obtained by taking central derivatives.Deriving expressions for other magnetic interactions such as asymmetric exchange interaction (of bulk-or interface origin), or magneto-crystalline anisotropies in the case of propagating waves in infinite layers works in an analogous and straightforward way.Therefore, these fields are not presented here.As the main result of this work is the calculation of the dynamic dipolar field in infinite layers using finite elements, it is enough to only consider exchange-, dipolar and interlayerexchange interaction in the following examples.

III. VALIDATION AND APPLICATIONS
In the remaining part of this paper, we want to validate our FEM dynamic-matrix approach for extended layers for a number of different examples.For this we implemented the developed numerical scheme into the TETRAX open-source micromagnetic modeling package 38 and will test it by calculating the spin-wave spectra in different mono-and bilayer systems.For our calculations, we adopt typical material parameters of the soft magnetic alloy Ni 80 Fe 20 as summarized in Tab.I.

A. External-field dependence of uniform modes (ferromagnetic resonance)
As a first example to validate our method, we consider a magnetic monolayer of thickness d.For this very simple case, the Dirichlet matrix in the plane-wave Fredkin-Koehler method has only four entries and is given by First, we only calculate the frequencies and spatial profiles of the spin-wave modes at k = 0 under an applied static external field parallel to the layer [see inset in Fig. 4(a)].At this wave number, k = 0, the magnetic precession is homogeneous within the layer plane but can still be inhomogeneous along the layer thickness, forming standing waves along the layer thickness.These modes are typically referred to as ferromagnetic-resonance (FMR) modes or perpendicular-standing spin waves (PSSWs) in common microwave-absorption experiments and can be denoted by an index n counting the number of nodal lines along the thickness [see Fig.
with h ext = H ext /M s being the unitless static external field and κ n = nπ/d being the wave number of the perpendicularstanding waves along the thickness (y) direction of the monolayer.
In Fig. 4(a), we show the oscillation frequencies of the different PSSWs n = 0, 1, 2, 3 in a permalloy monolayer of d = 150 nm thickness, in a field range between 0 and 60 mT, calculated with our dynamic-matrix approach implemented in TETRAX (solid lines), showing a perfect agreement with the theoretical prediction according to Eq. (26) (dashed lines).Note, that the obtained mode profiles along the thickness, shown in Fig. 4(b), are perfect unpinned sinusoidals.For the highest-order mode, n = 3, which exhibits the shortest wavelength along the thickness, a slight frequency-mismatch can be observed which originates from an underestimation of the exchange interaction, i.e., from insufficient accuracy when calculating the magnetization derivatives along the thickness of the layer (y direction).This, of course, could be improved simply by decreasing the characteristic length of the mesh.Here, the layer has been modeled on a line mesh with an average spacing of 1 nm between the nodes.

B. Spin-wave dispersion in thick films
Using the same geometry as in the previous section, we now calculate the dispersion of propagating spin waves (k = 0), starting with a thin layer of d = 10 nm thickness.Throughout this section, the monolayer is saturated in-plane by a constant external field of 20 mT.As the wave number k departs from zero, the frequency of the spin waves depends crucially on the orientation of their wave vector with respect to the equilibrium magnetization -a symmetry breaking which is introduced by the dipolar interaction.This can be seen for the two limiting cases of k m 0 and k ⊥ m 0 in Fig. 5(a).Commonly, the spin waves with k m 0 , which propagate parallel to the equilibrium magnetization, are referred to as backwardvolume magnetostatic waves (BVMSWs) due to the fact that, with increasing k, they are localized mainly to the volume of the layer and can, with increasing layer thickness d, exhibit a negative group velocity in certain regions of the wave-vector space.In contrast, the spin waves with k ⊥ m 0 , propagating perpendicular to the equilibrium, generally exhibit a much higher group velocity and, depending on the propagation direction, are localized to either surface of the layer.Hence, they are also referred to as magnetostatic surface waves (MSSWs).Solid and dashed lines correspond to results of TetraX and predictions of the zeroth-order perturbation theory of Kalinikos and Slavin 45 (Eq.21), respectively.(b) Sketch of the film with the used coordinates, the external field direction as well as the two main propagation direction, namely the backward-volume (BV) and surface waves (SW) geometry.Two exemplary mode profiles along the film thickness of the main directions are shown in (c).In (d) and (e), the dispersion relation of a 50 nm and 75 nm thin film is calculated and compared with the analytical predictions (dashed lines).The inset in (d) highlights the dipole-dipole mode hybridization between the first two branches in the SW geometry.In panel (f) the comparison of the dispersion, computed with three different numerical codes, for a 100 nm film shows a perfect agreement between the different numerical codes.The inset shows that the finite difference code, SWIIM, and our finite element code, TetraX, perfectly overlap even for the computationally most critical branch hybridization region.
From a theoretical point of view, the spin-wave propagation in thin films is most prominently described by the perturbation theory of Kalinikos and Slavin 45 (KS), in which dipolar fields are calculated in terms of magnetostatic Green's functions.For sufficiently thin layers, neglecting surface pinning and hybridization between different modes, the zeroth-order perturbation of KS provides explicit analytical expressions for the dispersion of the different modes n.For the two limiting cases (BVMSW and MSSW), the dispersion can be written as for k m 0 (BVMSW) and for k ⊥ m 0 (MSSW), with k 2 n = k 2 + κ 2 n being the square of the total wave vector and P nn being given by |k|d .
(27c) For a layer thickness of d = 10 nm, we see in Fig. 5(a) that our numerical calculations are in perfect agreement with the theoretical prediction made by the zeroth-order perturbation of KS in Eq. (27).With increasing layer thickness, the higherorder modes (n > 0) decrease in overall frequency due to decreasing confinement and, therefore, an overall decrease in exchange energy.At the same time, the dispersion curve of the zeroth MSSW mode (n = 0) acquires a much steeper slope (higher group velocity), as seen for a layer of thickness d = 10 nm in [Fig.5(d)].As soon as two modes cross they can, depending on their symmetry, share an avoided level crossing due to dipole-dipole hybridization, an effect, mostly present for the MSSW modes (red lines).For the layer of d = 50 nm thickness [Fig.5(d)], we see how the zeroth (n = 0) and the second (n = 1) MSSW modes are hybridized, as reflected by an avoided level crossing between them.Naturally, the zeroth-order theory of KS does not capture this feature.Apart from that, already at the thickness of 50 nm, a considerable deviation between our numerical calculations and the analytical theory can be observed, a trend which increases even further with thickness [Figs.5(d-f)].This, however, is to no surprise, as the zeroth-order KS theory does not consider the perturbation of the spatial mode profiles due to dipolar fields yet.Therefore, it leads erroneous results in thick layers, especially for the MSSW (k ⊥ m 0 ) modes as their mode profiles are strongly perturbed by internal dipolar fields.As the zeroth-order theory of KS is still widely used in many works even for thicker samples, it is worth noting that, strictly speaking, this theory is only applicable for very thin layers with thicknesses below the order of a couple of exchange lengths of the respective material (usually below 20 nm).For larger thicknesses, higher-order terms in the perturbation series need to be included, which is already welldescribed in the seminal work of KS.Recall that a dynamicmatrix approach, such as the one presented here, is a method relying on direct numerical diagonalization of the linearized equation of motion.Therefore, by design, it always provides the exact normal modes of the respective magnetic system (up to discretization errors) and, in principle, is applicable for all thicknesses.
In order to verify the correctness of our finite-element dynamic-matrix approach for thicknesses where Eq. ( 27) is not valid anymore, we compare our results to calculations performed on the same system using two different finitedifference (FD) codes: SWIIM, developed by Henry et al., 24 as well as a FD approach by Gallardo et al. 30 Both of these methods model the layer on a regular one-dimensional chain of constant spacing along the thickness.Finally, in Fig. 5(f), we show the dispersion of the lowest four MSSW modes (k ⊥ m 0 ) calculated for d = 100 nm using TETRAX as well as the two FD dynamic-matrix approaches, showing a perfect agreement.For visual clarity, the BVMSWs have been omitted for which the zeroth-order KS theory already provided a good approximation even for larger thicknesses.For completeness, a comparison between TETRAX and the two FD codes for all thicknesses between 10 and 75 nm is found in the supplementary material, also showing no discrepancies.
As another suitable display of how accurate dipolar fields are calculated in our method, as an inset in Fig. 5(f), we show a zoom-in on the extremely narrow avoided level crossing which appears due to dipole-dipole hybridization between the zeroth and the second MSSW mode.Here, we also obtain a perfect agreement with the finite-difference calculations.In conclusion, we have verified the correct calculation of dipolar and exchange fields in our numerical scheme of propagating spin waves in a monolayer.

C. Asymmetric spin-wave dispersion in antiferromagnetically-coupled bilayers
After we have validated the correctness of our method for a single magnetic monolayer, as a final example, we want to extend our consideration to a magnetic bilayer system.In particular, we consider two layers of the same thickness d, separated by a non-magnetic spacer of thickness s, as depicted in Fig. 6(a).The material parameters are the same as in the previous sections.For such a symmetric layer stack, the 4 × 4 Dirichlet matrix is given as When two magnetic layers are brought into proximity, their spin-wave spectra hybridize (via dynamic dipolar fields or possible dynamic interlayer-exchange fields).This leads, for example, to the mixing of the lowest modes of each layer either into an in-phase (acoustic) or an out-of-phase (optical) mode shown in Fig. 6(b).At k = 0, that is, for magnetic oscillations homogeneous within the bilayer plane, the optical and acoustic modes are degenerate as long as the layers are not coupled via interlayer exchange, J bl = 0.This degeneracy at k = 0 is lifted by a non-zero interlayer-exchange coupling, J bl = 0, as seen in Fig. 6(c) In case of two layers magnetized antiparallel to each other, spin waves propagating perpendicular to the two magnetizations are non-reciprocal [see Fig. 6(a)].Such a state can be stabilized by an interlayer-exchange coupling with negative sign, J bl < 0, (antiferromagnetic coupling) which, in our case, we set to As a consequence of the antiparallel layer alignment, counter-propagating waves with the same wavelength exhibit different frequencies, seen in Fig. 6(c).This nonreciprocity is purely of dipolar origin and a consequence of magnetochiral symmetry breaking in the pseudo charges generated by the dynamic magnetization.Next to antiparallel alignment of the magnetic layers, this dipolar symmetry breaking can also be introduced by surface curvature, 34,46,47 as observed in magnetic nanotubes, or by a chiral magnetic texture, as observed for the spin waves propagating along Bloch walls 48 in systems with perpendicular magnetic anisotropy.For the layer stack considered here, the spin-wave spectrum was studied theoretically and experimentally by Gallardo et al. in Refs.37,49.Finally, in Fig. 6(c), we compare the theoretical dispersion of the acoustic and optical mode in the considered system according to Gallardo et al. with the numerical calculations using our FEM dynamic-matrix approach implemented in TETRAX, for a layer stack with layer thickness d = 2 nm and spacing s = 2 nm.It is possible to see that our numerical scheme is in perfect agreement with the analytical theory of Ref. 37 in the case of thin layers.However, analogous to the previous section, as we increase the layer thickness d to 20 nm, in Fig. 6(d) we can see clear deviations between the theory and our numerical calculations.This again is due to the fact that the theory does not consider dipolar perturbations of the mode profiles, which, for large d, become inhomogeneous along the layer thickness.Because of this, in addition, we compare our results again with the FD difference dynamicmatrix approach by Gallardo et al., 49 which is also capable of modeling inhomogeneities along the layer thickness.As can be seen in Fig. 6(d), the correctness of our calculations with TETRAX is perfectly supported by the FD calculations.
Let us highlight here a major benefit of the plane-wave Fredkin-Koehler method that we use to calculate the dynamic The spin-wave asymmetry versus the propagation vector for a variety of interlayer spacers is summarized in panel (e).Finally, in (f), the maximum asymmetry as a function of the layer spacing is shown for layers with 2 nm thickness.For all panels, the dashed lines are the analytical predictions.
dipolar fields: As this method avoids for the non-magnetic material (here: the spacer/interlayer) to be part of the mesh, the computational effort is completely independent of the exact value of the spacer thickness s which solely appears as a parameter in the Dirichlet matrix Eq. (28).Not only is it easy to continuously vary the spacing s.It is also possible to take the limit of extremely large or small s without ever increasing the number of nodes in the mesh and, therefore, without increasing the computational effort.We note that the same is true when explicitly calculating the plane-wave dipolar tensors N(dip) k in matrix form, as done, for example, in the FD code SWIIM. 24In this case, too, only the magnetic material needs to be modeled (for the convolution with the Green's function), while any inter-cell spacings appear as parameters in the dipolar tensors.Recall, however, that this makes the dynamic-matrix and, therefore, the eigenvalue problem, dense.
To illustrate the flexibility with respect to changing the interlayer spacings in our FEM dynamic-matrix approach, we present how the dipole-induced dispersion asymmetry in the considered layer stack changes when varying the spacing s, for a fixed d = 2 nm.Technically, the sign and magnitude of the interlayer-exchange coupling will of course vary with the thickness s of the interlayer.Therefore, we will completely disregard it here, J bl = 0.This is reasonable in the sense that, for thin layers, interlayer-exchange coupling has no direct quantitative influence on the dispersion asymmetry.
In Fig. 6(e), we show the asymmetry |∆f | = |f (k) − f (−k)|, which, for the case of thin layers, has the same magnitude for the acoustic and optical mode.We start from s = 2 nm, which corresponds to the same case as shown in Fig. 5(c), and go up to s = 10 nm.All solid curves have been calculated using TETRAX in approximately the same amount of time (less than a minute).For visual clarity, the analytical results according to Ref. 37 are only shown for selected spacings.It is possible to see that, with increasing spacing s, the position of the maximum dispersion asymmetry shifts to lower wave numbers |k| while its maximum value decreases.As a figure of merit, in Fig. 6(f) we show the smooth transition of this maximum as a function of spacing s, again, showing perfect agreement between numerics and analytics.

IV. CONCLUSIONS
In summary, we have extended our finite-element dynamicmatrix approach for propagating spin waves for waveguides of arbitrary cross section to mono-and multilayers of arbitrary spacing and thickness.Therefore the dispersion relation for extended films can be calculated by using an 1D linetrace mesh only along the thickness of the ferromagnetic film.To do so, the previously presented Fredkin-Koehler method (also known as the hybrid finite-element / boundary element method) to solve the screened Poisson equation of propagating spin waves was extended for mono-and multilayers, allowing to compute the dipolar potential and related stiffness field in a very efficient manner.Remarkably, the obtained boundary matrix (or Dirichlet matrix) has exact elements, defined by analytical expressions.The major benefit of this method is that it avoids the discretization of any non-magnetic material while also keeping the resulting eigenvalue problem sparse.In particular, this means that the computational ef-fort is completely independent on the spacer thickness, which solely appears as a parameter in the Dirichlet matrix.Moreover, the resulting matrices only scale with the number of nodes n along the normal direction of the layers as O(n), whereas it scales as O(n 2 ) when explicitly calculating the dipolar tensors.This, provides our FEM approach with a comparably low arithmetic complexity.
Our method has been validated for a number of known systems using theoretical predictions as well as finitedifference implementations established previously.The presented method is readily implemented into the TETRAX 38 open-source micromagnetic modeling package offering for the magnonic community an easy and efficient calculation of spin-wave dispersions for various standard magnonic problems.

Figure 1 .
Figure 1.(a) Schematics of an infinitely extended waveguide with arbitrary cross section A. The magnified inset shows the angle attributed to a certain boundary node.(b) The real parts of the different lateral potentials in the plane-wave Fredkin-Koehler method are shown crossing the boundary ∂A of the magnetic sample (Figure adapted from Ref. 1).

Figure 2 .
Figure 2. (a) Schematics of a multilayer waveguide with finite width W composed of magnetic layers separated by non-magnetic spacers, both with arbitrary film thicknesses.Extending the lateral dimension to infinity will allow to model the layers using a one-dimensional line-trace mesh along the thickness (y direction) of each individual ferromagnetic layer only, as shown in (b).(c) The boundary integral in Eq. 12 can be split along the circumference of each cross section into integrals for the different edges with a fixed integration direction.The normal vector n is defined pointing outwards.

Figure 3 .
Figure 3. Schematics of two ferromagnetic layers separated by a nonmagnetic interlayer and exchange coupled due to RKKY interaction.Si and Sj are the surfaces of the bilayer system adjacent to the interlayer, while ci and cj are the "volumes" (being the length in the case of the 1D line-trace mesh) of the Wigner-Seitz cells associated to the boundary nodes.

Figure 4 .
Figure 4. (a) External-field dependence of the uniform mode and the perpendicular-standing waves along the thickness (n being the node number) in a 150 nm thick permalloy film, with the field applied in the plane of the film.A schematics of the magnetic film, the definition of the coordinate system and the equilibrium magnetic state is represented by the inset.In panel (b), the out-of-plane component of the magnetization (my) along the thickness for the different modes is shown.

FrequencyFigure 5 .
Figure 5. (a) Dispersion relation of the first spin-wave branch in a 10 nm thick film shown for the two main propagation directions.Solid and dashed lines correspond to results of TetraX and predictions of the zeroth-order perturbation theory of Kalinikos and Slavin45 (Eq.21), respectively.(b) Sketch of the film with the used coordinates, the external field direction as well as the two main propagation direction, namely the backward-volume (BV) and surface waves (SW) geometry.Two exemplary mode profiles along the film thickness of the main directions are shown in (c).In (d) and (e), the dispersion relation of a 50 nm and 75 nm thin film is calculated and compared with the analytical predictions (dashed lines).The inset in (d) highlights the dipole-dipole mode hybridization between the first two branches in the SW geometry.In panel (f) the comparison of the dispersion, computed with three different numerical codes, for a 100 nm film shows a perfect agreement between the different numerical codes.The inset shows that the finite difference code, SWIIM, and our finite element code, TetraX, perfectly overlap even for the computationally most critical branch hybridization region.

JFigure 6 .
Figure 6.(a) Schematics of a symmetric bilayer with thickness d and spacing s magnetized antiparallel to each other.The dispersion is calculated for spin waves propagating perpendicular to the static magnetization.The lowest two spin-wave branches are the homogeneous symmetric and antisymmetric modes, with mode profiles sketched in (b).In (c) and (d) the simulated dispersion relation of the modes in (b) in comparison with the analytical solution is shown for bilayers with 2 nm as well as 20 nm thickness and 2 nm separation.As expected and in accordance with previous micromagnetic simulations, for the larger film thickness the analytical results deviate from the simulation results using TETRAX.Instead, in (d), we compare our results also to the finite-difference dynamic-matrix results (crosses) of Gallardo et al. according to Ref.49.The spin-wave asymmetry versus the propagation vector for a variety of interlayer spacers is summarized in panel (e).Finally, in (f), the maximum asymmetry as a function of the layer spacing is shown for layers with 2 nm thickness.For all panels, the dashed lines are the analytical predictions.

Table I .
Parameters used for micromagnetic modeling.