Constrained multi-objective shape optimization of superconducting RF cavities considering robustness against geometric perturbations

High current storage rings, such as the Z-pole operating mode of the FCC-ee, require accelerating cavities that are optimized with respect to both the fundamental mode and the higher order modes. Furthermore, the cavity shape needs to be robust against geometric perturbations which could, for example, arise from manufacturing inaccuracies or harsh operating conditions at cryogenic temperatures. This leads to a constrained multi-objective shape optimization problem which is computationally expensive even for axisymmetric cavity shapes. In order to decrease the computation cost, a global sensitivity analysis is performed and its results are used to reduce the search space and redefine the objective functions. A massively parallel implementation of an evolutionary algorithm, combined with a fast axisymmetric Maxwell eigensolver and a frequency-tuning method is used to find an approximation of the Pareto front. The computed Pareto front approximation and a cavity shape with desired properties are shown. Further, the approach is generalized and applied to another type of cavity.

High current storage rings, such as the Z-pole operating mode of the FCC-ee, require accelerating cavities that are optimized with respect to both the fundamental mode and the higher order modes. Furthermore, the cavity shape needs to be robust against geometric perturbations which could, for example, arise from manufacturing inaccuracies or harsh operating conditions at cryogenic temperatures. This leads to a constrained multi-objective shape optimization problem which is computationally expensive even for axisymmetric cavity shapes. In order to decrease the computation cost, a global sensitivity analysis is performed and its results are used to reduce the search space and redefine the objective functions. A massively parallel implementation of an evolutionary algorithm, combined with a fast axisymmetric Maxwell eigensolver and a frequency-tuning method is used to find an approximation of the Pareto front. The computed Pareto front approximation and a cavity shape with desired properties are shown. Further, the approach is generalized and applied to another type of cavity.

I. INTRODUCTION
Accelerating cavities are metallic chambers with a resonating electromagnetic field that are used to impart energy to charged particles in many particle accelerators. Over the last few decades, a lot of research has been carried out on improving the shape and material properties of accelerating cavities. Many attractive features of superconducting radio frequency (RF) cavities, such as high intrinsic quality factor, have made them a favorable choice in applications where high continuous wave voltage is required [1]. The shape of the cavity, on the other hand, determines many figures of merit, such as the normalized peak electric and magnetic field on the surface of the cavity (E pk /E acc and B pk /E acc , respectively), geometric shunt impedance (R/Q), cell-to-cell coupling, etc. Therefore, based on the requirements of each specific accelerator, the shape of the cavity needs to be carefully optimized.
In many accelerators it is desirable to achieve a high accelerating field (E acc ) in order to improve efficiency and save on equipment cost. In superconducting RF cavities the maximum achievable E acc is limited by the maximum electric and magnetic field on the surface of the cavity. Thus, many cavity optimization methods have focused on minimizing E pk /E acc and B pk /E acc in order to provide room for increasing E acc [2][3][4][5][6]. In addition to the maximum achievable E acc , the surface losses of the cavity should be minimized, which can be achieved by maximizing G · R/Q (where G is the geometry factor). It was shown in [2] that high G · R/Q typically goes along with low B pk /E acc , so either of the two can be considered in the optimization. Since various conflicting objective functions are usually optimized simultaneously, a trade-off between these objective functions has to be considered. In [2,3] B pk /E acc is minimized, or G · R/Q maximized, while other figures of merit such as E pk /E acc , the wall slope angle, and the aperture radius of the cavity are kept fixed. In [4] a Pareto front between E pk /E acc and B pk /E acc was obtained by analyzing more than a thousand geometries. A geometry was then selected from this Pareto front based on the machine requirements (high gradient applications favor a smaller E pk /E acc , while low loss applications favor a smaller B pk /E acc ).
The optimization methods proposed in the literature mainly focus on optimizing the inner cells of a multicell cavity. The end half-cells are then optimized to get high field flatness in the cavity and ease the higher order mode (HOM) damping [2]. The optimization of singlecell cavities is slightly different than that of multi-cell cavities. The length of the inner cells of a multi-cell cavity is usually fixed to βλ/2, where λ is the wavelength of the fundamental mode (FM), and β is the ratio of the particle velocity to the speed of light (β ≈ 1 for particle velocities close to the speed of light). This restriction improves acceleration, since, as a particle traverses a cell, the direction of the electric field changes and the particle receives a force in the same direction. In single-cell cavities, the particle passes through only one cell, so there is no such restriction on the length of the cavity. Consequently, in single-cell optimization there is one more degree of freedom. Furthermore, since part of the field leaks into the beam pipe, in single-cell optimization the cell and the beam pipe have to be simulated together. This is in contrast to the optimization of the inner cells of multi-cell cavities, where only one cell, with appropri-arXiv:1905.13693v1 [physics.acc-ph] 31 May 2019 ii ate boundary conditions, is considered.
Many future accelerators, such as the Future Circular Collider (FCC), aim at colliding beams with unprecedented luminosities [7]. The FCC design study includes three colliders: a hadron collider (FCC-hh), a lepton collider (FCC-ee) and a lepton-hadron collider. The aim of FCC-ee is to study properties of Z, W and Higgs bosons, as well as top quark, with collision energies ranging from 90 GeV to 365 GeV. The high luminosity requirements demand a significant increase in the beam current of the machine, which in turn increases the HOM power deposited into the cavities by the traversing beam. Preliminary studies have suggested using single-cell cavities for FCC-hh and the Z-pole operating mode of FCC-ee (FCC-ee-Z) [8,9]. In the case of FCC-ee-Z, the main reasons for selecting a single-cell cavity are the issues related to the beam instability and high HOM power. In the FCC-ee-Z, the operating E acc is a few MV/m [8][9][10]. Reaching higher E acc is precluded due to limitations on the fundamental power coupler in providing high input power per cavity (e.g., for FCC-ee-Z this is around 1 MW per cavity). In such low E acc and high-current operations, minimizing surface peak fields should not be the primary goal. Instead, other figures of merit, in particular the HOM aspects, have to be taken into account right from the early design stages of the cavity. Enlarging the beam pipe radius is a common approach to untrap many HOMs and reduce the loss factor. However, even with an enlarged beam pipe radius, the first dipole band usually stays trapped and cannot propagate out of the cavity [11]. Additionally, the modes in the first dipole band usually have a large transverse impedance, which can cause transverse beam instability. Therefore, and in order to simplify their damping via HOM couplers, special attention should be given to the modes in the first dipole band.
In addition to the RF properties of the cavity, the cavity shape needs to be robust against geometric perturbations [12][13][14][15], which could, for example, arise from harsh operating conditions at cryogenic temperatures or manufacturing inaccuracies. Most importantly, the frequency of the FM has to remain at its nominal value in order to avoid excessive power being fed into the cavity for maintaining E acc , and the frequencies of the modes in the first dipole band should not be very sensitive to geometric changes in order to avoid unwantedly hitting a beam spectral line or harming their damping when coaxial HOM couplers are used.
The major contribution of this paper lies in finding a robust cavity shape while several properties of the FM and the first dipole band are optimized at the same time. For this purpose an optimization method for constrained multi-objective shape optimization of superconducting RF cavities is proposed. The focus of optimization is on axisymmetric single-cell cavities used in high-current accelerators, e.g., FCC-ee-Z, which unlike the multi-cell cavities have not been extensively studied before. As a part of the proposed optimization approach, a global sensitivity analysis is carried out in order to quantify the relative influence of each of the geometric parameters on the figures of merit of the cavity. The novel approach makes use of the most influential parameters, reduces the search space by discarding non-influential parameters, and provides valuable insights into the optimization of RF cavities.
The structure of the paper is as follows. The parameterization of the cavity cross section and the quantities of interest are described in section II. The results of the sensitivity analysis are shown in section III, and used for search space reduction in section IV. The constrained multi-objective optimization problem is defined in section V, the algorithm described in section VI, and the results presented in section VII. The approach described in sections II-VI is generalized to a different type of cavity in section VIII. Conclusions are drawn in section IX.

II. QUANTITIES OF INTEREST
Elliptical cavities are a commonly used cavity shape for the acceleration of particles with β ≈ 1. As shown in Fig. 1, the cross section of the cell of such cavities can be parameterized by seven variables, i.e., geometric parameters: R eq , L, A, B, a, b, and R i . The wall slope angle α is then determined by these geometric parameters, and the variable L pb determines the length of the beam pipe. In constrained multi-objective shape optimization the aim is to find the geometric parameters that best satisfy the given objectives and constraints. The objectives and constraints that need to be satisfied are determined by the requirements of the accelerator. The first constraint is to tune the frequency of the FM (f ), which is typically the TM 010 mode, to a desired value. Focusing on FCCee-Z, in this paper f is tuned to 400.79 MHz, which is the proposed value for both FCC-ee-Z and FCC-hh. In order to have a sufficient distance for decaying the FM leaked into the beam pipe, the beam pipe length L bp is set to the value of the wave length, i.e., to λ = 748 mm. The second constraint is α ≥ 90 • in order to avoid reentrant cavity shapes (due to the problems associated iii with chemical treatment).
In comparison with the middle-cells of multi-cell cavities, the electric field concentration around the iris region of single-cell cavities is lower because part of the electric field leaks into the beam pipe. Due to this reason and relatively low E acc of FCC-ee-Z, the surface peak fields are not the primary concern in the optimization. Lowering the surface losses helps to decrease the power released into the helium bath, which consequently reduces the amount of power required to maintain the cryogenic temperature. The surface losses of the FM are given by where R s is the surface resistance of the cavity, which depends on its material properties, V ||(r=0) the longitudinal voltage calculated along the longitudinal axis, and G and R/Q the geometry factor and the geometric shunt impedance of the FM, respectively. Both G and R/Q depend on the shape of the cavity [5]. Eq. (1) indicates that P c can be minimized by maximizing G · R/Q (the corresponding objective function is labeled as F 4 in Eq. (3)).
In addition to the properties of the FM, the properties of the HOM spectrum also have to be taken into account. The longitudinal loss factor is inversely proportional to the aperture radius R i of the cavity [16, p. 372]. Therefore, increasing the aperture radius lowers the longitudinal loss factor and consequently the HOM power deposited into the cavity by the beam. Larger aperture radius also helps to untrap the dangerous higher order monopole modes. As shown in Fig. 2, an aperture radius roughly above 145 mm untraps the monopole modes which typically have a large longitudinal geometric shunt impedance, i.e., the TM 011 and TM 020 modes. Therefore, a larger aperture radius is preferable (in Eq. (3), this is denoted by F 5 ). The first dipole band, however, stays trapped in the cavity, in particular the TE 111 mode whose frequency approaches the FM with enlarged R i .
The damping of the HOMs is usually done using coaxial couplers, waveguide (WG) couplers, beam pipe absorbers, or their combination. Coaxial and WG HOM couplers provide better damping of trapped modes as they can be placed close to the cell. In order to allow the propagation of the HOMs into the WG, the cutoff frequency of the first mode of the WG (the TE 01 mode) has to be between the FM and the first dipole mode. If the frequency of the first dipole mode is very close to f , a WG coupler with larger dimensions is required, which occupies more space in the cryomodule, in particular at 400.79 MHz. In that case the WG has to be wider in order to decrease the cutoff frequency of its TE 01 mode, and also longer in order to have a sufficient distance for decaying the FM leaked into the WG. Therefore, as another objective function, the distance between f and the frequency of the first dipole mode (f 1 ), which is typically the TE 111 mode, has to be maximized. Since f 1 is larger than f , this is equivalent to minimizing the nega-  [17]. An aperture radius Ri approximately above 145 mm helps to untrap the dangerous monopole HOMs. Note that the cavity in [17] is optimized such that the frequencies of the TE111 and TM110 modes are almost equal.
tive value f − f 1 (in Eq. (3), the corresponding objective function is labeled as F 1 ). Coaxial HOM couplers act like 3D resonant circuits that are optimized to have a notch at the FM and resonances (with a high transmission) at certain frequencies that require strong damping, such as the frequency of the TE 111 and TM 110 modes. A smaller difference between the frequencies of the first two dipole modes simplifies using a narrow band coaxial HOM coupler for damping them (such as the Hook-type coupler used in LHC cavities [18][19][20]). Therefore, another objective is to minimize the difference between the frequencies of the two trapped dipole modes (this is denoted by F 2 in Eq. (3)). Additionally, the sum of the transverse impedances of the first two dipole modes is minimized (F 3 in Eq. (3)). The following definition of the transverse impedance is used where k is the wave number, r 0 the offset from the axis, ω the angular frequency, and U the stored energy. The optimized cavity should be robust against geometric changes that might arise due to manufacturing inaccuracies or perturbations during operation. If the FM is detuned from its nominal value, additional power has to be fed into the cavity to maintain the same E acc . Thus, tuners are used to adjust the frequency to the desired value. Frequency tuning is usually carried out by changing the length of the cells via applying a longitudinal force to the cavity [1, p. 431]. A large tuning force should be avoided as it might plastically deform the cavity. Therefore, the cavities have to be fabricated with a high precision to avoid large tunings afterwards. On the other hand, a cavity fabricated with a higher precision is more expensive to fabricate and handle. As shown in Fig. 3, different cavity shapes can have different sensitivities with respect to geometric perturbations. E.g., cavities with a larger wall slope angle have a more robust f against changes in R eq (up to around 30 % difference between different shapes). FIG. 3. The local sensitivity of f with respect to Req plotted versus the wall slope angle α. Around 5'000 samples are studied. The frequency of each sample is tuned to 400.79 MHz. The local sensitivity is calculated using the forward difference method with the step h = 1 mm. Generally, a higher wall slope angle yields a more robust f against changes in Req.
The frequency of the trapped dipole modes could also vary due to perturbations or during the tuning of the frequency of the FM. Fig. 4 shows the local sensitivity of the frequency of the TM 110 mode with respect to R i against the geometric variable A. Sensitivity of the TM 110 mode against changes in R i can vary up to a factor of seven between different shapes. A change in the frequency of the dipole modes could harm their damping when coaxial HOM couplers are used. Coaxial HOM couplers are usually tuned to have a resonance at the frequency of the targeted HOMs and a mismatch between the frequency of the HOM and the coaxial HOM coupler could result in a poor damping of that mode. Therefore, in addition to the robustness of f , the robustness of the frequencies of the first two dipole modes is of importance. In Eq. (3), the objective functions corresponding to the minimization of the local sensitivities of f , f TE111 and f TM110 against geometric changes are denoted byF 6 ,F 7 andF 8 , respectively. To summarize, denoting by f 2 the frequency of the second dipole mode (which is typically the TM 110 mode) and the constrained multi-objective optimization problem considered in this paper can be written as In this formulation the local sensitivities of the FM and the dipole modes are calculated with respect to all seven geometric parameters, which entails a very expensive computation. In the following section a global sensitivity analysis is carried out in order to find the most influential geometric parameters on each of the frequencies. This information is then used to redefine objective functionsF 6 ,F 7 andF 8 , reduce the search space and, consequently, reduce the computational cost of the problem.

III. SENSITIVITY ANALYSIS
In order to determine which geometric parameters have the greatest influence on f , f TE111 and f TM110 , a variancebased global sensitivity analysis is performed [21,22]. The geometric parameters are considered to be independent, uniformly distributed random variables. The firstorder Sobol' indices and total Sobol' indices [23], representing the individual and total influences of these random variables on the variance of the quantities of inter-v est (QoIs), are computed using polynomial chaos (PC) expansion [22,[24][25][26][27]. The coefficients of the PC expansion are obtained non-intrusively, because non-intrusive methods allow the use of an existing solver as a black box, requiring only an evaluation of the QoIs in a set of either deterministic or random points. Since some of the deterministic points may correspond to infeasible cavity shapes, a random sample, i.e., a collection of random points, is used. According to [24], it is enough to use a sample of size where N is the number of geometric parameters and p the polynomial degree used. The values of Eq. (4), for a few relevant values of N and p, are given in Table I. For this analysis, the Uncertainty Quantification Toolkit (UQTk) [25,28] is used. The random sample is evaluated in parallel, taking into account only feasible cavity shapes with α ≥ 90 • . The changes in the sensitivity plots obtained for p = 3 (which needs 720 random points) are almost imperceptible, so only the case p = 2 (216 random points) is shown. The main sensitivities (i.e., first-order Sobol' indices), with p = 2 and with geometric parameters in the intervals shown in Tables II  and III, Table II. The polynomial degree is p = 2, so the size of the sample is 216 [cf. Table I].  Table III. The polynomial degree is p = 2, so the size of the sample is 216 [cf. Table I].
Sobol' indices are, by definition, normalized with respect to the total variance, so they (the first-order and higherorder indices) sum up to 1. Consequently, the sum of the first-order indices, representing the individual influences of the parameters, is at most 1. The higher-order indices represent mixed influences of the parameters, so the fact that the sum of the first-order indices, i.e., the height of the bars in the plots, is close to 1 indicates a low correlation between parameters. The geometric parameter R eq has the greatest influence on the frequency of the FM, f , (closely followed by A), so it is used to tune f to 400.79 MHz, i.e., to enforce the first constraint [cf. Eq. (3)]. This is explained in detail in the next section.

IV. SEARCH SPACE REDUCTION
According to Figs. 5 and 6, R eq has the greatest influence on f . Therefore, it is possible to use R eq to tune the frequency to 400.79 MHz. Specifically, for a point The QoIs are then computed for the cavity corresponding to d (the details are given in section VI B and Algorithm 2). The main sensitivities of the QoIs and R eq with respect to d 2 , . . . , d 7 , considering the wide (Table II) and narrow (Table III) intervals are shown in Figs. 7 and 8, respectively. The polynomial degree is again p = 2, but the number of variables is now 6, so only 140 training points are needed.
It can be seen in Figs. 5-8 that the influence of B on all of the QoIs is very low, and that b significantly influences only E pk /E acc . Therefore, since E pk /E acc is not part of any objective function in the optimization problem [cf. Eq. (3)], these two geometric parameters can be omitted. A natural way to do this is to set B = A and b = a, i.e., to consider circles instead of ellipses in the geometric parameterization of the cavity cross section [cf. Fig. 1].  need to be considered as well. However, due to the influence of A on R eq (Figs. 7 and 8) and their geometric connection (Fig. 1), in order to decrease the computation cost the local sensitivities with respect to A are omitted. Therefore, the constrained multi-objective optimization problem is the following [cf. Eq.
It is implied that, for each point I = (R i , L, A, a), the values of B and b will be set to B = A and b = a, and that f will be tuned to 400.79 MHz using R eq as described in Section IV.

A. Forward solver
In order to compute the values of the objective functions from Eq. (5) in a point d = (R eq , R i , L, A, B, a, b) [cf. Eq. (2)] time-harmonic Maxwell's equations with perfectly electrically conducting (PEC) boundary conditions (BC) are solved in the evacuated axisymmetric RF cavity parameterized by d. The mixed finite element method (FEM) leads to a generalized eigenvalue problem (GEVP) for each azimuthal mode number m ∈ N 0 [29,30]. For monopole and dipole modes, m = 0 and m = 1, respectively. If the cross section of the cavity vii is symmetric, as is the case for the single-cell elliptical cavity from Fig. 1, it is sufficient to solve time-harmonic Maxwell's equations for one half of it, once with PEC and once with perfectly magnetically conducting (PMC) BC on the cross-sectional symmetry plane (BC SP ).
To compute the properties of the TM 010 mode, the smallest eigenpair of the GEVP corresponding to m = 0 and PEC BC SP is found. This will be referred to as MAXWELL d (m,BC SP ) = MAXWELL d (0,PEC).
In order to compute the properties of the TM 110 and TE 111 mode, the smallest eigenpair of the GEVPs corresponding to m = 1 and PEC and PMC BC SP , respectively, is found. This will be referred to as MAXWELL d (1,PEC) and MAXWELL d (1,PMC).

B. Optimization algorithm
Since the minimizers of different objective functions are usually different points, the concept of dominance is used: a point I 1 = (R i,1 , L 1 , A 1 , a 1 ) dominates I 2 if it is not worse in any of the objectives, and it is strictly better in at least one objective. A point is called Pareto optimal if it is not dominated by any other point.
Because of conflicting objectives, the ability of an evolutionary algorithm (EA) to escape local optima, its suitability for parallelization, as well as good results in the area of particle accelerator physics [31][32][33][34] a multiobjective EA [35] is used to find an approximation of the set of Pareto optimal points, even though other methods, such as particle swarm optimization [36][37][38], ant colony optimization [39], simulated annealing [40], or artificial immune system [41] exist. The basic steps of an EA are shown in Algorithm 1.
In the case of the parameterization from Fig. 1, N = 7, the geometric parameters d 1 , . . . , d 7 are [cf. Eq. (2)] R eq , R i , L, A, B, a, and b, and j = 1, i.e., the frequency of the FM is tuned using d 1 = R eq . A design point for the CMOOP in Eq. (5), also called an individual in the context of an EA, is L, A, B, a, b).
The M individuals (M ∈ N) comprising the first generation are chosen randomly from the given intervals (Algorithm 1: line 1). In this paper, in the case of the singlecell elliptical cavity, these intervals are given in Table II. These M individuals are then evaluated (1:2), i.e., the corresponding objective function values are computed, as shown in Algorithm 2, in the following way.
First, an R eq [mm] ∈ [325, 375] [cf. Table II] is found such that the frequency of the FM of the corresponding cavity is f (d) = 400.79 MHz (2:1). This is done using the zero-finding method TOMS 748 [42]. Each evaluation of f requires finding the smallest eigenpair of the GEVP (which will be referred to as 'solving' the GEVP) corresponding to m = 0 and PEC BC SP . In case the cavity found this way is re-entrant, the second constraint is violated, so this individual is declared invalid and discarded from the population (2:2). Otherwise, the rest of the objective function values are computed (2:3-5).
In order to compute the properties of the TM 010 and TE 111 , two GEVPs need to be solved (2:3). In order to numerically compute the local sensitivity (i.e., the partial derivative) ∂f k /∂d l using the forward difference method, f k needs to be evaluated in the point d h,l = (d 1 , . . . , d l + h, . . . , d N ), i.e., another GEVP needs to be solved for the appropriate azimuthal mode number m f k and BC on the cross-sectional symmetry plane BC SP,f k . Once the first generation of the EA is evaluated, a predetermined number of cycles is performed, each resulting in a new generation (1:3-7). In each cycle, new individuals are created using crossover and mutation operators (1:4-5), and their objective function values are again computed (1:6) as shown in Algorithm 2. The new generation is then chosen to comprise approximately M fittest individuals (1:7).

A. Implementation and timings
The implementation of the optimization algorithm from the previous section is based on a combination of a massively parallel implementation of an EA [31,32], written in C++ and parallelized using MPI, with the axisymmetric Maxwell eigensolver [43]. For a point d, a mesh of half of the cross section of the corresponding cavity is created using the Gmsh [44] C++ API, and the resulting GEVPs are solved using the symmetric Jacobi-Davidson algorithm [45].
In (2:1) a cheap solve on a coarse mesh is performed first, in order to compute a good approximation for f , f coarse , (3-4 significant digits). Whenever f coarse is further away from 400.79 MHz than a given value ε (in this paper, ε = 1 MHz), the value f coarse is used in order to speed up the computation. Once the zero-finding method gets closer to 400.79 MHz than ε, a more expensive solve, on a much finer mesh, is performed in order to compute five significant digits of f . As a zero-finding method, the Boost 1 C++ library implementation of the TOMS 748 algorithm is used. Additionally, the zero-finding method is stopped as soon as five significant digits of f [MHz] match 400.79. When the search interval [325, 375] is first subdivided into three similarly-sized parts, (usually) using coarse solves, this entire approach requires, on average, only 2.2 fine solves in (2:1). The coarse solves use a mesh with around 20'000 triangles, and the fine ones around 500'000 triangles. On one core of the Intel Xeon Gold 6150 a coarse solve (creating a mesh, computing five smallest eigenpairs and the objective function values) takes around 2 s. A fine solve takes around 95 s (18 s for creating a mesh, 74 s for computing the eigenpairs, and 3 s for computing the objective function values). Additionally, in (2:3) two GEVPs on the same mesh need to be solved, so the mesh can be reused. Similarly, for each partial derivative that needs to be computed in (2:4), the solves are grouped in such a way as to avoid remeshing (eigenproblems corresponding to the same point, i.e., the same cavity shape, are solved consecutively).
To give an impression of the computation work for the entire optimization, using 108 processes on Euler IV (Euler cluster 2 of ETH Zurich, three Intel Xeon Gold 6150 nodes, each with 18 cores @ 2.7 GHz and cache size 24.75 MB) it took almost 15 h 19 min to compute 60 generations with M = 100 (around 30% of the evaluated individuals were discarded from the optimization).

B. Optimization results
The trade-off between the eight objective functions [cf. Eq. (5)] in the Pareto front approximation obtained in the 60-th generation of an optimization with M = 100 is illustrated in Fig. 9 using a scatter-plot matrix. The histogram at position (i, i) (i.e., on the main diagonal) shows the distribution of the objective function F i . The graph at position (i, j) where i < j (i.e., above the main diagonal) shows the values of F i (y axis) and F j (x axis) for the individuals in the last generation. The number in the i-th row and the j-th column where i > j (i.e., below the main diagonal) is the correlation coefficient between F i and F j . F 1 and F 2 are positively correlated because f 1 (in most cases f 1 is f TE111 ) is more sensitive to the geometric changes than f 2 , and the further it gets from f , the closer it gets to f 2 . Improving F 1 and F 2 leads to a higher transverse impedance (larger F 3 ) and a more sensitive f TE111 (larger F 7 ). The V-shaped curve of F 7 vs F 2 is created because the rise in F 2 (F 2 = |f 1 − f 2 |) corresponds to cases where f TE111 is higher than f TM110 . In order to simplify the damping of the dipole modes using coaxial HOM couplers, F 2 should be as close to zero as possible. A close-to-zero value of F 2 fixes F 7 to around 3.1 MHz/mm. However, in choosing a point from the Pareto front approximation, a low importance can be assigned to F 7 because the transverse impedance of the TE 111 mode is much smaller than that of the TM 110 mode.
All of the individuals, i.e., RF cavities, in the last generation are compared with four other storage ring singlecell cavities: CESR-B [46], HL-LHC [20], FCC HO18 [47] and FCC IC18 [17] in Table IV. Due to the conflicting nature of some of the objective functions, it is not feasible to surpass the other cavities in all objectives. However, all of these cavities are equally good in the Pareto sense.
TABLE IV. The value in i-th row, j-th column is the number of distinct cavities that have at least j objectives better than the cavity which corresponds to row i. The last (60-th) generation contains 95 distinct cavities. The design variable and objective function values, as well as QoIs, for a chosen cavity from the last generation, called FCC PR19 , are given in Table V. The shape of FCC PR19 is shown in Fig. 10. The chosen cavity outperforms the other cavities in five or six objective functions. Since the design of FCC-ee-Z demands the cavities to be operated at low E acc [8, p. 120], a small value of G · R/Q does not lead to a significantly large dynamical loss on the surface of cavity. Therefore, the value of F 4 is sacrificed in favor of other objective functions. The large iris radius of FCC PR19 ensures that the dangerous higher order monopole modes are untrapped and propagate out of the beam pipes. Compared to other shown cavities, the frequency of the FM of FCC PR19 is between 9% and 15% more robust against changes in R eq . The local sensitivity of FCC PR19 with respect to geometric parameters is given in Table VI. These numbers confirm that the assumptions made in section V are valid for FCC PR19 .  Table V.

VIII. GENERALIZATION
In this section the approach described in sections II-VI will be generalized and applied to a single-cell cavity with a symmetric cross section and half of its boundary defined as a complete cubic spline with horizontal end slopes. Such a parameterization is shown in Fig. 11 and referred to as 'single-cell spline cavity'. Related work includes the study of superconducting RF cavities whose boundary is a Bézier curve [48] or a non-uniform rational B-spline [49]. A summary of the approach described in sections II-VI and its application to the single-cell spline cavity is given in Fig. 12. Considering the intervals from Table VII, the main sensitivities are shown in Fig. 13, where it can be seen that the geometric parameter with the greatest influence on the frequency of the FM is y 4 . The main sensitivities in the case where, analogously to the use of R eq in section IV, y 4 is used to tune f to 400.79 MHz are shown in Sensitivity analysis (Fig. 13) Determine the most influential geometric parameter for f , dj (dj = y4) Sensitivity analysis when f is tuned using dj (Fig. 14) Determine the influential design variables for the chosen QoIs (f and fTM 110 ) and define the CMOOP (Eq. (6)) EA with constraint handling (Algorithms 1 and 2) Pareto front approximation and interesting individuals (Fig. 15) FIG. 12. A summary of Sections II-VI, applied to a different parameterization. The methods are red and the conclusions blue. The corresponding information for the single-cell spline cavity is referenced in parentheses. the elliptical cavity cannot be compared as their shapes are parameterized differently. Based on the information from Fig. 13 and Fig. 14 need to be taken into account, respectively. Due to a rather low value of R/Q ⊥,TE111 and the correlation observed between F 2 and F 7 in the previous section, we ignore the sensitivity of the TE 111 mode in the following optimization. Additionally, the influence of y 6 is small, so it can be omitted by, e.g., setting y 6 = 350 [cf. Ta , subject to f = 400.79 MHz. (6) It is implied that, for a point I = (L, y 1 , y 2 , y 3 , y 5 ), y 6 is 350 mm and f will be tuned to 400.79 MHz using y 4 .
The optimization problem is solved as described in sections VI B and VII A, and the shape of the chosen spline cavity shown in Fig.15. The dip in the contour could be removed by, e.g., adding an additional constraint to the optimization problem. That is, however, not the goal of this section.
The design variables, objective function values and QoIs for this cavity are given in Table VIII

IX. CONCLUSIONS
In this paper an algorithm for solving constrained multi-objective RF cavity shape optimization problems was proposed and applied to the problem of optimizing the shape of the superconducting RF cavity for the FCCee-Z. The shape of the cavity was optimized with respect to both the properties of the fundamental mode and the first dipole band, focusing in particular on robustness against geometric perturbations.
In order to decrease the computation cost, the results of a global sensitivity analysis were used to reduce the search space and define the objective functions of interest. A good single-cell elliptical cavity was found, and the algorithm generalized and applied to a different type of cavity. The proposed algorithm and its implementa-tion could also be applied to RF cavity shape optimization problems which take into account the properties of HOMs corresponding to arbitrary mode numbers.