# Uncertainty Quantification for Maxwell's Eigenproblem using Isogeometric Analysis

The electromagnetic field distribution as well as the resonating frequency of various modes in superconducting cavities used in particle accelerators for example are sensitive to small geometry deformations. The occurring variations are motivated by measurements of an available set of resonators from which we propose to extract a small number of relevant and independent deformations by using a truncated Karhunen-Loève expansion. The random deformations are used in an expressive uncertainty quantification workflow to determine the sensitivity of the eigenmodes. For the propagation of uncertainty, a stochastic collocation method based on sparse grids is employed. It requires the repeated solution of Maxwell's eigenvalue problem at predefined collocation points, i.e., for cavities with perturbed geometry. The main contribution of the paper is ensuring the consistency of the solution, i.e., matching the eigenpairs, among the various eigenvalue problems at the stochastic collocation points. To this end, a classical eigenvalue tracking technique is proposed that is based on homotopies between collocation points and a Newton-based eigenvalue solver. The approach can be efficiently parallelized while tracking the eigenpairs. In this paper, we propose the application of isogeometric analysis since it allows for the exact description of the geometrical domains with respect to common computer-aided design kernels, for a straightforward and convenient way of handling geometrical variations and smooth solutions.

## Authors

• 7 publications
• 2 publications
• 8 publications
• 61 publications
• ### Uncertainty Quantification for Geometry Deformations of Superconducting Cavities using Eigenvalue Tracking

The electromagnetic field distribution as well as the resonating frequen...
02/08/2018 ∙ by Niklas Georg, et al. ∙ 0

• ### Uncertainty quantification in imaging and automatic horizon tracking: a Bayesian deep-prior based approach

In inverse problems, uncertainty quantification (UQ) deals with a probab...
04/01/2020 ∙ by Ali Siahkoohi, et al. ∙ 0

• ### Isogeometric Analysis Simulation of TESLA Cavities Under Uncertainty

In the design of electromagnetic devices the accurate representation of ...
11/06/2017 ∙ by Jacopo Corno, et al. ∙ 0

• ### Fast approximation by periodic kernel-based lattice-point interpolation with application in uncertainty quantification

This paper deals with the kernel-based approximation of a multivariate p...
07/13/2020 ∙ by Vesa Kaarnioja, et al. ∙ 0

• ### Uncertainty Modeling and Analysis of the European X-ray Free Electron Laser Cavities Manufacturing Process

This paper reports on comprehensive efforts on uncertainty quantificatio...
06/21/2019 ∙ by Jacopo Corno, et al. ∙ 0

• ### Uncertainty Modeling and Analysis of the European XFEL Cavities Manufacturing Process

This paper reports on comprehensive efforts on uncertainty quantificatio...
06/21/2019 ∙ by Jacopo Corno, et al. ∙ 0

• ### Uncertainty quantification of microstructure variability and mechanical behaviour of additively manufactured lattice structures

Process-induced defects are the leading cause of discrepancies between a...
03/17/2021 ∙ by Nina Korshunova, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Motivation and Introduction

In rf resonators, the electromagnetic fields oscillate at certain frequencies with higher magnitude than at others. The advantageous operation of such devices at selected eigenfrequencies is technically profitable used for example to set up dedicated filter components or to provide field amplitudes, which are higher than the magnitude of the excitation. In the demanding context of charged particle accelerators for example, specifically designed resonators, which are also known as cavities, can efficiently accelerate bunched charged particle beams. Essentially, the acceleration is performed in the longitudinal direction only, although specific devices can also be developed to deflect the beam in a transverse direction. The shape of the resonating structure, as well as the walls materials, determine, to a large extent, the overall performance of the entire system. As a consequence, plenty of effort has been put in the development of proper simulation tools to computationally assist scientists and engineers during the complex design process.

Because of the limitations in the available computational resources, in the early days of such simulations the main focus was exclusively put on lossless two-dimensional (2D) eigenanalysis. Meanwhile, even expensive three-dimensional (3D) models including dissipative materials or leaky boundary conditions can be handled routinely. Due to the strong dependence of the solution with respect to the geometry, a huge amount of thought and work has been spent on the precise modeling of the resonator’s surface and the corresponding boundary conditions. Various discretization schemes have been applied to introduce proper computational elements ranging from a flexible irregular triangular mesh Halbach_1976aa ; Rienen_1985aa to an efficient structured rectangular variant Weiland_1984aa ; Weiland_1985aa . In both cases, special care has been put on the accurate description of the boundary surface to retain the specific convergence order of the underlying method. A further step in this direction can be found in the application of iga Hughes_2005aa , e.g. eigenvalue problems Hughes_2014aa , where the geometric modeling is naturally linked to the utilized cad kernels by using nurbs. Thus, no geometry-related modeling error is introduced if the cad model can be described by nurbs.

Another advantage of iga is the inter-element smoothness. Quantities of interest, for example the electromagnetic fields in the vicinity of the main axis which are important for subsequent particle simulations, can be obtained with high accuracy Corno_2016aa

. Furthermore, IGA promises faster convergence, i.e., less degrees of freedom are needed for the same accuracy, when compared to traditional polynomial approaches. As an illustration, we compare in Fig.

1 the accuracy of IGA discretizations of Maxwell’s eigenproblem on the unit square to Nédélec elements Monk_2003aa ; nedelec1980 . It can be clearly observed that higher regularity across elements leads to smaller numerical errors for smooth solutions. It shall be noted that, in both cases, there is no geometry-related modeling error. A similar comparison for a cylindrical cavity, i.e. a curved geometry, can be found in Corno_2016aa . In this paper, it is also demonstrated that iga is beneficial in terms of computational time.

The assembled eigenvalue formulation requires the solution of individual eigenpairs, which are also known as modes. They are characterized by a specific field distribution together with the corresponding eigenfrequency. Whereas the distribution of the electromagnetic field is connected to the determined eigenvector, the extracted eigenvalue specifies the resonating frequency and quality factor if available. In the lossy case, the modal analysis is obviously of complex nature, while in the simplified lossless approximation all calculations can be handled with real arithmetic only. One critical aspect for the manufacturing of the designed cavities is the specification of appropriate tolerances for the shape of the resonator to ensure a proper operability performance. While the sensitivity of the examined mode with respect to unavoidable geometric variations can be obtained with the help of the known sensitivity analysis, e.g.

Gorgizadeh_2018aa , a profound understanding of the corresponding distributions is obtained using the more expensive uq, see e.g. Xiao_2007aa ; Heller_2014aa . To our best knowledge, previous works handle only tiny deformations and represent them rather naively, e.g. using non-curved elements, deforming only one layer of elements or remesh the problem without considering the mesh sensitivity. All these drawbacks are elegantly circumvented by the workflow proposed in this paper.

To study the influence of geometric boundary variations on the distribution of the electromagnetic fields and the related eigenfrequencies, the introduction of proper shape parameters has to be performed. For each cavity, which is related to a specific shape, the respective eigensolutions can be calculated independently of each other. Whereas the assignment of the simulation results with respect to the examined modes is simple for non-overlapping functions, it becomes much worse for touching or crossing graphs. In those cases, a direct evaluation of the derivatives with respect to the shape parameter is not suitable anymore and reliable eigenvalue-tracking methods have to be applied instead. In the following chapters, the necessary steps to establish a uq workflow for electromagnetic fields in superconducting cavities using eigenvalue tracking will be formulated.

In the following section we state Maxwell’s eigenvalue problem and introduce the isogeometric discretization chosen for the approximation of the fields. In section 3 we present the uncertainty quantification methods used for the simulation and their relation to the geometry parametrization. Given the necessity to properly categorize the eigensolutions for each uncertain configuration, we present the eigenvalue tracking method in section 4. Finally, we show the application of the proposed scheme to the simple case of the pillbox cavity and to a real world example, i.e. the tesla cavity (see Fig. 1(a) and 1(b)). In the last section we give some concluding remarks.

## 2 Numerical Approximation of Maxwell’s Eigenproblem

### 2.1 Maxwell’s Eigenvalue Problem

Given a bounded and simply connected domain , with Lipschitz continuous boundary , the electromagnetic fields in the source-free, time-harmonic case, are governed by Maxwell’s equations:

 ∇×E =−jωμ0H in Ω (1) ∇×H =jωε0E in Ω ∇⋅ε0E =0 in Ω ∇⋅μ0H =0 in Ω,

where and denote the electric and magnetic field strength, respectively, and we assume the electric permittivity and magnetic permeability of vacuum. We further apply pec boundary conditions

 E×n =0 on ∂Ω (2) H⋅n =0 on ∂Ω.

Starting from Maxwell’s equations, one can derive the classical eigenvalue formulation by eliminating from (1): find the wave number and such that

 ∇×(∇×E) =k2E in Ω (3) ∇⋅E =0 in Ω E×n =0 on ∂Ω.

The standard variational formulation of (3) is: find with such that

 (∇×E,∇×v)=k2(E,v)∀v∈H0(curl;Ω), (4)

where we introduced the function space

of square-integrable vector fields with square-integrable

curl. For further information on function spaces in the context of Maxwell’s equations, the reader is referred to Monk_2003aa . The subscript in indicates vanishing tangential components on the boundary since we are imposing pec boundary conditions.

In order to be able to numerically solve (4) we deduce its discrete counterpart: following a Galerkin procedure, we introduce a sequence of finite-dimensional spaces and we restrict (4) to it: find with such that

 (∇×Eh,∇×v)=k2h(Eh,v)∀v∈Vh. (5)

Whatever the choice of the functional space , can be expressed in terms of a set of basis functions :

 Eh=Ndof∑j=1ejvj. (6)

The discrete solution of Maxwell’s eigenvalue problem is then obtained by solving the generalized eigenvalue problem

 Ke=k2hMe, (7)

where and are the stiffness and mass matrix for vacuum respectively, given by

 Kij=(∇×vj,∇×vi),Mij=(vj,vi). (8)

There exist multiple choices for the space and its basis although care must be taken to ensure a proper approximation of . For the fem case we refer the interested reader to Boffi_2010aa . In the following section we will follow Buffa_2010aa and introduce a curl-conforming iga discretization.

### 2.2 Isogeometric Analysis

iga is a technique that was introduced by Hughes et al. Hughes_2005aa with the aim of allowing the exact representation of geometries defined via cad software, independently of the level of refinement of the computational grid. It is clear that, since the geometry directly affects the resonant properties of rf cavities, the exact description of the simulation domain is strongly desired, even more so in the case of shape sensitivity analysis and uq.

In order to achieve this goal, the same classes of basis functions that are commonly used for geometry description in cad software, are used for the representation of solution vector fields. Such basis functions are the so called B-splines and nurbs. For the definition and construction of such basis, we refer the interested reader to de-Boor_2001aa ; Piegl_1997aa .

The use of iga in cavity simulation has been shown to be beneficial both in terms of accuracy and of overall reduction of the computational cost Corno_2016aa ; Corno_2015ac . For an overview of the applicability of iga to other electromagnetic problems, we refer to Bontinck_2017ag .

We assume that the domain can be exactly described through a smooth mapping obtained from the cad software, see Fig. 3:

 F0:^Ω→Ω, (9)

where is a parametric domain (e.g. the unit cube). We furthermore assume that is piecewise smoothly invertible.

We define the discretization space in a general way

 Vh :={v∈H(curl;Ω):^v=ι(v)∈^Vh} (10) ≡{v∈H(curl;Ω):v=ι−1(^v),^v∈^Vh},

where is a discrete space on the parametric domain and is a so-called pull-back function defined from the parametrization . Let be a basis for , the set is a basis for that can be substituted in (7)-(8) to obtain a finite-dimensional generalized eigenvalue problem.

We indicate as the standard space of B-Splines of degree along dimension , with regularity , constructed on the reference domain . Following Buffa_2010aa the discrete space is chosen as

 ^Vh:=Sp1−1,p2,p3α1−1,α2,α3×Sp1,p2−1,p3α1,α2−1,α3×Sp1,p2,p3−1α1,α2,α3−1. (11)

It can be proven that basis functions in this space obey to a discrete de Rham sequence, thus being suitable candidates for the discretization of .

The discrete space in the physical domain is obtained through a curl-conforming ‘Piola’ transformation Monk_2003aa :

 Vh:={v:(dF0)⊤(v∘F−10)∈^Vh}, (12)

where we indicate with the Jacobian of the parametrization.

### 2.3 Domain Parametrization

Given a B-Spline (or nurbs) basis of , the mapping is constructed using a set of control points in the physical space as:

 F=∑iPiBi. (13)

The points represent the so-called control mesh.

It is straightforward to show Piegl_1997aa that, given definition (13), cad objects are invariant under affine transformations such as translations, scalings and rotations, i.e. to apply the transformation to the object it is sufficient to apply it to the control mesh. Furthermore, is smooth with respect to the movement of the control points.

Let us now consider the domain and its parametrization to be given in terms of a -vector of parameters , i.e. becomes and the mapping

 Fδ:^Ω→Ω(δ). (14)

The shape variability can be reinterpreted as the variability of the control mesh:

 Fδ=∑iPi(δ)Bi, (15)

and by following the steps described in the previous section in the case of a parametrized domain, one derives the parametrized generalized eigenvalue problem:

 K(δ)e(δ)=k2h(δ)M(δ)e(δ) (16)

where denotes the numerical approximation of the eigenvalue for each parameter . Suppressing for brevity the dependency, we denote the discrete eigenpairs with where . In the following, we assume the parametric dependency of the discrete eigenpairs and their continuous counterparts on to be sufficiently smooth.

## 3 Uncertainty Quantification

In operation, the geometry of a particular TESLA cavity will differ from the design shown in Fig. 1(b): the manufacturing process, tuning and multiphysical effects as cooling, Lorentz detuning and microphoning will introduce geometry variations Xiao_2007aa ; Deryckere_2012aa ; Schmidt_2014aa ; Corno_2016aa . In the simplest case, these variations can be expressed in terms of the design parameters, see Fig. 1(b), i.e.,

 δdesign=[a1,a2,b1,b2,L,R%iris,Req]⊤ (17)

and lead to an explicit definition of the parametrized geometry mapping . However, in more realistic scenarios the variation of a particular geometry is unknown.

Let us assume that we have measured times the geometric variations and that each observation can be described by an -dimensional vector

 x(m)=[x(m)1,x(m)2,...,x(m)N]⊤,m=1,…,M. (18)

And let us further assume that all components

are realizations of normally distributed random variables

 Xn∼N(μn,σ2n),n=1,…,N, (19)

possibly mutually correlated. Expectation and variance are denoted by

and respectively and can be obtained from the measurement samples.

We propose the discrete Karhunen–Loève expansion (or equivalently principal component analysis

Xiu_2010aa ) to convert the observations into a set of linearly uncorrelated variables. uq methods can then be applied to investigate the sensitivity of the eigenmodes to those variables.

### 3.1 Karhunen–Loève

Let us collect all the observations in the matrix

 T=[x(1),x(2),…,x(M)]⊤, (20)

whose columns are the random variables and whose rows are the measured observations. We will denote by its covariance matrix

 C:=1M−1M∑m=1(x(m)−μ)(x(m)−μ)⊤, (21)

where is the sample mean

 μ=1MM∑m=1x(m). (22)

We search for an eigendecomposition of the covariance matrix of the type

 C=VΣV⊤, (23)

where the columns of are the orthonormal right eigenvectors of , and is a diagonal matrix whose elements are the corresponding eigenvalues. Without loss of generality, we will further assume the columns of to be ordered such that .

Now, only the most significant contributions are extracted by approximating

 C≈Ct:=VtΣtV⊤t, (24)

where and is a )-matrix with the corresponding eigenvectors.

Finally, the variations can be approximated by a low-dimensional representation DElia_2013aa

 ⎡⎢ ⎢⎣X1⋮XN⎤⎥ ⎥⎦≈⎡⎢ ⎢⎣μ1⋮μN⎤⎥ ⎥⎦+VtΣ1/2t⎡⎢ ⎢⎣δ1⋮δNt⎤⎥ ⎥⎦. (25)

that is based on normally distributed parameters . With respect to the original vector of parameters , may not have a physical meaning, but has a significantly smaller dimension. Reducing the size of the parameter space has a great impact in the efficiency of a stochastic collocation method.

Let us now consider the case of a nurbs domain specified by the mapping (13). Given , to describe the deformed geometry (14) we will firstly construct the

B-Spline interpolating curve mappings

de-Boor_2001aa

 Fμ =∑iPμiBi (26) Fj =∑iPjiBi∀j=1,…,Nt% , (27)

where interpolates the mean deformation vector and is an interpolating curve through the points given by the -th column of the matrix . Unique interpolations can be obtained by requiring additional properties as a certain polynomial degree and minimal curvature Piegl_1997aa . Both and the share the same basis functions as the geometry mapping and differ only in their control mesh. The control points defining the deformed geometry are then given by

 Pi(δ)=Pi+Pμi+Nt∑j=1δjPji. (28)

### 3.2 Stochastic Collocation

Quantities of interest are statistical measures about the eigenmodes with respect to the influence of the vector of random variables. In particular, we aim for expectation and variance of the -th eigenfrequency ,

 E[fj] =∫∞−∞fj(δ)ρ(δ)dδ (29a) var[fj] =∫∞−∞(fj(δ)−E[fj])2ρ(δ)dδ, (29b)

where

denotes the joint probability density of

.

In order to approximate (29) numerically, a stochastic collocation method Babuska_2010aa is applied which can be summarized as follows: based on the assumption that the parametric dependency is smooth, we construct a polynomial approximation based on Lagrangian interpolation

 fj(δ)≈Nk∑k=1fh,j(δk)lk(δ), lm(δn)={1if m=n0if m≠n (30)

where are multivariate Lagrange polynomials constructed by tensorization of univariate Lagrange polynomials and is the set of corresponding collocation points (tensor grid). The eigenfrequencies

 fh,j(δk)=kh,j(δk)2π√ε0μ0. (31)

are obtained by solving (16) for each set of discrete parameters . This requires rerunning the simulation times for slightly different geometries. No modifications of the code are necessary and hence the method can be considered as non-intrusive. Finally, inserting (30) into (29) yields, for the -th eigenfrequency:

 E[fj] ≈Nk∑k=1fh,j(δk)wk (32a) var[fj] ≈Nk∑k=1(fh,j(δk)−E[fj])2wk, (32b)

with weights

 wk=∫∞−∞lk(δ)ρ(δ)dδ. (33)

The integrals can be evaluated exactly using Gauss-Hermite quadrature for each independent parameter . The quadrature nodes are also used as collocation points .

The stochastic collocation method (32

) based on tensor grids is known to be very efficient for a small number of random variables but suffers strongly from the curse of dimensionality

Nobile_2008aa . The number of collocation points grows exponentially with the dimension of the parametric space . This problem can be mitigated by the sparse grid approach Bungartz_2004aa ,Smolyak_1963aa . Linear combinations of tensor grids with different degrees for different parameters lead to more efficient interpolation and cubature rules for a (moderately) large stochastic space Novak_1999aa . An illustration of stochastic collocation grids is given in Fig. 4.

In each of the collocation points, eigenvalue problem (7) is solved. The natural ordering of the eigenpairs according to their eigenvalue, however, might not be sufficient to ensure consistency since even small changes in the domain can cause eigenvalue crossing or separation/degeneration of modes. To be able to correctly place the frequencies for each parametric point , an eigenvalue tracking algorithm is applied.

## 4 Eigenvalue Tracking

The following procedure is adapted from Lui_1997aa . To be able to track the eigenvalues, we construct an algebraic homotopy between the quadrature points in the parameter space in such a way that:

 Kk,t :=(1−t)K(δ0)+tK(δk) (34a) Mk,t :=(1−t)M(δ0)+tM(δk) (34b)

where and is a given start configuration in the parameter space (e.g. corresponding to an unperturbed geometry). We further assume that at the modes are not degenerated and that there is no bifurcation from to . Other homotopies could be analogously defined, e.g. such that the eigenvalue problems for correspond to actual geometries. However, the algebraically motivated choice (34) allows the straightforward definition of the derivatives,

 K′k,t :=∂∂tKk,t=K(δk)−K(δ0) (35a) M′k,t :=∂∂tMk,t=M(δk)−M(δ0) (35b)

Now, all the quantities depending on the vector , can then be viewed as depending on the scalar quantity which makes the following one-dimensional analysis feasible. For each collocation point and eigenpair , the eigenvalue problem (7) can be rewritten as a linear system of equations with a normalization constraint using a suitable vector

 [Ktet−λtMtetc⊤et−1]=[00], (36)

where and . The subscripts and addressing eigenpair number and collocation point are suppressed for easier readability.

Let us consider a Taylor expansion at along

 [~et~λt]=[e0λ0]+(t−t0)[e′0λ′0] (37)

to obtain a first-order approximation for the eigenpair at . The derivatives of the eigenvalues and eigenvectors are obtained from (36) which is differentiated with respect to the parameter . The resulting linear system

 [K0−λ0M0−M0e0c⊤0][e′0λ′0]=[−K′0e0+λ0M′0e00] (38)

is then solved with respect to the eigenvalue and eigenvector derivatives. The predicted eigenpair can be used as an initial guess

 [e(0)tλ(0)t]:=[~et~λt] (39)

for the Newton-Raphson method in order to obtain the exact eigenpair . In the -th iteration of the Newton-Raphson method, the linear system of equations to be solved is

 [Kt−λ(i−1)tMt−Mte(i−1)tc⊤0][Δe(i)tΔλ(i)t]= −[Kte(i−1)t−λ(i−1)tMte(i−1)tc⊤e(i−1)t−1]. (40)

The increments and improve the prediction by assigning

 [e(i)tλ(i)t]:=[e(i−1)tλ(i−1)t]+[Δe(i)tΔλ(i)t]. (41)

As a termination criterion for the Newton-Raphson method, the norm of the residual of eigenvalue problem (36) can be considered as well as the norm of . The required number of Newton iterations can then be used for a stepsize control in order to increase the efficiency and robustness of the algorithm. The proposed tracking algorithm is illustrated by the flowchart in Fig. 5.

It shall be noted that the assumption of non-degeneracy of the eigenmodes at is made for simplicity only. The algorithm could be generalized in order to address the case of an intersection of eigenvalues at by using Ojalvo’s method ojalvo ; dailey for computing the derivatives instead of (38). The numerical efficiency of the presented tracking method could also be further improved, e.g. by using a simplified Newton method Deuflhard_2004aa , by employing a higher order Taylor approximation in the prediction step (37) or by tracking a subspace of eigenpairs and matching the eigenpairs based on a correlation coefficient Jorkowski_2018aa ; Yang_2007ab .

## 5 Application

In this section, we present the application of the eigenvalue tracking method to cavities with uncertain geometries. As a clear indication of the necessity of the tracking, we first analyze the test case of a cylindrical cavity with uncertain radius. Stochastic collocation is then applied to study the sensitivity of a tesla type cavity Aune_2000aa to eccentric deformation of its 9 cells. In both cases the discretization is carried out with iga using GeoPDEs and NURBS Octave packages Falco_2011aa ; nurbs .

### 5.1 Pillbox Cavity

As a numerical example with an available closed-form reference solution, we consider a pillbox (i.e. cylindrical) cavity with uncertain radius. We consider a fixed length

and we assume a uniformly distributed radius

. Second degree basis functions are chosen for a total of degrees of freedom. This choice is a compromise between increased regularity () and sparseness and condition of the system matrix. The pillbox cavity as well as the electric and magnetic field lines of the fundamental tm mode are shown in Fig. 6.

#### 5.1.1 Tracking

We use the proposed eigenvalue tracking technique from section 4 to track the eigensolutions over the parameter range . The exact eigenfrequencies are well known Jackson_1998aa . At a radius of approximately a transition of the fundamental mode from te to tm is observed.

To illustrate the relevance of eigenvalue tracking, the eigenvalues obtained by solving (16) at discrete sample points are shown in Fig. 6(a)

. Without tracking, the eigenvalues cannot be associated to specific modes without cumbersome post-processing procedures that typically require some heuristics

Brackebusch_2014aa .

For this first example, the algebraic homotopy (34) is replaced by the physical system matrices such that all eigenvalues computed along the homotopy path correspond to the solutions of a physical problem. This allows for the comparison with the closed-form solution with respect to the radius . In practice, one is probably not interested in the intermediate points. However, in this test case the derivatives cannot be straightforwardly computed and are then approximated by finite differences. Fig. 6(b) demonstrates the result of the tracked eigenvalues.

Remark: In section 4 we assumed, for simplicity, that the eigenmodes at are not degenerate. Obviously, this assumption is not fulfilled for this test case. However, the proposed algorithm can still be applied since the degeneracy of the respective eigenvalues is caused by symmetry of the geometry which is preserved for all parameter changes.

#### 5.1.2 Uncertainty Quantification

We determine the uncertainty in the 6 lowest eigenfrequencies at for the uniformly distributed radius using stochastic collocation based on Clenshaw-Curtis quadrature with collocation points. The eigenvalue problem for the first collocation point at is solved based on Arnoldi’s method Lehoucq_1996aa using the command eigs in Matlab R2017a. For each eigenpair the tracking method is applied for every remaining collocation point. The initial stepsize is chosen such that all points are reached using only one step. For this particular model, the Newton scheme always converges, using on average 2.2 iterations (max. 4) such that the stepsize is never reduced.

The first two statistical moments can be seen in Fig.

8. The relative errors in the numerical results with respect to the analytic reference solutions are below

(for expectations as well as standard deviations).

#### 5.1.3 Computational Cost

The tracking is necessary to obtain correct results since the eigenmodes need to be matched at each collocation point. However, as shown in this subsection, it is also numerically more efficient than naively solving the repeated eigenvalue problems independently of each other. The computations are performed on a workstation with an Intel(R) Xeon(R) CPU E5-2687W processor and 256 GB RAM. All time measurements are repeated times and averaged. It shall be noted that, for a fair comparison, in all eigenvalue solvers the occurring linear systems of equations were solved with the Matlab backslash operator \, i.e. sparse Gaussian elimination. Improvements are to be expected when reusing factorizations or preconditioners with appropriate iterative solvers.

As a reference, we solve the eigenvalue problem at each collocation point using a Matlab implementation of the Jacobi-Davidson algorithm Sleijpen_1996aa ; Fokkema_1998aa as well as the Matlab command eigs which internally calls the Arpack library Lehoucq_1998aa . Since we are interested in quantifying uncertainties for 10 eigenmodes and given that the ordering of eigenvalues may change, we calculate 20 eigenvalues to ensure that the corresponding eigenvalues are included at each collocation point. The Jacobi-Davidson algorithm solves the eigenproblems on average in iterations and while eigs takes on average for solving linear systems.

The computational cost of the proposed tracking approach is determined by the time needed for solving (38) in order to compute the derivatives and for solving the linear systems of the Newton scheme (40). Cost for assembly of the system matrices is excluded since the algebraic mapping (35) can be employed for uq which does not require any additional matrix assemblies. On average, the computational time needed for solving linear systems for each collocation point and eigenpair is .

#### 5.1.4 Shape morphing

In order to illustrate the ability of tracking the eigenmodes during more sophisticated deformations, we consider the shape morphing from the 1-cell TESLA cavity to a pillbox cavity (radius , length ), as illustrated in Fig. 9. The eigenmodes in each cavity are discretized using second degree basis functions and degrees of freedom.

For this test case, we consider on the one hand the algebraic homotopy (34) and on the other hand the physical system matrices corresponding to the depicted geometries . In both cases, we apply the proposed tracking algorithm for the first few eigenvalues of the TESLA cavity. Fig. 10 demonstrates the results of the tracked eigenvalues. It can be observed that both approaches match the same eigenmodes.

### 5.2 TESLA Cavity

A more challenging application is the simulation of a so-called tesla type cavity. This is an rf superconducting cavity for linear accelerators with a design accelerating gradient of . In particular, we focus on the ttf design of the cavities TTF_1995 which were installed at desy in Hamburg.

The TESLA cavity is a 9-cell standing wave structure of about in length, whose accelerating mode resonates at . The cavities are produced, one half-cell at a time, by deep drawing of pure niobium sheets TTF_1995 . The half cells are then connected together by ebm. The production process inherently introduces geometry imperfections that need to be taken into account in operation. Here, we focus on the misalignment of each of the 9 cells with respect to the ideal axis of the cavity. The misalignment can produce undesired effects on the field quality and on the resonating frequency.

We consider variables, corresponding to the displacement of the nine cell’s centres in the and direction and we collect approximately observations available from DESY measurements Gall_2013aa in the matrix (see Eq. 20). The deformations are normally distributed, but they also present cross correlations. To obtain uncorrelated variables for the uq process, and to reduce the size of the parameter space, the truncated Karhunen–Loève decomposition introduced in section 3.1 is applied. This allows for the reduction from to variables according to the relative information criterion

 ∑Nt=7i=1σi∑N=18i=1σi≥0.95. (42)

The collocation points for the numerical quadrature are taken from a Smolyak grid of level , with Gaussian points along each direction, for a total of points Nobile_2008aa . For a given collocation point it is possible to use Karhunen–Loève expansion and the curves and to build the deformation curve and thus the deformed geometry . In Fig. 13 the deformed domain is shown.

The measured eccentric deformations, e.g. shown in Fig. 11, cause frequency shifts in the fundamental passband below . These shifts are rather harmless and can be corrected during the initial compulsory tuning procedure Padamsee_2008aa . For the following numerical study, the magnitude of the deformations is artificially increased by a factor of 50 to demonstrate that our IGA workflow can cope with arbitrary deformations while linearization or ad-hoc mesh deformation will fail. The increased deformations lead to a frequency shift for the fundamental passband in the order of , a level of accuracy that can be easily achieved by the isogeometric method by choosing second degree basis functions and approximately degrees of freedom.

The results for the standard deviations of the modes in the first TM passband are depicted in Fig. 12. It is noticeable, that the last two modes in the passband are more than times more sensitive than the other ones. It shall be noted that, in practice, the frequency and the field distribution of the -th TM mode of all considered cavities are mechanically tuned on an automatic cavity tuning machine before installation. The impact of this process which should reduce the variability of the last modes in the first monopole passband has not been considered in this work and could be subject to future research.

## 6 Conclusion

In this work we investigate the application of uq methods for the evaluation of eigenvalue sensitivity. In particular we are interested in the eigenfrequencies of rf cavities since they are known to be highly sensitive to geometry deformations.

Given the paramount role of geometry parameterization to achieve good accuracy, we discretize Maxwell’s time-harmonic eigenvalue problem with iga. The properties of the isogeometric basis functions also allow for a straightforward description of the domain deformation.

Finally, given that the ordering of the eigenfrequencies can change for different geometry configurations, we apply an eigenvalue tracking method, that exploits the properties of the iga matrices, to ensure a correct classification of the modes in different collocation points. However, the tracking has the positive side effect of reducing the computational costs for solving the repeated eigenvalue problems.

## Acknowledgment

The work of J. Corno, N. Georg and S. Schöps is supported by the Excellence Initiative of the German Federal and State Governments and the Graduate School of Computational Engineering at Technische Universität Darmstadt. J. Corno’s and N. Georg’s work are also partially funded by the DFG grants SCHO1562/3-1 and RO4937/1-1, respectively.

The authors would like to thank M. Kellermann and A. Krimm for their contributions to the first implementation and DESY for providing access to the cavity database.

## References

• (1) K. Halbach, R. F. Holsinger, Superfish – a computer program for evaluation of RF cavities with cylindrical symmetry, Particle Accelerators 7 (1976) 213–222.
• (2) U. van Rienen, T. Weiland, Triangular discretization method for the evaluation of RF-fields in cylindrically symmetric cavities, IEEE Transactions on Magnetics 21 (6) (1985) 2317–2320.
• (3) T. Weiland, On the unique numerical solution of maxwellian eigenvalue problems in three dimensions, Particle Accelerators 17 (DESY-84-111) (1984) 227–242.
• (4) T. Weiland, Computer modelling of two-and three-dimensional cavities, IEEE Transactions on Nuclear Science 32 (5) (1985) 2738–2742.
• (5) T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (2005) 4135–4195.
• (6) J. Corno, C. de Falco, H. De Gersem, S. Schöps, Isogeometric simulation of Lorentz detuning in superconducting accelerator cavities, Computer Physics Communications 201 (2016) 1–7.
• (7) P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford University Press, Oxford, 2003.
• (8) J. C. Nedelec, Mixed finite elements in , Numer. Math. 35 (3) (1980) 315–341.
• (9) T. J. R. Hughes, J. A. Evans, A. Reali, Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems, Computer Methods in Applied Mechanics and Engineering 272 (0) (2014) 290–320.
• (10) M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, G. N. Wells, The FEniCS project version 1.5, Archive of Numerical Software 3 (100) (2015) 9–23.
• (11) D. Arnold, R. Falk, R. Winther, Finite element exterior calculus: from hodge theory to numerical stability, Bulletin of the American mathematical society 47 (2) (2010) 281–354.
• (12) 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.
• (13) M. Spink, D. Claxton, C. de Falco, R. Vázquez, The NURBS toolbox, http://octave.sourceforge.net/nurbs/index.html.
• (14) B. Aune, et al., Superconducting TESLA cavities, Physical Review Special Topics-Accelerators and Beams 3 (9) (2000) 092001.
• (15) S. Gorgizadeh, T. Flisgen, U. van Rienen, Eigenmode computation of cavities with perturbed geometry using matrix perturbation methods applied on generalized eigenvalue problems, Journal of Computational Physics 364 (2018) 347–364.
• (16) L. Xiao, C. Adolphsen, V. Akcelik, A. Kabel, K. Ko, L. Lee, Z. Li, C. Ng, Modeling imperfection effects on dipole modes in TESLA cavity, in: IEEE Particle Accelerator Conference (PAC) 2007, 2007, pp. 2454–2456.
• (17) J. Heller, T. Flisgen, C. Schmidt, U. van Rienen, Quantification of geometric uncertainties in single cell cavities for bessy vsr using polynomial chaos, in: C. Petit-Jean-Genaz, G. Arduini, P. Michel, V. R. Schaa (Eds.), IPAC2014: Proceedings of the 5th International Particle Accelerator Conference, 2014, pp. 415–417.
• (18) D. Boffi, Finite element approximation of eigenvalue problems, Acta Numerica 19 (2010) 1–120.
• (19) A. Buffa, G. Sangalli, R. Vázquez, Isogeometric analysis in electromagnetics: B-splines approximation, Computer Methods in Applied Mechanics and Engineering 199 (2010) 1143–1152. .
• (20) C. de Boor, A Practical Guide to Splines, Vol. 27 of Applied Mathematical Sciences, Springer, New York, 2001.
• (21) L. Piegl, W. Tiller, The NURBS Book, 2nd Edition, Springer, 1997.
• (22) J. Corno, C. de Falco, H. De Gersem, S. Schöps, Isogeometric analysis simulation of TESLA cavities under uncertainty, in: R. D. Graglia (Ed.), Proceedings of the International Conference on Electromagnetics in Advanced Applications (ICEAA) 2015, IEEE, 2015, pp. 1508–1511.
• (23) 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.
• (24) J. Deryckere, T. Roggen, B. Masschaele, H. De Gersem, Stochastic response surface method for studying microphoning and Lorentz detuning of accelerators cavities, in: ICAP2012: Proceedings of the 11th International Computational Accelerator Physics Conference, 2012, pp. 158–160.
• (25) C. Schmidt, T. Flisgen, J. Heller, U. van Rienen, Comparison of techniques for uncertainty quantification of superconducting radio frequency cavities, in: R. D. Graglia (Ed.), Proceedings of the International Conference on Electromagnetics in Advanced Applications (ICEAA) 2014, IEEE, 2014, pp. 117–120.
• (26) D. Xiu, Numerical Methods for Stochastic Computations: A Spectral Method Approach, Princeton University Press, 2010.
• (27) M. D’Elia, M. Gunzburger, Coarse-grid sampling interpolatory methods for approximating gaussian random fields, SIAM/ASA Journal on Uncertainty Quantification 1 (1) (2013) 270–296.
• (28)

I. Babuška, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM review 52 (2) (2010) 317–355.

• (29) F. Nobile, R. Tempone, C. G. Webster, A sparse grid stochastic collocation method for partial differential equations with random input data, SIAM Journal on Numerical Analysis 46 (5) (2008) 2309–2345.
• (30) H.-J. Bungartz, M. Griebel, Sparse grids, Acta numerica 13 (2004) 147–269.
• (31) S. Smolyak, Quadrature and interpolation formulas for tensor products of certain classes of functions, in: Soviet Math. Dokl., Vol. 4, 1963, pp. 240–243.
• (32) E. Novak, K. Ritter, Simple cubature formulas with high polynomial exactness, Constructive approximation 15 (4) (1999) 499–522.
• (33) S. H. Lui, H. B. Keller, T. W. C. Kwok, Homotopy method for the large, sparse, real nonsymmetric eigenvalue problem, SIAM Journal on Matrix Analysis and Applications 18 (2) (1997) 312–333.
• (34) I. U. Ojalvo, Gradients for large structural models with repeated frequencies, in: Aerospace Technology Conference and Exposition, Society of Automotive Engineers, Warrendale, PA, 1986.
• (35) R. L. Dailey, Eigenvector derivatives with repeated eigenvalues, AIAA journal 27 (4) (1989) 486–491.
• (36) P. Deuflhard, Newton methods for nonlinear problems: affine invariance and adaptive algorithms, Springer, Berlin, 2004.
• (37) P. Jorkowski, R. Schuhmann, Mode tracking for parametrized eigenvalue problems in computational electromagnetics, in: 2018 International Applied Computational Electromagnetics Society (ACES) Symposium, 2018.
URL http://www.aces-society.org/conference/Denver_2018/
• (38) D. Yang, V. Ajjarapu, Critical eigenvalues tracing for power system analysis via continuation of invariant subspaces and projected arnoldi method, IEEE Transactions on Power Systems 22 (1) (2007) 324–332.
• (39) J. D. Jackson, Classical Electrodynamics, 3rd Edition, Wiley and Sons, New York, 1998.
• (40) K. Brackebusch, T. Galek, U. van Rienen, Automated mode recognition algorithm for accelerating cavities (2014).
• (41) R. B. Lehoucq, D. C. Sorensen, Deflation techniques for an implicitly restarted arnoldi iteration, SIAM Journal on Matrix Analysis and Applications 17 (4) (1996) 789–821.
• (42) G. L. G. Sleijpen, H. A. van der Vorst, A Jacobi-Davidson iteration method for linear eigenvalue problems, SIAM Journal on Matrix Analysis and Applications 17 (2) (1996) 401–425.
• (43) D. R. Fokkema, G. L. Sleijpen, H. A. Van der Vorst, Jacobi–davidson style qr and qz algorithms for the reduction of matrix pencils, SIAM journal on scientific computing 20 (1) (1998) 94–125.
• (44) R. B. Lehoucq, D. C. Sorensen, C. Yang, Arpack User’s Guide: Solution of Large-Scale Eigenvalue Problems With Implicityly Restorted Arnoldi Methods (Software, Environments, Tools), Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 1998.
• (45) D. Edwards, et al., Tesla test facility linac-design report, DESY Print March (1995) 95–01.
• (46) P. Gall, V. Gubarev, S. Yasar, A. Sulimov, J. Iversen, A database for the european XFEL, in: Proc. 16th Int. Conf. on RF Superconductivity, SRF (Paris, France, 23–27 September), 2013, pp. 205–7.
• (47) H. Padamsee, J. Knobloch, T. Hays, RF Superconductivity for Accelerators, Wiley, 2008.