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

Stefan Kurz<EMAIL_ADDRESS>Sebastian Schöps<EMAIL_ADDRESS>darmstadt.de Gerhard Unger<EMAIL_ADDRESS>Felix Wolf<EMAIL_ADDRESS>darmstadt.de Technische Universität Darmstadt, Centre of Computational
Engineering  Graz University of Technology, Institute of Applied Mathematics

###### Abstract

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.

###### keywords:

IGA , BEM , Eigenvalue Problem , Cavity Problem

###### MSC:

[2010] 34L16, 35P30, 65N38, 65D07

††journal: arXiv

##  1 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].

Figure 1: 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].

A bottleneck with these classical approaches has always been the
representation of the geometry, which often limits the achievable accuracies,
cf. Figure 1. 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.

##  2 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].

Figure 2: Isogeometric boundary element methods enable the computation
directly on the CAD representation. Graphics from [5].

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
isogeometric 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 ${\mathcal{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
${\mathcal{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.

##  3 The Eigenvalue Problem

For a compact and simply connected Lipschitz domain $\Omega\in\mathbb{R}^{3}$
the electromagnetic fields in the source-free case are governed by the
equations

$\displaystyle{\boldsymbol{\operatorname{curl}}}{}{(\mathbf{E})}$
$\displaystyle=-i\omega\mu_{0}\mathbf{H}$ $\displaystyle\text{in }\Omega$ (1)
$\displaystyle{\boldsymbol{\operatorname{curl}}}{}{(\mathbf{H})}$
$\displaystyle=i\omega\varepsilon_{0}\mathbf{E}$ $\displaystyle\text{in
}\Omega$ $\displaystyle\operatorname{div}{}{(\varepsilon_{0}\mathbf{E})}$
$\displaystyle=0$ $\displaystyle\text{in }\Omega$
$\displaystyle\operatorname{div}{}{(\mu_{0}\mathbf{H})}$ $\displaystyle=0$
$\displaystyle\text{in }\Omega,$

assuming the time-harmonic case. As usual, $\mathbf{E}$ and $\mathbf{H}$
denote the electric and magnetic field [18], respectively. In the case of the
cavity problem, the electric permittivity $\varepsilon_{0}$ and the magnetic
permeability $\mu_{0}$ are those of vacuum. Moreover, since superconducting
alloys can be modelled as perfect electric conducting, we assume the boundary
conditions on $\Gamma\coloneqq\partial\Omega$ to be given by

$\displaystyle\mathbf{E}\times\mathbf{n}$ $\displaystyle=0$
$\displaystyle\text{on }\partial\Omega$ (2)
$\displaystyle\mathbf{H}\cdot\mathbf{n}$ $\displaystyle=0$
$\displaystyle\text{on }\partial\Omega.$

By eliminating $\mathbf{H}$ from (1) one can then derive the classical cavity
problem [19].

######  Problem 1 (Cavity problem).

Find the wave number $k:=\omega\sqrt{\mu_{0}\epsilon_{0}}\in\mathbb{R}$ and
$\mathbf{E}\neq 0$ such that

$\displaystyle{\boldsymbol{\operatorname{curl}}}{}{\bigl{(}{\boldsymbol{\operatorname{curl}}}{}{(\mathbf{E})}\bigr{)}}$
$\displaystyle=k^{2}\mathbf{E}$ $\displaystyle\text{in }\Omega$ (3)
$\displaystyle\operatorname{div}{}{(\mathbf{E})}$ $\displaystyle=0$
$\displaystyle\text{in }\Omega$ $\displaystyle\mathbf{E}\times\mathbf{n}$
$\displaystyle=0$ $\displaystyle\text{on }\Gamma.$

##  4 The Electric Field Integral Equation

Before we introduce the variational formulation of Maxwell’s eigenvalue
problem, we need to introduce some notation.

###  4.1 Function Spaces Related to Maxwell’s Equations and the EFIE

Let $\Omega$ be a Lipschitz domain. By $H^{0}(\Omega)$ we denote the usual
square integral functions $L^{2}(\Omega)$, and we write $H^{s}(\Omega)$ for
the usual Sobolev spaces of higher regularity $s>0$, cf. [20]. Their vector-
valued counterparts
$\boldsymbol{H}^{s}(\Omega)=\big{(}H^{s}(\Omega)\big{)}{}^{3}$ are denoted by
bold letters. We now define the spaces

$\displaystyle\boldsymbol{H}^{s}({\boldsymbol{\operatorname{curl}}}{},\Omega){}$
$\displaystyle{}=\\{\boldsymbol{f}\in\boldsymbol{H}^{s}\colon{\boldsymbol{\operatorname{curl}}}{}\boldsymbol{f}\in\boldsymbol{H}^{s}(\Omega)\\},$
$\displaystyle\boldsymbol{H}^{s}({\boldsymbol{\operatorname{curl}}}{}^{2},\Omega){}$
$\displaystyle{}=\\{\boldsymbol{f}\in\boldsymbol{H}^{s}({\boldsymbol{\operatorname{curl}}}{},\Omega)\colon{\boldsymbol{\operatorname{curl}}}{}\,{\boldsymbol{\operatorname{curl}}}{}\boldsymbol{f}\in\boldsymbol{H}^{s}(\Omega)\\},\text{
and}$ $\displaystyle\boldsymbol{H}^{s}(\operatorname{div}{},\Omega){}$
$\displaystyle{}=\\{\boldsymbol{f}\in\boldsymbol{H}^{s}\colon\operatorname{div}{}\boldsymbol{f}\in
H^{s}(\Omega)\\},$

dropping the index $s$ in the case of $s=0$ for convenience. For
$\mathbf{n}_{\boldsymbol{x}_{0}}$ denoting the outwards directed unit normal
at $\boldsymbol{x}_{0}\in\Gamma$, the _rotated tangential trace_

$\boldsymbol{\gamma}_{\mathrm{t}}(\boldsymbol{f})\coloneqq\lim_{\boldsymbol{x}\to\boldsymbol{x}_{0}}\boldsymbol{f}(\boldsymbol{x})\times\mathbf{n}_{\boldsymbol{x}_{0}}$

is well defined for smooth $\boldsymbol{f}.$ Note that
$\boldsymbol{n}_{\boldsymbol{x}_{0}}$ is well-defined for almost all
$\boldsymbol{x}_{0},$ cf. [20, p. 96]. We will equip the rotated tangential
trace with the superscript int if the limit is taken from within $\Omega$, and
with the superscipt ext if taken from
$\mathbb{R}^{3}\setminus\overline{\Omega},$ omitting the notation if the
mapping properties are clearly stated. The operator
$\boldsymbol{\gamma}_{\mathrm{t}}$ is extended to a weak setting through
density arguments. We can now define
$\boldsymbol{H}_{\times}^{-1/2}(\operatorname{div}{}_{\Gamma},\Gamma)=\boldsymbol{\gamma}_{\mathrm{t}}\big{(}\boldsymbol{H}({\boldsymbol{\operatorname{curl}}}{},\Omega)\big{)}.$
By definition this renders the rotated tangential trace
$\boldsymbol{\gamma}_{\mathrm{t}}\colon\boldsymbol{H}({\boldsymbol{\operatorname{curl}}}{},\Omega)\to\boldsymbol{H}_{\times}^{-1/2}(\operatorname{div}{}_{\Gamma},\Gamma)$
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 $\Gamma$ by using the anti-symmetric pairing

$\displaystyle\langle\boldsymbol{\nu},\boldsymbol{\mu}\rangle_{\times}\coloneqq\int_{\Gamma}(\boldsymbol{\nu}\times\mathbf{n})\boldsymbol{\mu}\mathrm{d}\Gamma,$
(4)

cf. [21, Def. 1].

###  4.2 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_

${\boldsymbol{V}_{\kappa}}{}\colon\mathbf{H}^{-1/2}_{\times}(\operatorname{div}{}_{\Gamma},\Gamma)\to\mathbf{H}(\mathbf{{\boldsymbol{\operatorname{curl}}}{}}^{2},\Omega)$

via the Helmholtz fundamental solution

$\displaystyle
u^{*}_{\kappa}(\boldsymbol{x},\boldsymbol{y})=\frac{e^{-i\kappa\|\boldsymbol{x}-\boldsymbol{y}\|}}{4\pi\|\boldsymbol{x}-\boldsymbol{y}\|}$
(5)

as

$\displaystyle{\boldsymbol{V}_{\kappa}}{}(\boldsymbol{\mu})(\boldsymbol{x})=\int_{\Gamma}u^{*}_{\kappa}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\mu}(\boldsymbol{y})\mathrm{d}\Gamma_{\boldsymbol{y}}+\frac{1}{\kappa^{2}}\boldsymbol{\operatorname{grad}{}}_{\boldsymbol{x}}\int_{\Gamma}u^{*}_{\kappa}(\boldsymbol{x},\boldsymbol{y})(\operatorname{div}{}_{\Gamma}\circ\boldsymbol{u})(\boldsymbol{y})\mathrm{d}\Gamma_{\boldsymbol{y}}.$

The mapping properties are known, see [21, Thm. 5], where it is also shown
that the image of ${\boldsymbol{V}_{\kappa}}{}$ is divergence free. With the
help of this operator, one can show the following.

######  Lemma 2 (Single Layer Representation, [21, Thm. 6]).

If
$\mathbf{E}\in\mathbf{H}(\mathbf{{\boldsymbol{\operatorname{curl}}}{}}^{2},\Omega)$
is solution to the electric wave equation (3) on $\Omega$, then it can be
represented via

$\displaystyle\mathbf{E}={\boldsymbol{V}_{\kappa}}{}\big{(}{\mathbf{\gamma}}_{\mathbf{t}}^{\mathrm{int}}\mathbf{(}{\boldsymbol{\operatorname{curl}}}{}\;\mathbf{E})\big{)}.$

We now define the _Maxwell single layer operator_
${\boldsymbol{\mathcal{V}}_{\kappa}}{}=\boldsymbol{\gamma}_{\mathrm{t}}^{\mathrm{int}}\circ{\boldsymbol{V}_{\kappa}}{}$.
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:

######  Problem 3 (Variational Eigenvalue Problem).

Find a non-zero
$\boldsymbol{j}\in\boldsymbol{H}_{\times}^{{-1/2}}(\operatorname{div}{}_{\Gamma},\Gamma)$
and $\kappa>0$ such that

$\displaystyle\langle{\boldsymbol{\mathcal{V}}_{\kappa}}{}\boldsymbol{j},\boldsymbol{\mu}\rangle_{\times}=0$
(6)

holds for all
$\boldsymbol{\mu}\in\boldsymbol{H}_{\times}^{{-1/2}}(\operatorname{div}{}_{\Gamma},\Gamma)$.

The problem in (6) is a nonlinear eigenvalue problem with respect to the
eigenvalue parameter $\kappa$, since $\kappa$ occurs nonlinearly in the
fundamental solution (5) which builds the kernel of the single layer boundary
integral operator ${\boldsymbol{\mathcal{V}}_{\kappa}}{}$. The eigenvalue
problem formulations (3) and (6) are equivalent in the following sense:

######  Lemma 4 (Equivalence of eigenvalue problems).

Let $\kappa\in\mathbb{R}$ and $\kappa>0$.

  1. 1.

If $(\kappa,\mathbf{E})$ is an eigenpair of (3), then
$(\kappa,({\mathbf{\gamma}}_{\mathbf{t}}^{\mathrm{int}}\mathbf{(}{\boldsymbol{\operatorname{curl}}}{}\;\mathbf{E}))$
is an eigenpair of (6).

  2. 2.

If $(\kappa,\boldsymbol{j})$ is an eigenpair of (6), then
$(\kappa,{\boldsymbol{V}_{\kappa}}{}\boldsymbol{j})$ is an eigenpair of (3).

###### Proof.

Assertion a) has been already shown by the derivation of the boundary integral
formulation (6). For assertion b) note that
${\boldsymbol{V}_{\kappa}}{}\boldsymbol{j}$ is a solution of Maxwell’s
equation in $\Omega$ [21, Sect. 4] and that
$0={\boldsymbol{\mathcal{V}}_{\kappa}}{}\boldsymbol{j}=\boldsymbol{\gamma}_{\mathrm{t}}^{\mathrm{int}}{\boldsymbol{V}_{\kappa}}{}\boldsymbol{j}$.
It remains to show that ${\boldsymbol{V}_{\kappa}}{}\boldsymbol{j}\neq 0$ in
$\Omega$ which follows form the unique solvability of the exterior problem and
the jump relations of the single layer potential on the boundary $\Gamma$, see
[16, Prop. 2.1(ii)]. ∎

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.

##  5 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\geq 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

$\Xi^{p}_{m}=\big{[}\underbrace{0,\cdots,0}_{p+1\text{
times}},1/2^{m},\cdots,(2^{m}-1)/2^{m},\underbrace{1,\cdots,1}_{p+1\text{
times}}\big{]}.$

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\coloneqq 2^{-m}.$ We proceed to define B-spline
bases via the well-known recursion formula

$\displaystyle b_{i}^{p}(x)$
$\displaystyle=\frac{x-\xi_{i}}{\xi_{i+p}-\xi_{i}}b_{i}^{p-1}(x)+\frac{\xi_{i+p+1}-x}{\xi_{i+p+1}-\xi_{i+1}}b_{i+1}^{p-1}(x),$

anchored for $p=0$ by the locally constant functions
$b_{i}^{0}(x)=\chi_{[\xi_{i},\xi_{i}+1)},$ cf. [25, Sec. 2.2]. We define the
spline space $S^{p}_{m}$ as the span of the B-spline basis defined on
$\Xi^{p}_{m}.$

In reference domain, i.e., on $\square\coloneqq(0,1)^{2}$, we can now define
the divergence conforming spline space
$\boldsymbol{\mathbb{S}}^{p}_{m}(\square)$ as done in [24, Sec. 5.5] via

$\boldsymbol{\mathbb{S}}^{p}_{m}(\square)\coloneqq(S^{p}_{m}\otimes
S^{p-1}_{m})\times(S^{p-1}_{m}\otimes S^{p}_{m}),$

for $\otimes$ denoting the tensor product and $\times$ denoting the Cartesian
product.

###  5.1 Geometry and Discretisation in the Physical Domain

Let a _patch_ $\Gamma$ be given by the image of $\square$ under an invertible
diffeomorphism $\boldsymbol{F}\colon\square\to\Gamma\subseteq\mathbb{R}^{3}$.
For $\Omega$ being a Lipschitz domain, define a _multipatch geometry_ to be a
compact, orientable two-dimensional manifold $\Gamma=\partial\Omega$ invoked
via $\bigcup_{0\leq j<N}\Gamma_{j}$ by a family of patches
$\\{\Gamma_{j}\\}_{0\leq j<N},$ $N\in\mathbb{N}$. The family of
diffeomorphisms
$\\{\boldsymbol{F}_{j}\colon\square\hookrightarrow\Gamma_{j}\\}_{0\leq j<N}$
will be called _parametrisation_.

We assume the $\Gamma_{j}$ to be disjoint and that for any _patch interface_
$D$ of the form
$D=\partial\Gamma_{j_{0}}\cap\partial\Gamma_{j_{1}}\neq\emptyset$ we require
the continuous extensions of the parametrisations $\boldsymbol{F}_{j_{0}}$ and
$\boldsymbol{F}_{j_{1}}$ onto $\overline{\square}=[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

$\displaystyle\boldsymbol{F}_{j}(x,y)\coloneqq\sum_{0\leq
j_{1}<k_{1}}\sum_{0\leq
j_{2}<k_{2}}\frac{\boldsymbol{c}_{j_{1},j_{2}}b_{j_{1}}^{p_{1}}(x)b_{j_{2}}^{p_{2}}(y)w_{j_{1},j_{2}}}{\sum_{i_{1}=0}^{k_{1}-1}\sum_{i_{2}=0}^{k_{2}-1}b_{i_{1}}^{p_{1}}(x)b_{i_{2}}^{p_{2}}(y)w_{i_{1},i_{2}}}.$

To extenf the definition of $\boldsymbol{\mathbb{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
$\boldsymbol{F}_{j}\colon\square\to\Gamma_{j}$ and
$\boldsymbol{f}\in\boldsymbol{H}^{-1/2}_{\times}(\operatorname{div}{}_{\Gamma},\Gamma)$
that allow for point evaluations, its explicit form is given by

$\displaystyle\iota_{\boldsymbol{F}_{j}}(\boldsymbol{f})(\boldsymbol{x})$
$\displaystyle\coloneqq\eta(\boldsymbol{x})(\boldsymbol{\mathrm{d}}\boldsymbol{F}_{j})^{-1}(\boldsymbol{f}\circ\boldsymbol{F}_{j})\big{)}(\boldsymbol{x}),$
$\displaystyle\boldsymbol{x}\in\square.$

Therein, $\eta$ denotes the _surface measure_

$\eta(\boldsymbol{x})\coloneqq\|{\partial_{x_{1}}\boldsymbol{F}_{j}({\boldsymbol{x}})\times\partial_{x_{2}}\boldsymbol{F}_{j}({\boldsymbol{x}})}\|,$

and $\boldsymbol{\mathrm{d}}\boldsymbol{F}_{j}$ denotes the Jacobian of
$\boldsymbol{F}_{j}.$ Although it is not readily invertible, an inverse exists
in the sense of a mapping from the two-dimensional tangent space of
$\Gamma_{j}$ to vector fields on $\square,$ 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

$\boldsymbol{\mathbb{S}}^{p}_{m}(\Gamma)\coloneqq\left\\{\boldsymbol{f}\in\boldsymbol{H}_{\times}^{{-1/2}}(\operatorname{div}{}_{\Gamma},\Gamma)\colon\iota_{\boldsymbol{F}_{j}}(\boldsymbol{f}|_{\Gamma_{j}})\in\boldsymbol{\mathbb{S}}^{p}_{m}(\square)\text{
for all }0\leq j<N\right\\}.$

It has been analysed in [15, Def. 10], where it has been shown that it enjoys
quasi-optimal approximation properties w.r.t. energy norm of the EFIE.
Specifically, the spline space satisfies estimates of the kind

$\displaystyle\min_{\boldsymbol{f}_{h}\in\boldsymbol{\mathbb{S}}^{p}_{m}(\Gamma)}\|{\boldsymbol{f}-\boldsymbol{f}_{h}}\|_{\boldsymbol{H}_{\times}^{{-1/2}}(\operatorname{div}{}_{\Gamma},\Gamma)}\leq
Ch^{s+1/2}\|\boldsymbol{f}\|_{\boldsymbol{H}^{s}_{\mathrm{pw}}(\operatorname{div}{}_{\Gamma},\Gamma)},\qquad
0\leq s\leq p,$ (7)

for all
$\boldsymbol{f}\in\boldsymbol{H}^{s}_{\mathrm{pw}}(\operatorname{div}{}_{\Gamma},\Gamma)$,
cf. [15, Thm. 3], where the space
${\boldsymbol{H}^{s}_{\mathrm{pw}}(\operatorname{div}{}_{\Gamma},\Gamma)}$ is
equipped with the norm

$\|\boldsymbol{f}\|_{\boldsymbol{H}^{s}_{\mathrm{pw}}(\operatorname{div}{}_{\Gamma},\Gamma)}\coloneqq\sum_{0\leq
j<N}\big{\|}{{\iota_{\boldsymbol{F}_{j}}(\boldsymbol{f}|_{\Gamma_{j}})\big{\|}_{\boldsymbol{H}^{s}(\operatorname{div}{},\square)}}}.$

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 $\boldsymbol{j}_{h}\in\boldsymbol{\mathbb{S}}^{p}_{m}(\Gamma)$
and $\kappa_{h}$ such that

$\displaystyle\langle{\boldsymbol{\mathcal{V}}_{\kappa}}{}\boldsymbol{j}_{h},\boldsymbol{\mu}\rangle_{\times}=0$
(8)

holds for all $\boldsymbol{\mu}\in\boldsymbol{\mathbb{S}}^{p}_{m}(\Gamma)$.

In its discrete form, this problem induces a linear system

$\displaystyle{\boldsymbol{\mathrm{V}}_{\kappa}}{}\mathbf{j}_{h}=0,$ (9)

which we choose to assemble as in [11, 27]. Therein,
${\boldsymbol{\mathrm{V}}_{\kappa}}{}$ can be interpreted as a matrix valued
function dependent on $\kappa.$

The well-posedness of this discrete problem is closely related to the
following criteria, cf. [28, Sec. 3].

  * (C1)

There exists a continuous splitting
$\boldsymbol{H}^{-1/2}_{\times}(\operatorname{div}{}_{\Gamma},\Gamma)=\boldsymbol{W}\oplus\boldsymbol{V}$
such that the bilinear form
$\langle{\boldsymbol{\mathcal{V}}_{\kappa}}{}\cdot,\cdot\rangle_{\times}$ is
stable and coercive on $\boldsymbol{V}\times\boldsymbol{V}$ and
$\boldsymbol{W}\times\boldsymbol{W}$, and compact on
$\boldsymbol{V}\times\boldsymbol{W}$ and $\boldsymbol{W}\times\boldsymbol{V}$.

  * (C2)

$\boldsymbol{\mathbb{S}}^{p}_{m}(\Gamma)$ can be decomposed into a sum
$\boldsymbol{\mathbb{S}}^{p}_{\boldsymbol{m}}(\Gamma)=\boldsymbol{W}_{h}\oplus\boldsymbol{V}_{h}$
of closed subspaces of
$\boldsymbol{H}^{-1/2}_{\times}(\operatorname{div}{}_{\Gamma},\Gamma)$.

  * (C3)

$\boldsymbol{W}_{h}$ and $\boldsymbol{V}_{h}$ are stable under complex
conjugation.

  * (C4)

Both $\boldsymbol{W}_{h}\subseteq\boldsymbol{W}$, as well as the so-called
_gap-property_

$\displaystyle\sup_{\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}}\inf_{\boldsymbol{v}\in\boldsymbol{V}}\frac{\|{\boldsymbol{v}-\boldsymbol{v}_{h}}\|_{\boldsymbol{H}^{-1/2}_{\times}(\operatorname{div}{}_{\Gamma},\Gamma)}}{\|{\mathbf{v}_{h}}\|_{\boldsymbol{H}^{-1/2}_{\times}(\operatorname{div}{}_{\Gamma},\Gamma)}}\stackrel{{\scriptstyle
h\rightarrow 0}}{{\longrightarrow}}0$ (10)

hold.

These properties have been proven for the isogeometric discretisation [11,
Thm. 3.9].

###  5.2 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 ${\boldsymbol{\mathcal{V}}_{\kappa}}{}$ satisfies a
$T$-Gårding’s inequality with respect to the splitting in (C1) [11], and that
the mapping
$\mathbb{C}\setminus\\{0\\}\ni\kappa\mapsto{\boldsymbol{\mathcal{V}}_{\kappa}}{}$
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 semi-simple eigenvalues
$\kappa$:

$|\kappa-\kappa_{h}|=\mathcal{O}(\delta_{m,p}^{2}),$ (11)

where

$\delta_{m,p}:=\sup_{\begin{subarray}{c}\boldsymbol{j}\in\mathrm{ker}{\boldsymbol{\mathcal{V}}_{\kappa}}{}\\\
\|\boldsymbol{j}\|_{\boldsymbol{H}^{-1/2}_{\times}(\operatorname{div}{}_{\Gamma},\Gamma)}=1\end{subarray}}\inf_{\boldsymbol{j}_{h}\in\boldsymbol{\mathbb{S}}^{p}_{m}(\Gamma)}\|\boldsymbol{j}-\boldsymbol{j}_{h}\|_{\boldsymbol{H}^{-1/2}_{\times}(\operatorname{div}{}_{\Gamma},\Gamma)}.$
(12)

The quadratic convergence order with respect to $\delta_{h,p}$ follows from
the fact that for any eigenfunction $\boldsymbol{j}_{\text{adj}}$ of the
adjoint eigenvalue problem there exists an eigenfunction $\boldsymbol{j}$ of
the eigenvalue problem ${\boldsymbol{\mathcal{V}}_{\kappa}}{}\boldsymbol{j}=0$
such that $\overline{\boldsymbol{j}}=\boldsymbol{j}_{\text{adj}}$, see [16,
Lem. 3]. The estimate (7) together with (11) and (12) yields the final
estimate

$\displaystyle|\kappa-\kappa_{h}|=\mathcal{O}(h^{2(p+1/2)}),$ (13)

for sufficiently smooth surface current densities $\boldsymbol{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
$\boldsymbol{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 $\boldsymbol{j}$ is known to be smooth [14], and even for
other, non-trivial geometries this can be argued.

##  6 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 $\mathbf{T}\colon M\to\mathbb{C}^{m\times m}$ be holomorphic on some
domain $M\subset\mathbb{C}$. Let $\mathbf{T}^{h}$ denote the Hermitian
transpose of $\mathbf{T}$ and $\mathbf{T}^{\prime}$ denote its usual
transpose. We want to solve the nonlinear eigenvalue problem

$\mathbf{T}(z)\mathbf{v}=0,\quad\mathbf{v}\in\mathbb{C}^{m},\quad\mathbf{v}\neq\mathbf{0},\quad
z\in M.$

If $z$ is no eigenvalue of $\mathbf{T}$, we find by the condition above that
the kernel of $\mathbf{T}(z)$ is trivial, thus $\mathbf{T}(z)$ has full rank
and is invertible.

Fix some domain $D\subseteq M$ with boundary $\partial D,$ and assume there to
be $k$ eigenvalues $(\lambda_{j})_{j\leq 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 $\mathbf{T}$ as above and holomorphic $f\colon M\to\mathbb{C}$ it holds
that

$\displaystyle\frac{1}{2\pi i}\int_{\partial
D}f(z)\mathbf{T}(z)^{-1}\mathrm{d}z=\sum^{k}_{j=1}f(\lambda_{j})\mathbf{v}_{h}\mathbf{w}_{h}^{h},$

where $\mathbf{v}_{j}$ and $\mathbf{w}_{j}$ are left and right eigenvectors
corresponding to $\lambda_{n}$ that are normalized according to

$\mathbf{w}_{j}^{h}\mathbf{T}^{\prime}(\lambda_{j})\mathbf{v}_{j}=1$

for all $j\leq k$.

For our purposes, where we set $\mathbf{T}(z)=\mathbf{V}_{z}$ as a matrix-
valued function, cf. (9), the contour integral method can be reduced to the
following theorem.

######  Theorem 8 (Linearisation of Eigenvalue Problems, [12, Thm. 3.1]).

Suppose that $\mathbf{T}\colon M\to\mathbb{C}^{m\times m}$ is holomorphic and
has only simple eigenvalues $(\lambda_{j})_{j\leq k}$ in some connected
subdomain $D\subset M$. Then there exists a diagonalizable matrix $\mathbf{B}$
which can be computed from evaluations of $\mathbf{T}$ such that $\mathbf{B}$
that has the same eigenvalues as the eigenvalue problem under consideration
within $D$.

The proof of this theorem is constructive and corresponds to the following
algorithm, on which our implementation is based.

######  Algorithm 9 (Contour Integral Method, [12, p. 15]).

Let $\mathbf{T}$ be given as above, $\\{t_{n}\\}_{n\leq N}$ be a
discretisation of some boundary $\partial D$ as above, $\delta,\epsilon>0$ and
$\ell<\text{size}(\mathbf{T}(z))$.

Data: $\\{t_{n}\\}_{n\leq N},\mathbf{T},\delta,\epsilon,\ell$

Result: $(\lambda_{j})_{j\leq k}$

1 kfound $\leftarrow$ false;

2 $m\leftarrow$ size($\mathbf{T}$);

3 while  _kfound $==$ false _ do

4 $\hat{\mathbf{V}}\leftarrow\operatorname{RandomFullRank}(m\times\ell)$;

5
$\mathbf{A}_{0}\leftarrow\frac{1}{iN}\sum_{j=1}^{N}\mathbf{(}T(t_{j}))^{-1}\hat{\mathbf{V}}$;

$({\mathbf{V}},{\mathbf{\Sigma}},{\mathbf{W}}^{h})\leftarrow$
SVD$({\mathbf{A}}_{0})$ ;

// $\mathbf{\Sigma}=\text{diag}(\sigma_{1},\dots,\sigma_{\ell})$

6 $k\leftarrow j,$ where
$\sigma_{1}\geq\dots\geq\sigma_{j}>\delta>\sigma_{j+1}\approx\dots\approx\sigma_{\ell}\approx
0$;

7 if  _$k <\ell$_ then

8 kfound $\leftarrow$ true;

9 $\mathbf{V}_{0}\leftarrow{\mathbf{V}}(1\colon m,1\colon k)$;

10 $\mathbf{W}_{0}\leftarrow{\mathbf{W}}(1\colon\ell,1\colon k)$;

11 ${\mathbf{\Sigma}}_{0}\leftarrow\mathbf{\Sigma}$;

12 else

13 $\ell\leftarrow\ell+1$;

14

15$\mathbf{A}_{1}\leftarrow\frac{1}{iN}\sum_{j=1}^{N}t_{j}\mathbf{T}(t_{j})^{-1}\hat{\mathbf{V}}$;

16
$\mathbf{B}\leftarrow\mathbf{V}_{0}^{h}\mathbf{A}_{1}\mathbf{W}_{0}\mathbf{\Sigma}_{0}^{-1}$;

17

18$(\lambda_{j})_{j\leq k}\leftarrow{eigs}(B)$;

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 $\ell_{0}$
of the number of eigenvalues is known such that one can choose
$\ell=\ell_{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 $\mathbf{A}_{1}$ and $\mathbf{A}_{0}$
are assembled simultaneously, since the most expensive operation of the
algorithm is evaluating $\mathbf{T}(t_{j})^{-1}\hat{\mathbf{V}}$. In general,
$k$ and $\ell$ will be small, so storing both $\mathbf{A}_{j}$ poses no
issues. Since smooth contours should be used, exclusively, one can expect the
trapezoidal rules for the assembly of the $\mathbf{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 $\mathbf{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 $\mathbf{V}$ and $\mathbf{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.

##  7 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 open-source [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.

11/21/41/8$10^{-12}$$10^{-9}$$10^{-6}$$10^{-3}$$10^{0}$$h$error$p=1$$p=2$$p=3$$p=4$$\mathcal{O}(h^{3})$$\mathcal{O}(h^{5})$$\mathcal{O}(h^{7})$$\mathcal{O}(h^{9})$

12345678910$10^{-10}$$10^{-7}$$10^{-4}$$10^{-1}$$10^{2}$$i$$\sigma_{i}$$h=1$$h=1/2$$h=1/4$$h=1/8$

Figure 3: On the left the minimal difference of the simultaneously computed
first and second eigenvalue $\lambda_{\mathrm{cim}}$ of the cube to their
analytical solution $\lambda_{i}$ via
$\min_{i=0,1}|\lambda_{\mathrm{cim}}-\lambda_{i}|$. On the right the computed
singular values for the example with $p=3.$

As a first example, we investigate the first two eigenvalues of the unit cube,
where an analytical solution is given by the eigenvalues
$\lambda_{\mathrm{ana},0}=\pi\sqrt{2}$ of multiplicity three and
$\lambda_{\mathrm{ana},1}=\pi\sqrt{3}$ of multiplicity two. The ellipse was
defined as

$\sin(t)+i\cdot 0.05\cdot\cos(t)+5.0,\quad\text{ for }t=[0,2\pi),$

where the contour integrals were evaluated, again, with $N=25.$ The error
visualised corresponds to

$\mathrm{error}=\frac{1}{5}\sum_{0\leq
j<5}\min_{i=0,1}|\lambda_{j}-\lambda_{\mathrm{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 ${\mathcal{O}}(h^{2p+1})$.

11/21/41/8$10^{-11}$$10^{-9}$$10^{-7}$$10^{-5}$$10^{-3}$$10^{-1}$$h$error$p=1$$p=2$$p=3$$p=4$$\mathcal{O}(h^{3})$$\mathcal{O}(h^{5})$$\mathcal{O}(h^{7})$$\mathcal{O}(h^{9})$

12345678910$10^{-11}$$10^{-8}$$10^{-5}$$10^{-2}$$10^{1}$$i$$\sigma_{i}$$h=1$$h=1/2$$h=1/4$$h=1/8$

Figure 4: On the left the mean absolute difference of the computed first
eigenvalue of the sphere to the analytical solution. On the right the computed
singular values for the example with $p=3.$

$2$$2.2$$2.4$$2.6$$2.8$$3$$-4$$-2$$0$$2$$4$$\cdot 10^{-2}$real partimaginary
part

$4$$4.5$$5$$5.5$$6$$-4$$-2$$0$$2$$4$$\cdot 10^{-2}$real partimaginary
partAnalytical eigenvalueComputed eigenvalueQuad. pts. $t_{j}$

Figure 5: The setup for the contour integral method. Sphere (left) and cube
(right), both computed with $N=25$, $p=2$, $h=1/4$.

###  7.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\cdot\sin(t)+i\cdot 0.05\cdot\cos(t)+2.5,\quad\text{ for }t=[0,2\pi).$

The trapezoidal rule for lines 5 and 15 of the algorithm was chosen with
$N=25$. A closed-form solution $\lambda_{\mathrm{ana}}$ is known [16, Sec.
5.1.2] and given as the first root of the spherical Bessel function of the
first kind, cf. Figure 5. It is an eigenvalue of multiplicity three, thus
three non-zero singular values in $\Sigma$ are expected. This behaviour is
reflected by the numerical examples perfectly. The error in Figure 4 refers to
the average error of all three computed eigenvalues, i.e.,

$\mathrm{error}=\frac{1}{3}\sum_{0\leq
j<3}|\lambda_{j}-\lambda_{\mathrm{ana}}|.$

The convergence behaviour w.r.t. $h$ matches the orders
${\mathcal{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\times$) 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.

IGA-BEM | volume-based IGA
---|---
$p$ | $m$ | DOFs | error | $p$ | SD | DOFs | error
3 | 1 | 192 | 2.66e-05 | 3 | 4 | 4350 | 3.345e-05
3 | 2 | 432 | 4.63e-07 | 3 | 8 | 20450 | 3.45e-07
3 | 3 | 1200 | 2.12e-09 | 3 | 14 | 84560 | 1.10e-08
4 | 1 | 300 | 1.62e-06 | 4 | 4 | 6944 | 2.36e-06
4 | 2 | 588 | 3.05e-08 | 4 | 8 | 27280 | 3.94e-09
4 | 3 | 1452 | 3.89e-11 | 4 | 13 | 84560 | 7.88e-11
5 | 1 | 432 | 8.97e-08 | 5 | 4 | 10408 | 1.88e-07
5 | 2 | 768 | 2.04e-09 | 5 | 8 | 35484 | 7.33e-11
5 | 3 | 1728 | 7.17e-12 | 5 | 11 | 69600 | 2.1498e-12
Table 1: Error and degrees of freedom (DOFs) for the approximation of the
smallest eigenvalue of the sphere for different polynomial degrees $p$,
refinement levels $m$ and number of subdivisions (SD). CST Microwave Studio
2018
---
elements | SD | DOFs | error
1089 | 6 | 21597 | 1.8638e-07
5292 | 10 | 101997 | 5.4204e-09
9191 | 12 | 175818 | 1.9524e-09
Table 2: Error and degrees of freedom (DOFs) for the approximation of the
smallest eigenvalue of the sphere with CST Microwave Studio 2018. Settings are
Mesh: Tetrahedral, Desired accuracy: 1e-12, Solver order: 3rd (constant),
Curved elements: up to order 5 (user defined)

###  7.2 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 $\kappa\approx 26.5$. We choose the contour

$0.5\cdot\sin(t)+i\cdot 0.01\cdot\cos(t)+26.5,\quad\text{ for }t=[0,2\pi)$

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
$\lambda_{\mathrm{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.

11/21/41/8$10^{-10}$$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$h$error$p=1$, Bembel
Error$p=1$, CST Error$p=2$, Bembel Error$p=2$, CST Error$p=3$, Bembel
Error$p=3$, CST Error$p=4$, Bembel Error$p=4$, CST Error Figure 6: Eigenvalue
problem of the one-cell TESLA cavity solved for the first eigenvalue. Error
has to be understood as the relative difference to the result of a
computation. Bembel Reference refers to the error of our method to a reference
computation with $h=1/8$ and $p=5$ of our own implementation, and the CST
Reference to the error w.r.t. a reference solution obtained via CST Microwave
Studio 2018.

##  8 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
${\mathcal{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.

> Acknowledgement: The work of Gerhard Unger is supported by the Austrian
> Science Fund (FWF), project P31264. The work of Felix Wolf is supported by
> DFG Grants SCHO1562/3-1 and KU1553/4-1, the Excellence Initiative of the
> German Federal and State Governments and the Graduate School of
> Computational Engineering at TU Darmstadt.

## References

  * [1] B. Aune, R. Bandelmann, D. Bloess, B. Bonin, A. Bosotti, M. Champion, C. Crawford, G. Deppe, B. Dwersteg, D. A. Edwards, H. T. Edwards, M. Ferrario, M. Fouaidy, P.-D. Gall, A. Gamp, A. Gössel, J. Graber, D. Hubert, M. Hüning, M. Juillard, T. Junquera, H. Kaiser, G. Kreps, M. Kuchnir, R. Lange, M. Leenen, M. Liepe, L. Lilje, A. Matheisen, W.-D. Möller, A. Mosnier, H. Padamsee, C. Pagani, M. Pekeler, H.-B. Peters, O. Peters, D. Proch, K. Rehlich, D. Reschke, H. Safa, T. Schilcher, P. Schmüser, J. Sekutowicz, S. Simrock, W. Singer, M. Tigner, D. Trines, K. Twarowski, G. Weichert, J. Weisend, J. Wojtkiewicz, S. Wolff, K. Zapfe, Superconducting tesla cavities, Phys. Rev. Accel. Beams 3 (2000) 092001.
  * [2] K. Halbach, R. F. Holsinger, SUPERFISH – a computer program for evaluation of RF cavities with cylindrical symmetry, Part. Accel. 7 (1976) 213–222.
  * [3] U. van Rienen, T. Weiland, Triangular discretization method for the evaluation of RF-fields in cylindrically symmetric cavities, IEEE T. Magn. 21 (6) (1985) 2317–2320.
  * [4] T. Weiland, On the unique numerical solution of Maxwellian eigenvalue problems in three dimensions, Part. Accel. 17 (227–242) (1985).
  * [5] F. Wolf, Figures: CAD/CAE and Mesh Generation (2020).  doi:10.6084/m9.figshare.11590698.v1.
  * [6] N. Georg, W. Ackermann, J. Corno, S. Schöps, Uncertainty quantification for Maxwell’s eigenproblem based on isogeometric analysis and mode tracking, Comput. Method. Appl. M. 350 (2019) 228–244.
  * [7] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Method. Appl. M. 194 (2005) 4135–4195.
  * [8] A. Buffa, G. Sangalli, R. Vázquez, Isogeometric methods for computational electromagnetics: B-spline and T-spline discretizations, J. Comput. Phys. 257, Part B (0) (2013) 1291–1320.
  * [9] Z. Bontinck, J. Corno, H. De Gersem, S. Kurz, A. Pels, S. Schöps, F. Wolf, C. de Falco, J. Dölz, R. Vázquez, U. Römer, Recent advances of isogeometric analysis in computational electromagnetics, ICS Newsletter (International Compumag Society) 3, preprint available: arXiv:1709.06004 (2017).
  * [10] J. Corno, C. de Falco, H. De Gersem, S. Schöps, Isogeometric simulation of Lorentz detuning in superconducting accelerator cavities, Comput. Phys. Commun. 201 (2016) 1–7.
  * [11] J. Dölz, S. Kurz, S. Schöps, F. Wolf, Isogeometric boundary elements in electromagnetism: Rigorous analysis, fast methods, and examples, SIAM J. Sci. Comput. 41 (5) (2019) B983–B1010.
  * [12] W.-J. Beyn, An integral method for solving nonlinear eigenvalue problems, Linear Algebra Appl. 436 (10) (2012) 3839–3863, special Issue dedicated to Heinrich Voss’s 65th birthday.
  * [13] J. Asakura, T. Sakurai, H. Tadano, T. Ikegami, K. Kimura, A numerical method for nonlinear eigenvalue problems using contour integrals, JSIAM Letters 1 (2009) 52–55.
  * [14] C. Wieners, J. Xin, Boundary element approximation for Maxwell’s eigenvalue problem, Math. Methods Appl. Sci. 36 (18) (2013) 2524–2539.
  * [15] A. Buffa, J. Dölz, S. Kurz, S. Schöps, R. Vázquez, F. Wolf, Multipatch approximation of the de Rham sequence and its traces in isogeometric analysis, Numer. Math. 144 (1) (2020) 201–236.
  * [16] G. Unger, Convergence analysis of a Galerkin boundary element method for electromagnetic eigenvalue problems, Technical report 2017/2, Institute of Computational Mathematics, Graz University of Technology (2017).
  * [17] A. Buffa, J. Rivas, G. Sangalli, R. Vázquez, Isogeometric discrete differential forms in three dimensions, SIAM J. Numer. Anal. 49 (2) (2011) 818–844.
  * [18] J. D. Jackson, Classical Electrodynamics, 3rd Edition, Wiley and Sons, New York, 1998.
  * [19] P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford University Press, Oxford, 2003.
  * [20] W. McLean, Strongly elliptic systems and boundary integral equations, Cambridge University Press, Cambridge,United Kingdom, 2000.
  * [21] A. Buffa, R. Hiptmair, Galerkin boundary element methods for electromagnetic scattering, in: M. Ainsworth, P. Davies, D. Duncan, B. Rynne, P. Martin (Eds.), Topics in computational wave propagation, Springer, Berlin Heidelberg, 2003, pp. 83–124.
  * [22] J. Dölz, H. Harbrecht, S. Kurz, S. Schöps, F. Wolf, A fast isogeometric BEM for the three dimensional Laplace- and Helmholtz problems, Comput. Method. Appl. M. 330 (Supplement C) (2018) 83–101.
  * [23] A. Buffa, G. Sangalli, R. Vázquez, Isogeometric analysis in electromagnetics: B-splines approximation, Comput. Method. Appl. M. 199 (2010) 1143–1152.
  * [24] L. Beirão da Veiga, A. Buffa, G. Sangalli, R. Vázquez, Mathematical analysis of variational isogeometric methods, Acta Numer. 23 (2014) 157–287.
  * [25] L. Piegl, W. Tiller, The NURBS Book, 2nd Edition, Springer, Berlin Heidelberg, 1997.
  * [26] A. F. Peterson, Mapped vector basis functions for electromagnetic integral equations, Synthesis Lectures on Computational Electromagnetics 1 (1) (2006) 1–124.
  * [27] J. Dölz, S. Kurz, S. Schöps, F. Wolf, A numerical comparison of an isogeometric and a parametric higher order Raviart-Thomas approach to the electric field integral equation, IEEE T. Antenn. Propag. 68 (1) (2020) 593–597.
  * [28] A. Bespalov, N. Heuer, R. Hiptmair, Convergence of the natural hp-BEM for the electric field integral equation on polyhedral surfaces, SIAM J. Numer. Anal. 48 (2010) 1518–1529.
  * [29] M. Halla, Regular Galerkin approximations of holomorphic T-Gårding operator eigenvalue problems, ASC Report 4, TU Wien (2016).
  * [30] O. Karma, Approximation in eigenvalue problems for holomorphic Fredholm operator functions. I, Numer. Funct. Anal. Optim. 17 (3-4) (1996) 365–387.
  * [31] O. Karma, Approximation in eigenvalue problems for holomorphic Fredholm operator functions. II. (Convergence rate), Numer. Funct. Anal. Optim. 17 (3-4) (1996) 389–408.
  * [32] A. Imakura, L. Du, T. Sakurai, Relationships among contour integral-based methods for solving generalized eigenvalue problems, Jpn. J. Ind. Appl. Math. 33 (3) (2016) 721–750.
  * [33] R. Mennicken, M. Möller, Non-Self-Adjoint Boundary Eigenvalue Problems, 1st Edition, Vol. 192 of Mathematics Studies, Elsevier, 2003.
  * [34] J. Dölz, H. H., S. Kurz, M. Multerer, S. Schöps, F. Wolf, Bembel – the C++ boundary element engineering library, official website, www.bembel.eu. Date of access August 20, 2019.
  * [35] J. Dölz, H. Harbrecht, S. Kurz, M. Multerer, S. Schöps, F. Wolf, Bembel: The fast isogeometric boundary element C++ library for Laplace, Helmholtz, and electric wave equation, Technical report, arxiv:1906.00785 (2019).
  * [36] F. Wolf, Fork of Bembel, will be made available upon acceptance of the manuscript.
  * [37] J. Benoît, G. Guennebaud, Eigen3 C++ linear algebra template library, official website, eigen.tuxfamily.org. Date of access August 20, 2019.
  * [38] R. Vázquez, A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0, Comput. Math. Appl. 72 (3) (2016) 523–554.
  * [39] C. de Falco, A. Reali, R. Vázquez, GeoPDEs: A research tool for isogeometric analysis of PDEs, Advances in Engineering Software 42 (2011) 1020–1034.