Freeform shape optimization of a compact DC photo-electron gun using isogeometric analysis

Compact DC high-voltage photo-electron guns are able to meet the sophisticated demands of high-current applications such as energy recovery linacs. A main design parameter for such sources is the electric field strength, which depends on the electrode geometry and is limited by the field emission threshold of the electrode material. In order to minimize the maximum field strength for optimal gun operation, isogeometric analysis (IGA) can be used to exploit the axisymmetric geometry and describe its cross section by non-uniform rational B-splines, the control points of which are the parameters to be optimized. This computationally efficient method is capable of describing CAD-generated geometries using open source software (GeoPDEs, NLopt, Octave) and it can simplify the step from design to simulation. We will present the mathematical formulation, the software workflow, and the results of an IGA-based shape optimization for a planned high-voltage upgrade of the DC photogun teststand Photo-CATCH at TU Darmstadt. The software builds on a general framework for isogeometric analysis and allows for easy adaptations to other geometries or quantities of interest. Simulations assuming a bias voltage of -300 kV yielded maximum field gradients of 9.06 MV/m on the surface of an inverted insulator electrode and below 3 MV/m on the surface of the photocathode.


I. INTRODUCTION
Advanced applications of electron accelerators such as energy recovery linacs (ERLs) [1,2] require beams with high current and small emittance, therefore placing sophisticated demands on electron sources. State-of-theart DC high-voltage photo-electron guns are promising candidates for meeting these requirements [3,4]. The electrostatic design for this type of source, in light of optimizing the beam properties, has been discussed for many decades [5,6]. For example, there exists extensive research focused on the optimization of beam parameters depending on electron bunch parameters [7]. For instance, the electrode geometry was optimized for beam emittance in [8], using a set of parameters that describe a few key geometric features. In contrast, this paper is dedicated to optimizing the freeform shape of the electrode in terms of CAD basis functions to minimize the electric field strength, which has a crucial impact on field emission and thus still represents a major design problem depending on the specific geometry of the setup.
Low-level field emission can have a significant negative impact on the vacuum conditions within the gun * peter.foerster@tu-darmstadt.de and may severely degrade beam quality and operational lifetime [9]. High-level field emission can cause extensive damage to both electrode and insulator, necessitating repair or even replacement of the components. However, a high bias voltage is desired to provide sufficient initial acceleration for the beam and to minimize the emittance in spite of space charge effects. Since common negative bias voltages of DC photo-electron guns are in the range of −100 kV to −500 kV [10][11][12][13], the combination of such high voltages with a suitable electrode geometry and material poses a great challenge for the design of compact guns. The decisive limiting factor is the field emission threshold of the electrode material, imposing a maximum electric field strength on the geometric design. While increasing the curvature of the electrode surface reduces the field maximum, the overall size of the electrode is limited since the surface area susceptible to field emission should be kept small. Furthermore, a larger surface area also downgrades the vacuum conditions [14] and lastly, a larger electrode requires a larger vacuum chamber, which can be impractical due to cost and space constraints. A promising approach is the so-called inverted insulator geometry gun (IIGG) design [10,15], which significantly reduces the size of the electrode by placing the high-voltage insulator inside the vacuum chamber.
Practical experience shows that unavoidable material impurities and limitations in machining may cause sig- nificant variations in the field emission threshold. It is therefore paramount to keep the maximum electric field strength of the design well below the threshold. For stainless steel (1.4429 ESU), a commonly used electrode material, the threshold estimate from operational observations is 10 MV m −1 [16]. Other available materials, such as niobium, titanium, and molybdenum possess a higher threshold for field emission [17,18], but are more expensive and more difficult to machine. Common electrode designs range from simple spherical and cylindrical forms to more complex geometries like the T-shaped design used at JLab [15]. At TU Darmstadt a test facility for Photo-Cathode Activation, Test, and Cleaning using atomic-Hydrogen, Photo-CATCH, which is also dedicated to DC photo-electron gun research and development, has been established recently [19]. It uses an axisymmetric IIGG, featuring a two-part electrode consisting of a main electrode body and an extendable lift for photocathode loading [20]. An upgrade from −60 kV to −300 kV bias voltage has been envisioned and is currently under development. In order to meet design constraints concerning available space and chamber size, an adaptation and optimization of the existing geometry is necessary. The important components of the planned design are shown in Figure 1.
A key limitation of the design optimization process is the manual input and adaption of shapes based on simulations that must be repeated accordingly. An automation of these steps is desired in order to accelerate and simplify the design process. This leads to numerical shape optimization. Since the spatial description of the electric field inside the gun follows (the electrostatic approximation of) Maxwell's partial differential equations (PDEs), the shape optimization problem is PDEconstrained [21]. Furthermore, there commonly is no closed-form solution available for complex geometries, so the PDE is solved numerically, for example, by finite ele-ments [22]. PDE-constrained optimization is well known in the computational electromagnetics community, see the textbook [23] and references therein. Particularly in the context of electron guns, several design workflows to optimize their geometry have been proposed in the last decades [24][25][26][27][28]. However, all of them belong to the class of parameter-based optimization, i.e., the designer has to create a template which contains the design variables describing the geometry, e.g., width, height, and radius. This restricts the design space and is an inconvenient manual effort. On the other hand, computer aided design (CAD) tools allow freeform shapes in terms of splines and non-uniform rational basis splines (NURBS) [29,30].
Numerical shape optimization uses the parameters of these NURBS as the degrees of freedom (DOFs) and thus allows for an improved balance between design freedom and ease of implementation. Both parameter and shape optimization may also be used to describe the shape and position of holes in the geometry, however neither is able to introduce new ones. This requires a further generalization and leads to topology optimization; this however is not of interest for our application. An illustration of the different types of geometric design optimization is given in Figure 2.
There are additional important differences between our method and previous workflows. In the approaches cited above, each geometry realization is discretized separately, which requires rather fine spatial resolutions to avoid numerical errors due to remeshing ('mesh noise') and may again require additional manual intervention. According to Sandia Labs about 75 % of the simulation time in research laboratories is spent on modeling, parameterization, mesh generation, as well as pre-and post-processing [31]. Yet another distinction should be made regarding the quality of the resulting field solutions. The current state of the art in the accelerator community are low order finite element codes, e.g., POISSON [32]. However, even higher-order classical finite element codes, such as CST [33] yield noisy fields due to the lack of global regularity, see [34,Figure 4]. This is cumbersome for particle tracking and either needs smoothing or dedicated (symmetry preserving, mixed element) meshing. To avoid these problems, this paper proposes a spline-based shape optimization workflow using isogeometric analysis (IGA) w h Figure 2: Different types of design optimization. From left to right: parameter, shape, and topology optimization. [35], which integrates finite element analysis into the conventional NURBS-based CAD design workflow and allows for integrated particle tracking. IGA-based optimization is well established in many communities, but less explored in electromagnetism. However, [36,37] applied IGA-based optimization to accelerator magnets (without tracking) and more recently [38] suggested a freeform optimization workflow based on shape calculus for rotating electric machines; more references can be found in the survey article [39]. The paper is structured as follows: after this introduction, section II gives a short summary on CAD geometry handling and introduces splines. The following section III introduces the electric field problem, its weak formulation and discretization. Then section IV formulates the optimization problem and introduces numerical methods for its solution, and section V discusses the results for the particular gun in the context of Photo-CATCH. Finally, the paper closes with conclusions and an outlook.

II. SPLINES AND GEOMETRY
CAD models are essentially represented by B-splines [29] and NURBS [40], since they can exactly describe circular objects, allow local smoothness control, and give an intuitive definition of freeform curves and surfaces by so-called control points [41].
i=1 of a one-dimensional B-spline space S p α of degree p and regularity α may be constructed from a knot vector Ξ = (ξ 1 , ξ 2 , . . . , ξ n ) ∈ [0, 1] n , ξ 1 ≤ ξ 2 ≤ · · · ≤ ξ n using the Cox-de Boor algorithm [42] The knot vector uniquely determines the basis and its properties, including smoothness and the like. The knots need not be unique and the multiplicity m j of a knot value ξ j determines the continuity of the basis in that knot to be C p−mj . Furthermore, a knot vector is said to be open if its first and last knot each have multiplicity p + 1. For geometry modeling this usually is the case, since it makes the curve interpolatory in these knots. It also leads to a distinction between the first and last, and the internal knots. The latter influence the shape of the basis splines, as they represent the interfaces between each of the polynomial pieces (or elements) that make up the splines.

B. Geometry description
Given a set of control points {P i } N1 i=1 ⊂ R 3 , a threedimensional B-spline curve is described by a linear combination of the basis functions This representation is convenient for shape optimization for multiple reasons. For one, the uniqueness of the basis for a given knot vector leads to an interpretation of the control points as giving the curve its shape. As a consequence, changes in the coordinates of the control points directly translate to changes in the shape of the curve and most importantly they do so smoothly. An exemplary curve along with the corresponding basis is shown in Figure 3a. The extension of (1) to the bivariate case follows from choosing bases of S p1,p2 α1,α2 and a control mesh, given by an ordered set of N 1 × N 2 control points P i,j . A B-spline surface S P is then defined via and also volumetric (trivariate) mappings V P can be defined analogously. For the construction and handling of the bi-and trivariate geometry descriptions we make use of the free NURBS package [43].

C. Refinement
There are several approaches to refine an existing Bspline basis {B i,p } N1 i=1 . One is degree elevation, whereby the polynomial degree p of the basis functions is increased. To preserve the continuity of the original curve, the multiplicity of each knot ends up being increased alongside the degree. Furthermore each element, i.e., each polynomial piece, gains new control points equal to the increase in degree and the positions of all control points are recomputed such that the shape and parameterization of the curve are maintained. Figure 3b shows an example of the process.
A second refinement strategy is given by knot insertion. Here, an arbitrary internal knot is added to the knot vector. This does not impact the degree of the basis, however the total number of basis functions is still increased and the continuity of the basis is reduced in the new knot. For each inserted knot a new control point is added as well and again the positions of all control points are determined in a way to keep the shape of the curve intact. An illustration of the process is given in Figure 3c. Note that after inserting a knot at 0.5 the basis becomes C 0 continuous in that point, as expected. For a more comprehensive treatment of the geometry descriptions and refinement strategies we refer to [35].

III. FIELD FORMULATION AND DISCRETIZATION
Let Ω be the computational domain of the electron gun with boundary ∂Ω. In the absence of space charges, the electric field strength E = E(x), x ∈ Ω within the gun is described by the electrostatic subset of Maxwell's equations [44] ∇ × E = 0 and ∇ · (εE) = 0 in Ω, where the permittivity is given by Figure 4. Here ε 0 and ε ins are the permittivities of empty space and the insulator respectively. We assume that the domain is given by a multipatch spline mapping from the reference domainΩ = (0, 1) 3 to the physical domain, that is, Ω = Ω(P) in terms of control points P, see [45]. Introducing the electric scalar potential φ by E = −∇φ yields the boundary value problem [44, Section 1.7] along with the Dirichlet boundary conditions φ = φ Di on Γ Di (i = 0, 1, 2), see Figure 4.

A. Weak formulation
Exploiting the axisymmetry of the configuration we may restrict our analysis to Ω 2D , i.e., the ρ-z-plane. Let V = H 1 (Ω 2D ) denote the space of square-integrable functions with square-integrable gradients [22]. Following the Ritz-Galerkin approach, we deduce the weak form of (2) as: find φ ∈ V D such that and the approximated electric field strength follows from E h = −∇φ h . The linear system of equations reads where we only consider the N dof unconstrained coefficients as degrees of freedom, i.e., for 1 ≤ i, j ≤ N dof and include the coefficients known due to boundary conditions in the right-hand side . Please note, until now the basis functions have not yet been specified. The next section will propose to use Bsplines instead of the more common finite-element-type hat functions [35].

B. Isogeometric analysis
The main idea of IGA is to use B-splines or NURBS not only for the geometry description but also to represent the solution. This enables the solution of numerical problems, as the one defined above, on computational domains without introducing a geometric modeling error. Moreover, the thereby obtained geometry parameterization lends itself very nicely towards shape optimization, since it offers an intuitive set of degrees of freedom which immediately deform the underlying mesh. Finally, the use of high order B-spline basis functions in (5) guarantees rapid convergence and a high continuity of the solution [35].
→ Ω 2D denote the bivariate geometry mapping from the reference domainΩ 2D := (0, 1) 2 to the physical domain, from which eventually a 3D description is obtained by revolution [46]. Assuming that S P is piecewise smoothly invertible, we may define the approximation space V h using a gradient preserving transformation HereV h is a discrete space on the parametric domain, for which we elect to use the space of B-splines with degree p i and continuity α i along dimension i, denoted by S p1,p2 α1,α2 in the presented two-dimensional case. We have implemented the weak form and the discretization by IGA within the free software GeoPDEs [47].

IV. SHAPE OPTIMIZATION
The overall aim is to optimize the geometry of the electron gun to achieve two, possibly competing, goals. Firstly, we want to minimize the maximum electric field strength on the electrode surface, and secondly, we want to ensure a proper beam. Since the regions of the geometry that are relevant for field emission are far away from the cathode and the beam axis, we investigate the two problems separately. We optimize the shape first and then perform particle tracking to determine the beam properties of the optimized geometry.
Only the shape of the electrode is relevant for the geometry, i.e., the boundary Γ D1 in Figure 4. Furthermore, as can be seen in Figure 5a, it makes sense to restrict our attention to the domain Ω 2D opt , as indicated in Figure 4. The degrees of freedom for the optimization are given by the positions of the control points P of the curve C P (ρ, z) describing that part of Γ D1 , which intersects with Ω 2D opt . On a further note, the volume of the electrode may not exceed some fixed value V c due to space and weight considerations. For this, let V el (P) denote the volume of the electrode in dependence on Γ D1 (P), as characterized in Figure 1.
We allow geometries from an admissible set where ≤ is to be read componentwise. A accounts for constraints on the coordinates of the control points in terms of upper P i , and lower bounds P i , for example, to avoid intersections. The optimization problem is finally obtained as subject to where x = (ρ, z) is a position in the ρ-z-plane and the inner optimization, max E h (x; P) 2 , is approximated by a discrete maximum over a set of sample points which are used for the numerical quadrature of (6). The function f track denotes quantities of interest from the particle tracking, as defined in section V C, and tol describes associated bounds that ensure functionality. The above formulation of the optimization problem neglects two aspects which may also be critical for the electrostatic design of electron guns: the high-voltage cable conductor and the electric field magnitude at the triple point, where electrode, insulator and vacuum meet. The position of the triple point is highlighted in Figure 5 by the pop-outs. Including the high-voltage cable will have little effect on the maximum field on the electrode surface, but it may significantly influence the field at the triple point. To showcase the flexibility and effectiveness of our approach, we also add a term to (7) that aims to minimize the field at the triple point x tp leading to the objective function where w is a weighting factor to balance the two terms and U = U(x tp ) is a neighborhood of x tp , from which N U sample points are taken to approximate the value at the triple point. It should be noted that the field becomes infinite at the triple point due to the sharp corner of the geometry, see Figure 4 and [44, Section 2.11]. Therefore, we only evaluate the discrete representation of the field in the sample points, such that this issue can be mitigated. The results in section V indicate that this improves the design, but the results must still be treated with care. The question of which optimization algorithm to employ for solving the given problem is determined by the lack of smoothness of the min max problem, the unavailability of derivatives, and the nature of the constraints. In this work, we use a two step process consisting of the successive application of a global, followed by a local, optimization algorithm. The global algorithm (ISRES) is an evolution strategy based on a stochastic ranking to balance the objective function with a constraint based penalty function [48]. It does not need to compute or estimate derivatives of the objective function and at the same time is able to handle arbitrary nonlinear constraints, thus it meets our requirements. However the associated computational effort is comparatively high, as is to be expected with global optimization in general and evolutionary algorithms in particular. The local algorithm (COBYLA) works by creating linear approximations of both the objective and constraint functions via interpolating their evaluations at the vertices of a simplex [49]. This method again meets our criteria of not needing to compute derivatives of the objective function and being able to deal with nonlinear constraints. Even more importantly, the computational cost of this local algorithm is much less compared to that of the global one. We can therefore select a smaller tolerance on the change of the objective function over consecutive iterations and still obtain results in a shorter period of time.
For either algorithm we make use of the freely available implementations from the NLopt package [50]. Derivatives, alternative formulations, or approximations of the optimization problem (7) may allow for more sophisticated algorithms. In particular, one may look for a 'smoother' objective function that avoids the discrete maximum, or aim for convexity of the optimization problem.

V. NUMERICAL RESULTS
Based on the abstract formulation of the optimization problem given in (7), let us discuss the specific choices Original curve Initial guess Control points Figure 6: Original curve and the least squares fit serving as the initial shape for the optimization.
for the electron gun shown in Figure 1. We begin the optimization procedure with a B-spline curve of degree p = 7 without any internal knots. This is equivalent to a simple polynomial of degree 7, however both in terms of the presented isogeometric setting, and also for later refinements, it makes sense to interpret it in the B-spline context. The curve is a reasonable compromise between design freedom, simplicity, and the desire to obtain a smooth and manufacturable solution. The control points of the initial guess are determined by a least squares fit of the original 'flat' design, as shown in Figure 6, and the exact parameters can be found in [51]. In order to keep the overall geometry intact, the first and last control points are both fixed in their original positions over the course of the optimization.

A. Optimization results
Two successive optimization cycles are performed, as described in section IV. The first one uses a global optimization algorithm (ISRES) with a relative tolerance of 10 −3 on the objective function, and the second utilizes a local algorithm (COBYLA) with a relative tolerance of 10 −4 . The volume constraint is set at V c = 625 cm 3 , based on the assumptions that the insulator assembly can support a maximum weight of 5 kg and a stainless steel (type: 1.4404) electrode is used. The bounds for the admissible set can be found in [51].
The resulting shapes are shown in Figure 7 and for comparison, we also include the optimized curve obtained by using the modified objective function (9). As the highvoltage cable is still missing from the model, we only use the local algorithm for the second formulation, since the results are most likely not reliable for the final design. Nonetheless we find that the shapes look similar and the larger bulge at the back of the electrode for the modified objective is suitable for shielding the triple point. From this point onward we refer to the curve obtained via COBYLA, compare Figure 7, as the optimized shape if not explicitly stated otherwise and it will serve as the starting point for further analyses. Regarding the computational effort, the global optimization algorithm took about a week to find a solution satisfying our strict numerical tolerances, however it is possible to lower this number significantly by choosing a larger value for the tolerance or electing to only perform a local optimization. In contrast, the local algorithm only required computation times of around 7 hours to find a sufficiently accurate solution.
The electric field solutions corresponding to the original and both optimized curves are depicted in Figure 5. We observe a clearly visible reduction of the maximum field strength, and in addition, the change in the electric field magnitude along the electrode appears to be smoother for the optimized geometries. For the solution based on the isogeometric technique described in section III B, the open source package GeoPDEs is used [47]. The B-spline space is chosen as S 3, 3 2,2 and each of the parametric domains, of the patches indicated in Figure 4, is divided into n sub = 16 elements per coordinate direction using uniform knot insertion. For verification, both the original and optimized geometries are imported into CST Studio Suite 2019 and the field problem is solved using their adaptive mesh refinement with a tolerance of 10 −4 , based on a discretization with second order tetrahedral elements. The Dirichlet boundary conditions, as marked in Figure 4, are chosen as Γ D0 = 0 V, Γ D1 (P) = −300 kV, and Γ D2 = 1 kV.
The numerical values of the objective functions and the volume constraint, for the original and both optimized geometries respectively, are listed in Table I, where is introduced for brevity; refers to the used code. We observe a significant reduction in the maximum electric field strength, such that it falls well below the desired ISRES COBYLA Triple point 10 MV m −1 for the optimized electrode. The volume constraint is also fulfilled at 618 cm 3 even though the initial shape had violated the requirement (630 cm 3 ). Finally, it can be seen that the results of our code ('IGA') and CST's EM Studio ('CST') are in good agreement. Apart from the electrode, the maximum electric field strength on the cathode surface is also of interest. The numerical solution gives a value of 2.99 MV m −1 (E IGA c ) for the optimized geometry, well below the 3.9 MV m −1 that are documented for the former Jefferson Lab FEL gun that was routinely operated at −320 kV [52]. This value is expected to yield a sufficiently low energy spread [53], and preliminary results are shown in section V C. A more thorough integration of a particle tracking software into the shape optimization process will allow further improvement upon this value.
Another important quantity is the magnitude of the electric field at the so-called triple point, where electrode, insulator, and vacuum meet. Closeups of the field surrounding the triple point are shown in Figure 5. Our simulations predict a value of 3.27 MV m −1 (E IGA tp ) for the optimized geometry, a significant increase compared to the 2.55 MV m −1 for the original geometry. The studies in [54] suggest that one should aim for field strengths below 1 MV m −1 at the triple point. We see that the modified formulation (9) significantly reduces the field strength at this point. However, it may be necessary to further optimize the design and include additional shielding to minimize the field gradient at this critical point, as shown in [55]. Such measures could also influence the potential distribution and electric field magnitude along the outer insulator at the back of the electrode. Numerical results for these quantities can be seen in Figure 8 and show a nonlinear behavior. Adapting the design to linearize the field strength along the insulator surface may improve performance and reduce the chance of electrical breakdown at high voltages [12,54].
Lastly, the maximum value of the field gradient on the surface of the anode ring is 5.63 MV m −1 (E IGA ar ) according to our computations. It is also possible to further reduce this value, since the shape of the anode ring was not optimized in this work. The field magnitudes, at all critical points and for both the original and optimized geometries, are listed in Table I. Looking at the results, it can be seen that while decreasing the field strength at the triple point, the second formulation (9) leads to a higher gradient on the electrode surface.
We conclude our study of the geometry by looking at the convergence of the optimized parameters with respect to degree elevation and knot insertion. As discussed in section II, both refinement types add control points to an existing curve which increases the number of degrees of freedom. In the case of knot insertion the solution space is expanded even further, since the continuity of the basis in the new knot values is reduced, thus allowing a reduced continuity of the curve. The results in terms of the maximum electric field strength and the volume of the electrode are shown in Figure 9. In the case of degree elevation, the degree of the curve is continually increased by 1, i.e., p ∈ {7, 8, 9, 10}. For knot insertion, the inter- vals of the underlying knot vector are repeatedly halved by inserting additional knots, i.e., Ξ 0 = (0 1×7 , 1 1×7 ), , while the degree is kept constant at p = 7. The corresponding optimization cycles are carried out with COBYLA. One can clearly observe a correlation between the number of control points N opt and the quality of the solution. For this example, the solutions based on knot insertion seem to make better use of the available volume when compared to the ones from degree elevation, however this may simply be due to a local optimum.

B. Smoothness of IGA solutions
In section I and section III B we mentioned the higher global smoothness of the discrete fields when using IGA instead of classical finite elements. This is especially relevant for tracking applications, as the quality of the particle trajectories directly depends on the quality and properties of the electric field. Even if the tracking tool only supports pointwise data import; this is, for example, the case for ASTRA [56], regularity can be reconstructed if high order interpolation is used. ASTRA can make use of higher order polynomial interpolation internally, such that fields and first derivatives with respect to the space coordinates are continuous functions, e.g., when computing space charge effects [56,Sections 4.3,4.4 and 6.9]. To illustrate our point, Figure 10 shows a comparison between an IGA based solution using GeoPDEs and a linear finite element solution from CST. We use E x , E y to denote the x-or y-components of the electric field strength respectively and look at the field close to the beam axis. Again we employ the adaptive refinement tool of CST with a tolerance of 10 −4 to obtain a 'realistic' comparison. The results show a distinct difference in favor of IGA, since the unstructured mesh of the FE method, combined with the low order basis functions, leads to noisy fields.
As the shape optimization is only carried out in Ω 2D opt , one may assume that the tracking results are not severely affected by it, and therefore we did not consider the last constraint (8) within the optimization. Nonetheless, the following section shows a preliminary investigation of the emission process and electron acceleration to ensure a solid gun performance with the optimized geometry.

C. Particle tracking
Aside from fulfilling the previously discussed optimization criteria, the electric field should also be suitable for the initial acceleration of electrons emitted from the photocathode. In order to verify this, we perform simulations using the well-established particle tracking software AS-TRA [56]. Once the particle trajectories are computed, it is possible to evaluate statistical quantities that give insight into the gun performance. Two of these quantities of interest are the root mean square (RMS) beam width x rms ∈ R Nz and the related normalized transverse RMS emittance x ∈ R Nz , both in the x-and y-direction, where N z is the number of discrete points along the z-axis where the trajectories are known.
In this context, Figure 11 shows the initial macroparticle distribution in space for N p = 2 11 particles, which was obtained from a measurement using the DataRay BeamMap2 Beam Profiler. The laser spot has an oval shape with RMS radii of r x,rms = 0.41 mm and r y,rms = 0.72 mm. The corresponding data is also available at [51]. The emission times of the particles are drawn from a normal distribution with mean 0 s and standard deviation 5 ps. This represents a practical compromise between accuracy and simplicity, however there exists extensive work on the details of the emission process and the bunch profile in time [57]. The thermal emittance, representing the minimal emittance of a photoemission electron source, depends on material properties and the illuminating wavelength [58], and therefore on the photocathode material. In order to conduct a more general simulation, we thus assume the particles to have no initial momentum. The total bunch charge is estimated y (mm) Figure 11: Initial spatial distribution of the N p = 2 11 macroparticles to be emitted from the cathode. The positions are sampled from a measurement of the diode laser in use at Photo-CATCH.
at 100 fC, corresponding to the planned operation of the gun at Photo-CATCH. The gun is expected to produce a continuous waveform beam with a current of 300 µA at a repetition rate of 3 GHz, which is optimized for the operational parameters of the superconducting electron accelerator S-DALINAC at TU Darmstadt [59,60]. In addition to the already described parameters we initially choose a total of N p = 2 11 macroparticles, a time step of about ∆t = 0.244 ps for the Runge-Kutta integrator, and we set the grid for the electric field to be equidistant with n x = n y = 16 points (∆x = ∆y = 0.156 mm) in the transverse directions and n z = 256 points (∆z = 0.547 mm) in the longitudinal direction. The space charge computation is performed on a grid with n r = 64 (∆r = 0.039 mm) radial and n l = 64 (∆l = 2.188 mm) longitudinal cells. We definex rms and x to be reference solutions computed with refined parameters. More specifically, we look at two further simulations: The first uses half of the original time step ∆t, twice the number of grid points and cells, as well as double the number of macroparticles N p . The second one is the aforementioned reference, which again halves or doubles the parameter values. We then consider errors defined by in the computed statistical quantities with respect to the selected parameters. The corresponding numerical results are collected in Table II. We observe changes below 5 % between the reference and the first refinement step, indicating reliable results. The computational effort for the tracking simulations has two parts: The creation of fieldmaps based on the numerical field solution, and the actual particle tracking procedure. For the initial parameters, the fieldmap computation took about one hour and the tracking algorithm required around ten minutes to complete a single run. Thus the fieldmap computation clearly is the bottleneck in terms of integrating particle simulations into the optimization procedure. However, it may be possible to obtain sufficiently accurate maps with significantly less effort and only verify the results afterwards using a more accurate simulation. The computation times for the refined simulations can be estimated by a linear scaling, i.e., doubling the number of longitudinal grid cells or particles roughly doubles the execution time as well. It should also be noted that the creation of fieldmaps is a highly parallelizable task, such that considerable speedup could be achieved by a more optimized implementation.
The numerical solutions for x rms and y rms , interpreted as functions of z, are shown in Figure 12. Also included is the bunch length z rms , which may be defined analogously. It can be observed that the optimization slightly  worsens these beam parameters, but not by an intolerable amount. In contrast we even find a small improvement in the emittance values. Still, this indicates that future work could achieve yet better results by also considering the shape of the electrode leading into the photocathode. Our quantities of interest, evaluated at the chamber exit, are listed in Table III. They match the values measured at operating photo-electron guns [11,12], indicating the reliability of our simulations.
xrms yrms z ρ Figure 13: RMS beam widths in relation to part of the geometry, compare Figure 4. The values are scaled up by a factor of 10 to indicate that the focussing effect of the electric field is sufficiently strong.
Lastly, Figure 13 shows that the beam fits well inside the 20 mm aperture of the anode ring, even when considering an upscaling of the values by a factor of 10 to include a safety margin. In addition to the values collected in Table III, an RMS energy spread of ∆E rms = 48.4 eV is observed at the chamber exit for a field strength of 2.99 MV m −1 on the photocathode surface at −300 kV bias voltage. Comparing these values to the 84 eV with 2.5 MV m −1 at −200 kV bias voltage reported in [53], further supports the validity of our results. It should still be possible to improve on these values, since the laser shape that was used for the simulations was unprocessed after emission from the laser diode. Moreover, the procedure presented in this work focused solely on reducing the maximum electric field strength on the electrode surface and did not include the full shapes of the electrode or anode ring in the optimization. Taking their effects on the beam parameters into account could therefore provide additional improvements. As mentioned before, the initial momentum of the emitted particles needs to be included as well, which is expected to increase both RMS emittance and RMS energy spread. Specific optimization of the emission properties is the focus of ongoing work.

VI. CONCLUSION AND OUTLOOK
A successful IGA-based shape optimization of a DC high-voltage photo-electron gun was performed. The maximum electric field strength of the optimized geometry was computed as 9.06 MV m −1 , which is well below the field emission threshold of 10 MV m −1 . This constitutes a reduction in the maximum field gradient by more than 25 % compared to the initial design prior to optimization. Furthermore, the optimized electrode complies with the weight and volume restrictions, and the maximum field strength on the cathode surface is determined to be 2.99 MV m −1 , which allows for a sufficiently low energy spread of the electron beam. The procedure was carried out for an electrode voltage of −300 kV and an anode voltage of 1 kV, with a fixed anode cathode gap of 80 mm. Some beam parameters of the resulting geometry, namely the RMS beam widths and length, the normalized transverse and longitudinal RMS beam emittances, and the RMS energy spread were investigated using the particle tracking software ASTRA. Preliminary results were found to be in agreement with measurements at operating guns and exhibited values suitable for the setup at Photo-CATCH.
The work presented here focused solely on the optimization of the maximum electric field strength. In the future, the same procedure may be applied to the entire electrode, and the anode ring also, optimizing both their shapes and the anode-cathode distance. This includes fully coupling the shape optimization with a particle tracking software as well, to optimize the emission properties of the gun and investigate the influence of the electrode geometry on beam characteristics for ERL-typical bunch charges of 100 pC and above. In this context, it could also prove useful to allow for a direct evaluation of the spline basis functions within the tracking code, to make full use of the increased accuracy and smoothness without sacrificing computational efficiency.