Solving Maxwell's Eigenvalue Problem via Isogeometric Boundary Elements and a Contour Integral Method

We solve Maxwell's eigenvalue problem via isogeometric boundary elements and a contour integral method. We discuss the analytic properties of the discretisation, outline the implementation, and showcase numerical examples.


Motivation
The development and construction of particle accelerators is arguably one of the most time and money consuming research projects in modern experimental physics. A cause for this is that essential components are not available off the shelf and must be manufactured uniquely tailored to the design specification of the planned accelerator.
One of the most performance critical components are so-called cavities, resonators often made out of superconducting materials in which electromagnetic fields oscillate at radio frequencies. The resonant fields are then used to accelerate bunches of particles up to speeds close to the speed of light. The geometry, and consequently the resonance behaviour of these structures is vital to the overall performance of the accelerator as a whole.
Due to the expensive (e.g. superconducting) materials and vast amounts of manual labor that are needed in the manufacturing of these devices, the design of cavities has become its own area of research, cf. [1] and the sources cited therein. Consequently, the development of simulation tools specifically for this purpose became an important part of the related scientific advancement, see e.g. [2,3,4].
A bottleneck with these classical approaches has always been the representation of the geometry, which often limits the achievable accuracies, cf.  : For a finite element computation a mesh is generated from the boundary data available from the design framework. Afterwards, a volume mesh is created. This introduces geometrical errors and limits the obtainable order of convergence to that of the geometry representation. Graphics from [5].
However, high accuracies are desired such that the initial design and its simulation are not the weakest link within the manufacturing pipeline. As an example, manufacturers alone are interested in the simulation of deformation effects, the so-called Lorentz detuning, which are dependent on a relative error margin of roughly 10 −7 , [1, Tab. II]. There are also other, more advanced applications which have such high demands on accuracy. One is presented by Georg et al. [6], who show that even higher accuracies than those already achievable are required to simulate eccentricities.

Introduction
Nowadays, isogeometric analysis [7] has been established as the method of choice when dealing with such high demands on accuracy w.r.t. geometry representation. Isogeometric methods are well understood for the case of electromagnetism [8,9] and a corresponding finite element approach has already been applied to Maxwell's eigenvalue problem [10].
The boundary-only representations in modern CAD frameworks, as well as the demand for highly accurate simulation techniques, suggests the transition to isogeometric boundary element methods for these types of problems, cf. Figure 2. While the use of boundary element methods promises to both reduce the number of degrees of freedom w.r.t. the element size h drastically, and at the same time double the rate of convergence for the point evaluation of the solution in the domain [11,Cor. 3.11 & Rem. 3.12], the system matrices become densely pouplated and the corresponding eigenvalue problem becomes non-linear. However, through the recent introduction of a new family of eigenvalue solvers known as contour integral methods [12,13], this problem was mitigated.
A first approach to a boundary element eigenvalue problem via the contour integral method was investigated in [14], however, neither with a higher-order approach, a discussion of the related convergence theory, or within the isogeo- metric setting. Moreover, a comparison of volume-based and boundary element methods has not been discussed within the literature.
This article is build on top of the recent mathematical results of [11,12,15,16] and aims at advancing the solution of eigenvalue problems via boundary elements by discussing the convergence analysis of the isogeometric discretisation of the eigenvalue problem. In this, we prove that a convergence order O(h 2p+1 ) for the eigenvalue aproximation can be achieved for discretisations with mesh-size h and ansatz functions of order p. For a corresponding volume-based IGA only a convergence order O(h 2p ) can be expected [17].
The organisation of this document is straight forward. Section 3 will introduce the cavity problem based on Maxwell's equations. Afterwards, Section 4 will show how this problem can be rephrased as a boundary integral equation, the well-known electric field integral equation, and we will discuss the equivalence of both formulations. We will then discuss our discretisation scheme, where we will review our isogeometric discretisation of the boundary integral equation in Section 5, making theoretical predictions regarding the convergence behaviour. Our approach to the arising non-linear eigenvalue problem will then be discussed in Section 6. In Section 7 a selection of numerical examples will be presented, and finally we will briefly conclude our findings in Section 8.

The Eigenvalue Problem
For a compact and simply connected Lipschitz domain Ω ∈ R 3 the electromagnetic fields in the source-free case are governed by the equations assuming the time-harmonic case. As usual, E and H denote the electric and magnetic field [18], respectively. In the case of the cavity problem, the electric permittivity ε 0 and the magnetic permeability µ 0 are those of vacuum. Moreover, since superconducting alloys can be modelled as perfect electric conducting, we assume the boundary conditions on Γ := ∂Ω to be given by By eliminating H from (1) one can then derive the classical cavity problem [19].

The Electric Field Integral Equation
Before we introduce the variational formulation of Maxwell's eigenvalue problem, we need to introduce some notation.

Function Spaces Related to Maxwell's Equations and the EFIE
Let Ω be a Lipschitz domain. By H 0 (Ω) we denote the usual square integral functions L 2 (Ω), and we write H s (Ω) for the usual Sobolev spaces of higher regularity s > 0, cf. [20].
dropping the index s in the case of s = 0 for convenience. For n x x x0 denoting the outwards directed unit normal at x x x 0 ∈ Γ, the rotated tangential trace is well defined for smooth f f f . Note that n n n x x x0 is well-defined for almost all x x x 0 , cf. [20, p. 96]. We will equip the rotated tangential trace with the superscript int if the limit is taken from within Ω, and with the superscipt ext if taken from R 3 \ Ω, omitting the notation if the mapping properties are clearly stated. The operator γ γ γ t is extended to a weak setting through density arguments. We can now define H H H (div Γ , Γ) surjective, and one can prove that it is continuous, compare [21].
We can now reformulate Maxwell's eigenvalue problem as a variational problem exclusively on the boundary Γ by using the anti-symmetric pairing cf. [21, Def. 1].

Recasting the Eigenvalue Problem
Any solution of the electric wave equation can be derived via a boundary integral formulation, which we will review within this section. The version reviewed here resembles the one presented in [21,Thm. 6].
We define the Maxwell single layer potential The mapping properties are known, see [21,Thm. 5], where it is also shown that the image of V V V κ is divergence free. With the help of this operator, one can show the following.
Lemma 2 (Single Layer Representation, [21,Thm. 6]). If E ∈ H(curl curl curl 2 , Ω) is solution to the electric wave equation (3) on Ω, then it can be represented via We now define the Maxwell single layer operator Applying the rotated tangential trace to the identity in Lemma 2 allows us to recast the eigenvalue problem (3) as a variational problem w.r.t. the duality pairing (4). The underlying boundary integral equation reads as follows: The problem in (6) is a nonlinear eigenvalue problem with respect to the eigenvalue parameter κ, since κ occurs nonlinearly in the fundamental solution (5) which builds the kernel of the single layer boundary integral operator V V V κ . The eigenvalue problem formulations (3) and (6) are equivalent in the following sense: Lemma 4 (Equivalence of eigenvalue problems). Let κ ∈ R and κ > 0.
Proof. Assertion a) has been already shown by the derivation of the boundary integral formulation (6). We want to mention that the eigenvalue problem (6) has in addition to the real eigenvalues also non-real eigenvalues which correspond to the eigenvalues of the exterior eigenvalue problem. We refer to [16] for an analysis of this kind of eigenvalues.

A Brief Review of Isogeometric Analysis
This section is devoted to a brief review of isogeometric analysis as required for its utilisation in the context of boundary element methods for electromagnetism. We follow the lines of [22,11], which in their turn are based on the works of Buffa et al. [23] aimed at finite element discretisations.
Let p ≥ 0 and m > 0 be integers. While all of the results reviewed in this article are applicable to so-called p-open locally quasi-uniform knot vectors, cf. [24, Assum. 2.1], we assume all knot vectors to be of the form This is introduced only for notational convenience and in agreement with our numerical examples. We will then refer to m as the refinement level and define the mesh size h := 2 −m . We proceed to define B-spline bases via the well-known recursion formula We define the spline space S p m as the span of the B-spline basis defined on Ξ p m . In reference domain, i.e., on := (0, 1) 2 , we can now define the divergence conforming spline space S S S p m ( ) as done in [24, Sec. 5.5] via for ⊗ denoting the tensor product and × denoting the Cartesian product.

Geometry and Discretisation in the Physical Domain
Let a patch Γ be given by the image of under an invertible diffeomor- For Ω being a Lipschitz domain, define a multipatch geometry to be a compact, orientable two-dimensional manifold Γ = ∂Ω invoked via 0≤j<N Γ j by a family of patches {Γ j } 0≤j<N , N ∈ N. The family of diffeomorphisms {F F F j : → Γ j } 0≤j<N will be called parametrisation.
We assume the Γ j to be disjoint and that for any patch interface D of the form D = ∂Γ j0 ∩ ∂Γ j1 = ∅ we require the continuous extensions of the parametrisations F F F j0 and F F F j1 onto = [0, 1] 2 to coincide.
Within the framework of isogeometric analysis parametrisations will usually be given by tensor products of non-uniform rational B-splines (NURBS), i.e., functions of the form .
To extenf the definition of S S S p m to the physical domain, we resort to an application of the so-called Piola transformation [26,Chap. 6].
For the case of geometry mappings F F F j : (div Γ , Γ) that allow for point evaluations, its explicit form is given by Therein, η denotes the surface measure and d d dF F F j denotes the Jacobian of F F F j . Although it is not readily invertible, an inverse exists in the sense of a mapping from the two-dimensional tangent space of Γ j to vector fields on , as has already been discussed in [11,Sec. 3.2]. We remark that it needs not be computed, since all computations can be reduced to computations in the reference domain. Explicit formulae are easily derived, see e.g. [26,Sec. 6.3].
We can now introduce the globally divergence conforming space It has been analysed in [15,Def. 10], where it has been shown that it enjoys quasioptimal approximation properties w.r.t. energy norm of the EFIE. Specifically, the spline space satisfies estimates of the kind for We remind of the relation h = 2 −m due to the assumption of uniform refinement on the knot vectors. These discrete spaces give rise to the discrete analogue of Problem 3.
Problem 5 (Discrete Eigenvalue Problem). Find a non-zero j j j h ∈ S S S p m (Γ) and κ h such that holds for all µ µ µ ∈ S S S p m (Γ).
In its discrete form, this problem induces a linear system which we choose to assemble as in [11,27]. Therein, V V V κ can be interpreted as a matrix valued function dependent on κ.
The well-posedness of this discrete problem is closely related to the following criteria, cf. [28,Sec. 3].

(C1) There exists a continuous splitting H H H
hold.

Numerical Analysis
The convergence of the eigenvalues and eigenfunctions of the Galerkin eigenvalue problem (8) can be shown using abstract results of [29,16] and [30,31]. Crucial for the convergence analysis are the above listed criteria (C1)-(C4), that V V V κ satisfies a T -Gårding's inequality with respect to the splitting in (C1) [11], and that the mapping C \ {0} κ → V V V κ is holomorphic [16,Lem. 2]. In [29,16] sufficient conditions for the convergence of conforming Galerkin methods are specified for eigenvalue problems for holomorphic T -Gårding operator-valued functions. According to [16,Lem. 5] and [29, Lem. 2.7] the criteria (C1)-(C4) allow the application of the convergence theory established in [30,31] to the Galerkin eigenvalue problem (8). From the comprehensive convergence results presented in [30,31,29,16] we only want to state the error estimate for semisimple eigenvalues κ: where δ m,p := sup j j j∈kerV V Vκ j j j The quadratic convergence order with respect to δ h,p follows from the fact that for any eigenfunction j j j adj of the adjoint eigenvalue problem there exists an eigenfunction j j j of the eigenvalue problem V V V κ j j j = 0 such that j j j = j j j adj , see [16,Lem. 3]. The estimate (7) together with (11) and (12) yields the final estimate for sufficiently smooth surface current densities j j j.
Remark 6 (Convergence on non-smooth geometries). In this sense sufficiently smooth needs to be understood in terms of patchwise regularity as explained in [15]. In general, for surface current densities on non-smooth geometries this smoothness assumption will not be fulfilled, since j j j may admit singularities at corners and edges. However, for the specific densities of resonant modes within interior cavities this smoothness assumption will often be fulfilled. For the cube an analytical representation of j j j is known to be smooth [14], and even for other, non-trivial geometries this can be argued.

The Contour Integral Method
This chapter will give a short summary of the contour integral method as introduced by Beyn [12], without any non-trivial modifications. Note that there are alternative approaches by other authors, where the first publication [13] seems to go back to 2009, compare also [32]. After this short review, we will discuss three numerical examples that utilise this method to solve Problem 8.
Let T : M → C m×m be holomorphic on some domain M ⊂ C. Let T h denote the Hermitian transpose of T and T denote its usual transpose. We want to solve the nonlinear eigenvalue problem If z is no eigenvalue of T, we find by the condition above that the kernel of T(z) is trivial, thus T(z) has full rank and is invertible.
Fix some domain D ⊆ M with boundary ∂D, and assume there to be k eigenvalues (λ j ) j≤k in the interior of D, such that all eigenvalues are simple. We remark that the entire approach was also generalised to that case of non simple eigenvalues [12]. As a consequence of the version of the Keldysh Theorem stated in [33] one can proof the following result, on which this method is based.
Theorem 7 (Contour Integral Theorem, [12,Thm. 2.9.]). For T as above and holomorphic f : M → C it holds that where v j and w j are left and right eigenvectors corresponding to λ n that are normalized according to For our purposes, where we set T(z) = V z as a matrix-valued function, cf. (9), the contour integral method can be reduced to the following theorem. The proof of this theorem is constructive and corresponds to the following algorithm, on which our implementation is based.
At first the algortihm might seem prohibitively expensive due to the application of a singular value decomposition. However, since the number of eigenvalues k will be small and often a reasonable upper estimate 0 of the number of eigenvalues is known such that one can choose = 0 . Then the complexity of the singular value decomposition becomes negligible in comparison to solving the linear system in lines 5 and 15 of Algorithm 9.
Moreover, in an actual implementation A 1 and A 0 are assembled simultaneously, since the most expensive operation of the algorithm is evaluating T(t j ) −1V . In general, k and will be small, so storing both A j poses no issues. Since smooth contours should be used, exclusively, one can expect the trapezoidal rules for the assembly of the A j to converge exponentially with respect to N . Thus, the bottleneck in terms of accuracy of the entire scheme is the accuracy in which X represents the bilinear form of the EFIE.
Remark 10 (Obtaining the Eigenvectors). Note that through this algorithm we not only obtain the eigenvalues of Problem 3, but also the coefficients of the corresponding eigenfunctions in the form of the matrices V and W. Note moreover that the number of non-zero singular values reflects the number of solutions within D and can be used for verification of an implementation if the number of solutions is known from analytical representations or measurements.

Numerical Examples
In the following we will discuss some numerical experiments showcasing the application of the contour integral method to the isogeometric boundary element method. The particular implementation used is Bembel, which is available opensource [34,35]. A branch containing the code to recreate all of the presented numerical examples is also available [36].
Since no quasi-optimal preconditioners for the isogeometric discretisation of the electric field integral equation are known, iterative solvers yield suboptimal runtimes. Thus, the following numerical experiments will use a dense matrix assembly together with the partially pivoted LU decomposition of the linear algebra library Eigen [37] to solve the arising systems. The utilised higher-order approaches yield systems small enough for this approach to be more than feasible.
As a first example, we investigate the first two eigenvalues of the unit cube, where an analytical solution is given by the eigenvalues λ ana,0 = π √ 2 of multiplicity three and λ ana,1 = π √ 3 of multiplicity two. The ellipse was defined as where the contour integrals were evaluated, again, with N = 25. The error visualised corresponds to error = 1 5 0≤j<5 min i=0,1 |λ j − λ ana,i |.  As one can see in Figure 3, the multiplicity of the eigenvalues is reflected perfectly by the non-zero singular values, i.e., all eigenvalues have been recognised. Moreover, we have the theoretically obtainable convergence order of O(h 2p+1 ).

Comparison to IGA-FEM and Commercial Tools
As a second example, we compute the first eigenvalue for the sphere, cf. Figure 4 for the results. The contour was chosen as the curve 0.5 · sin(t) + i · 0.05 · cos(t) + 2.5, for t = [0, 2π). The convergence behaviour w.r.t. h matches the orders O(h 2p+1 ) predicted by (13) once more. We have also computed approximations of the smallest eigenvalues of the unit sphere with the volume-based IGA software package GeoPDEs 3.0 [38,39]. In Table 1 we showcase results for different polynomial degrees and refinements the number of degrees of freedom and the reached accuracy for the IGA-BEM and the volume-based IGA. One can observe that for the IGA-BEM a notably fewer (at least 20×) number of degrees of freedom are necessary for each polynomial degree in order to reach the same accuracy as for the volume-based IGA. The difference to commercial tools, in this case CST Microwave Studio 2018, is even more pronounced, cf. Table 2. However, the matrices for the IGA-BEM are dense and the eigenvalue problem non-linear.

An Industrial Application: TESLA Cavity
As a third example we discuss an eigenvalue computation of the one-cell TESLA geometry as shown in Figures 1 and 2, with the results depicted in Figure 6. For this example no analytical solution is known, but experience dictates a resonant frequency around κ ≈ 26.5. We choose the contour 0.5 · sin(t) + i · 0.01 · cos(t) + 26.5, for t = [0, 2π) and N = 8. As a reference solution we utilise the result of a computation with p = 5 and h = 1/8. This reference solution was compared to a computation   with CST Microwave Studio 2018. The set solver order was 3rd (constant) and the mesh was generated with 200 771 curved tetrahedral elements of order 5. CST yields the solution of 1.27666401260 GHz. Our reference solution of λ cim = 26.75690023 corresponds to a frequency of 1.276664064 GHz. This results in a relative error of 4.018e-08. Thus our experiments are in good agreement with those of the CST software. However, note that in Figure 6 one can clearly see stagnation w.r.t. the CST Solution on higher-order computations and small h. This suggests that the Bembel reference solution is more accurate, provided the convergence of the contour integral approach behaves as observed in the previous numerical experiments. The order of convergence for the cavity is not as pronounced as in the examples with analytical solution, which, in light of the estimate (13) and Remark 6 is most likely due to the reduced regularity of the corresponding eigenfunction. Either way, one can clearly see an increased accuracy in higher-order approaches.

Conclusion
Overall, the contour integral method yields exceptional accuracies in conjunction with our isogeometric boundary element method. Judging from the asymptotic behaviour predicted in (13) and observed in the numerical examples, the accuracy of the combination of isogeometric boundary element method and the contour integral method promises higher orders of convergence to the correct solution compared to currently implemented volume based approaches, since these will not benefit from the convergence order of O(h 2p+1 ). This means that, for the same computational resources, the maximum reachable accuracy of the IGA-BEM with contour integration is higher compared to many volume-based methods. However, due to the long time spend solving the systems, this means that the IGA-BEM with contour integral offers this higher accuracy only in exchange for runtime, until efficient preconditioning strategies are available, thus making the utilisation of fast methods and iterative solvers viable for this type of application.