Global Sensitivity Analysis on Numerical Solver Parameters of Particle-In-Cell Models in Particle Accelerator Systems

Every computer model depends on numerical input parameters that are chosen according to mostly conservative but rigorous numerical or empirical estimates. These parameters could for example be the step size for time integrators, a seed for pseudo-random number generators, a threshold or the number of grid points to discretize a computational domain. In case a numerical model is enhanced with new algorithms and modelling techniques the numerical influence on the quantities of interest, the running time as well as the accuracy is often initially unknown. Usually parameters are chosen on a trial-and-error basis neglecting the computational cost versus accuracy aspects. As a consequence the cost per simulation might be unnecessarily high which wastes computing resources. Hence, it is essential to identify the most critical numerical parameters and to analyze systematically their effect on the result in order to minimize the time-to-solution without losing significantly on accuracy. Relevant parameters are identified by global sensitivity studies where Sobol' indices are common measures. These sensitivities are obtained by surrogate models based on polynomial chaos expansion. In this paper, we first introduce the general methods for uncertainty quantification. We then demonstrate their use on numerical solver parameters to reduce the computational costs and discuss further model improvements based on the sensitivity analysis. The sensitivities are evaluated for neighbouring bunch simulations of the existing PSI Injector II and PSI Ring as well as the proposed Daedalus Injector cyclotron and simulations of the rf electron gun of the Argonne Wakefield Accelerator.


Introduction
Numerical models in scientific research disciplines are usually extremely complex and computationally intensive. A common feature to all is the dependency on numerical model parameters that do not represent an actual property of the underlying scientific problem. They are either prescribed in the source code and, thus, hidden to the user or can be chosen at runtime. Examples are the seed for pseudo random number generators, an error threshold in a convergence criterion, the number of grid points in mesh-based models or the step size in time integrators.
Latter two are often chosen to satisfy memory constraints or time limits. When applying new algorithms and modelling techniques the sensitivity of such input values on the response is usually studied by varying only a single parameter. While this captures the main influence of the tested parameter, possible correlations with other parameters are missed. A remedy is the evaluation of Sobol' indices [1] which are variance-based global sensitivity measures to express both individual and correlated parameter influences. Instead of Monte Carlo estimates, these quantifiers can easily be obtained by surrogate models based on polynomial chaos expansion (PCE). Uncertainty quantification (UQ) based on PCE is generally used in many areas of scientific computing and modelling [2,3,4,5] and several frameworks exist such as [6,7,8]. The coefficients of the truncated PCE that are required to determine Sobol' indices are mostly computed using the projection or regression method [9]. Since some numerical parameters are limited to integers, the former method is not applicable.
In this paper we study the sensitivity of adaptive mesh refinement (AMR) and multi bunch [10,11] parameters in Particle-In-Cell (PIC) simulations of high intensity cyclotrons where we use the new AMR capabilities of OPAL (Object Oriented Parallel Accelerator Library) [12] as presented in [13]. We further explore the sensitivity of a rf electron gun model with respect to the number of macro particles, the energy binning and the time step. The sensivitities are evaluated using ordinary least squares and Bayesian compressive sensing (BCS) [14,15]. Both methods are part of the uncertainty quantification toolkit UQTk [16,17]. The results are cross-checked using Chaospy [7] together with the orthogonal matching pursuit [18,19] regression model of the Python machine learning library scikit-learn [20,21].
Although we apply UQ on numerical solver parameters, it is a general method used for example in [5] to evaluate the sensitivities and to predict the quantities of interest due to uncertainties in physical parameters.
The paper is organised as follows: In Sec. 2.1 we elaborate the physical applications and their numerical modelling. An introduction to UQ is given in Sec. 3. The experimental setup and its results are presented in Sec. 4 and Sec. 5, respectively. In Sec. 6 are the final conclusions.

High Intensity Cyclotrons
Cyclotrons are circular machines that accelerate charged particles (e.g. protons) or ions (e.g. 2 ). Depending on the particle species and delivered energy, these machines find different applications ranging from isotope production [22,23] and neutron spallation [24] to cancer treatment [25,26]. An example that provides a beam for neutron spallation is the High Intensity Proton  [27,28] is the first acceleration stage delivering a 60 MeV/amu H + 2 beam to the Superconducting Ring Cyclotron (DSRC) with extraction energy of 800 MeV/amu. Since cyclotrons are isochronous, i.e. the magnetic field is increased radially in order to keep the orbital frequency constant, i.e. the revolution time per turn is energy-independent and, thus, the bunches lie radially on axis. A sketch showing five bunches denoted as circles on adjacent turns is depicted in Fig. 1a. As shown in [11], a small turn separation causes interactions between neighbouring bunches which yields to more halo particles (cf. Fig. 1b). In order to resolve these effects, the open source beam dynamics code OPAL [12] got recently enhanced with AMR capabilities [13] which adds more complexity to the numerical model. The influence of the AMR solver parameter settings on the statistical measures of the particle bunches is yet unknown and a too conservative AMR regrid frequency worsens the time-to-solution. Furthermore, the applied energy binnnig technique [10] to fulfill the electrostatic assumption (cf. Sec. 2.1.1) increases the computational costs considerably. Hence, the goal of this study is to quantify the impact of AMR solver and energy binning parameters in order to improve further computational investigations of these bunch interactions.
(a) Five bunches evolving radially on axis due to the isochronism of the cyclotron. The origin of the coordinate system denotes the center of the machine. The orbit radius r of each bunch is proportional to its energy E1 < E2 < · · · < E5.
(b) Separation of a particle bunch into core (blue) and halo (red) particles. In the PSI Ring the overall loss is on the order of 10 −4 which corresponds to a beam intensity of about 2 µA, i.e. all particles outside of approximately 3σ of a Gaussian distribution with standard deviation σ.

Numerical Model
The numerical model of neighbouring bunch simulations in OPAL, as presented in [11], is based on [10]. Due to the energy difference of the particle bunches on neighbouring turns a single transformation into the particle rest frame does not fully satisfy the requirements of the electrostatic assumption to solve Poisson's equation. Instead, each of the N macro particles is assigned to an energy bin b due to its momentum βγ according to where the binning parameter η is a measure of the energy spread. In each time step the force on a particle exerted by all others is the sum of the electric field contributions of each energy bin b evaluated in the appropriate rest frame of the particles obtained by a Lorentz transform with the proper relativistic factor γ b . The algorithm is summarised in Alg.

RF Electron Gun Model
To study the effect of energy binning we further use the example of the Argonne Wakefield Accelerator (AWA) [29,30,31] facility, an experiment setup for beam physics studies and accelerator technology developments. The facility is equipped with a photocathode rf electron gun that emits high intensity electron beams at high accelerating gradients ( 1 MV/m). Due to the high gradients the electrostatic approximation is invalidated and, hence, energy binning is necessary. In OPAL, we model the particle emission by p y = p total sin (ϕ) sin (θ) with x ∈ [0, 1] and θ ∈ [0, π] randomly sampled [12].

Quantities of interest
In accelerator physics interesting quantities of interest (QoI) also denoted as observables, of the co-moving frame are the rms beam size and the projected emittance that describes the phase space volume per dimension. The bracket · represents the moment. In order to quantify halo (cf. Fig. 1b), i.e. the tails of a particle distribution, we use two statistical definitions for bunched beams by [32,33], the spatial-profile parameter and the phase-space halo parameter with Eq. (4) and fourth order invariant (cf. also [34]) In case the bunch has uniform density Eq. (5) is zero due to the constant 15/7. An important quantity of interest in the rf electron gun model is the energy spread

Non-Intrusive Uncertainty Quantification
In contrast to intrusive UQ, non-intrusive UQ uses the computational model as a black box. In this paper we only give a short overview following the description and notation of [5,9,1,35,36].
A detailed introduction to UQ in general is found in e.g. [37,38].

Surrogate Model based on Polynomial Chaos Expansion
The PC-decomposition originates from [39], where a random variable of a Gaussian distribution is represented as a series of multivariate Hermite polynomials of increasing order. Later, this description was generalized by [40] to other probability measures and their corresponding orthogonal polynomials.
Let a multivariate polynomial Ψ α (ξ) of dimension d ∈ N >0 and multiindex α ∈ N d be defined with orthogonal univariate polynomials {ψ j } d j=1 . The response of a model m(x) with random input vector x ∈ Ω 1 × · · · × Ω d , where Ω j ∀j{1, . . . , d} denotes the sample space of the j-th random variable, can then be represented as PDF of ξ j Polynomial Basis Support Table 1: Examples of the Wiener-Askey polynomial chaos of random variables ξ j with appropriate probability density function (PDF) [40].
where the basis of a random input component is determined by its probability distribution The truncation scheme is not clearly defined. A common rule, which is also used here, is the so-called total order truncation that keeps all multiindices α for which ||α|| 1 ≤ p. This yields a number of multiindices. Three other schemes are explained in [35]. In the next sections we describe methods to compute the coefficients c αi of Eq. (11).

Projection Method
The (spectral) projection method computes the coefficients of Eq. (11) making use of the orthogonality of the basis functions, i.e. Ψ αi (ξ)Ψ αj (ξ) = 0 with ∀i = j. Thus, the PC coefficients are given by While the denominator is evaluated by analytic formulas (see examples in the appendix of [9]), the numerator is computed by Gaussian quadrature integration where integration points, i.e. high fidelity model m(x) evaluations, are required.

Linear Regression Method
The coefficients of Eq. (11) can also be computed with regression-based methodŝ with regularization parameters λ 1 , λ 2 ≥ 0. The minimization problem is called ordinary least squares if λ 1 = λ 2 = 0, elastic net [43] if λ 1 , λ 2 > 0, ridge regression [44] (or Tikhonov regularization) if only λ 1 > 0 and Lasso [45] if only λ 2 > 0. In matrix form the problem readŝ . . , c αp ∈ R P ×1 and system matrix A ∈ R N ×P . In case of λ 2 = 0, the coefficients of Eq. (16) are obtained in closed form bŷ with P × P identity matrix I. In contrast to the projection method (cf. Sec. 3.2), this method does not require a fixed number of samples N . However, [9] gives an empirical optimal training sample size of with P defined in Eq. (12).

Orthogonal Matching Pursuit
The matching pursuit (MP) is a greedy algorithm developed by [18] which was enhanced by [19] to obtain better convergence. This improved method is called orthogonal MP (OMP).
In terms of PCE, the algorithm searches a minimal set of non-zero coefficients to represent the model response, i.e.ĉ = arg min c ||y − Ac|| where N c is a user-defined maximum number of non-zero coefficients [46]. The vectors and matrices are defined according to Eq. (16). It is an iterative procedure where in the (i + 1)th step a new coefficient vector c i+1 is found that has one non-zero value more than the previous vector c i and that maximizes the inner product between the remaining residual r i = |y −ŷ i | and the model estimateŷ i+1 . In OMP, the new modelŷ i+1 in Eq. (20) is obtained by a linear combination of the previous modelŝ and together with the orthogonality constraint

Bayesian Compressive Sensing
As stated in [14,15], the linear regression model Eq. (15) can be interpreted in a Bayesian manner, i.e.
with posterior distribution p (c|D), likelihood p (D|c), prior p (c) and evidence p (D) of training . The likelihood is assigned a Gaussian noise model with variance σ 2 . It is a measure of how well the high fidelity model is represented by the surrogate model Eq. (11). In order to favour a sparse PCE solution, a Laplace prior is chosen. Using Eq. (25) in the maximum a posteriori (MAP) estimate for c, i.e.
arg max An iterative algorithm to obtain the coefficients is described in [15]. It requires a user-defined stopping threshold ε that basically controls the number of kept basis terms, with more being skipped the higher the value is. The overall method is known as Bayesian Compressive Sensing (BCS) [14,15].

Sensitivity Analysis
Sobol' indices [1] are good measures of sensitivity since they provide information about single and mixed parameter effects. In addition to these sensitivity measures, there are various other methodologies such as Morris screening [48]. A survey is presented in [49] on the example of a hydrological model. Instead of Monte Carlo, Sobol' indices are also easily obtained by surrogate models based on PCE as discussed in the following subsections.

Sobol' Sensitivity Indices
In [1] Sobol' proposed global sensitivity indices that are calculated on an analysis of variance with mean and for k = i 1 , . . . , i s and s = 1, . . . , d. Since Eq. (30) holds, the components of Eq. (28) are mutually orthogonal. Therefore, the total variance of Eq. (28) is that can also be written as where Based on Eq. (32) and Eq. (33), the Sobol' indices are defined as The first order indices S i are also known as main sensitivities. They describe the effect of a single input parameter on the model response. The total effect of the i-th design variable on the model response, proposed by [51], is the sum of all Sobol' indices that include the i-th index, i.e.

Sobol' Indices using Polynomial Chaos Expansion
Instead of Monte Carlo techniques, Sobol' indices can be estimated using surrogate models based on PCE since the truncated expansion can be rearranged like Eq. (28). The Sobol' estimates are then given by [ where and varianceD The main and total sensitivities are computed bŷ

Confidence Intervals using Bootstrap
In this subsection we briefly outline the computation of confidence intervals for the estimates of Sobol' indices using the bootstrap method [52]. In the context of PCE, the bootstrap method has already been applied in [53], where it is referred to as bootstrap-PCE (or bPCE and bootstrap sample mean

High Intensity Cyclotrons
In order to study the effect of AMR solver parameters and energy binning in neighbouring bunch simulations we perform sensitivity experiments with three different high intensity cyclotrons, the PSI Ring [55], the PSI Injector II [56] and the DAEδALUS Injector Cyclotron (DIC) [27,28]. We always accelerate 5 particle bunches with 10 6 particles.      L z = 32 grid points. The final number of emitted macro particles is given by

RF Electron Gun Model
where the particle multiplication factor p f is an integer. This parameter basically controls the number of particles per grid cell and, hence, the noise of the PIC model.

Surrogate Model Selection
In order to avoid overfitting we proceed like [35] where the truncation order of the PC expansion and the settings of the regression models are chosen such that the relative l 2 error between the surrogatem(x) and high fidelity model m(x) of the training and validation set is approximately of equal magnitude. As an additional error measure we also compare the relative l 1 error

Results
The estimated sensitivities are obtained from PC surrogate models where we use either ordinary least squares (OLS) and Bayesian compressive sensing (BCS) of UQTk [16,17]

High Intensity Cyclotrons
In order to study the effect of the input parameters we evaluate the sensitivities of the halo parameters Eq. (5) and Eq. (6) with respect to the center bunch of the 5 bunches (cf. Fig. 1a).
The initial kinetic energy of the center bunch in the different models is approximately 98 MeV, 25 MeV and 17 MeV for the PSI Ring, PSI Injector II and DIC, respectively. As shown in Fig. 4,   Fig. 6 correlates with the decrease of the standard deviation in Fig. 3. It is best observed for h x at around 215 • or H x between 195 • and 255 • . In contrast to the PSI Injector II, the DIC also strongly depends on the regrid frequency.
It has an average main sensitivity of approximately 60 % for h x . The parameters also exhibit more correlations as observed between the main and total sensitivities (cf. Tab. 8). The standard deviations for the DIC are one order of magnitude smaller than for the PSI Injector II, causing the confidence intervals to increase as illustrated in Fig. 10. This effect is even stronger in the PSI Ring where Coulomb's repulsion is less dominant and the halo parameters are smaller (cf. Fig. 11) compared to the PSI Injector II. The standard devation is in the order of 10 −4 which is why the confidence intervals in Fig. 14

RF Electron Gun Model
In order to approximate the high fidelity model we use PC surrogate models of second order where the BCS method uses a tolerance of ε = 10 −9 and the OMP method is stopped once 7 non-zero coefficients are found. In Fig. 16

Conclusions
In this paper we discussed uncertainty quantification based on polynomial chaos expansion and four numerical methods to compute the polynomial coefficients. Beside a cheap surrogate model that mimics the high fidelity model, it has the additional benefit to easily evaluate the Sobol' sensitivity indices. As demonstrated in this paper, this technique is not only suitable to gain knowledge about the sensitivity of physical parameters (e.g. initial beam properties) on the quantities of interst but also numerical parameters of computer codes. Since some tested numerical parameters might be limited to integers, the projection method to obtain the polynomial coefficients is, however, not applicable. Instead, regression-based methods, Bayesian Compressive Sensing or Orthogonal Matching Pursuit and others have to be applied. A further difficulty with numerical parameters is a fair random sampling. In some cases (cf. Fig. 2) a straightfor-ward, uniform sampling of the parameter yields biased input data and, hence, may induce wrong conclusions. To circumvent, we perform a stratified sampling that guarantees a well-balanced distribution.
The sensitivity studies of the three high intensity cyclotrons show that the sensitivity results are different among accelerators of the same type. While the AMR threshold is the most important parameter in the PSI Injector II with a sensitivity of about 90 %, the regrid frequency is relevant in the DAEδALUS Injectory Cyclotron (DIC), too. Large bootstrap confidence intervals of the Sobol' indices are due to small variances in the quantities of interest in the order of O(10 −4 ). In such a case no reliable statements based on the sensitivity estimates can be done.
In contrast to our intuition the standard deviation of the halo parameters remain pretty constant throughout one turn in the PSI Ring and PSI Injector II and the considered three turns in the DIC. Nevertheless, these findings give rise to computational savings. Without losing significantly on accuracy (cf. Fig. 3, Fig. 7 and Fig. 11), energy binning can be totally switched off for these cyclotrons. In addition, this reduces the amount of AMR hierarchy updates which reduces the time-to-solution even further since the operators of the adaptive multigrid solver to solve Poisson's equation do not need to be set up in every time step. To illustrate this, we take the benchmark example in [13] that solves Poisson's equation 100 times using a three level AMR hierarchy with a base level of 576 3 grid points. The benchmark running on 14 400 CPU (Central Processing Unit) cores shows that the matrix setup due to AMR regriding takes up 42.15 % computing time. A reduction of the regrid frequency by a factor 10 yields a speedup of 7.10 in the matrix setup timing. In our UQ samples, the speedup between the cheapest and most expensive model is at least 2.0 and at most 3.3.
Another interesting case we have studied is the rf electron gun model of the AWA. Relevant parameters for this model are the energy binning N E and the time step ∆t. The particle multiplication factor p f , that basically controls the number of particles per grid cell, is only important for the longitudinal beam size. Although N E and p f have together an average main sensitivity of 86 %, ∆t is the dominating parameter close to the cathode. An adaptive energy binning and time step scheme that is based on Sobol' sensitivity indices is therefore a possible future enhancement to reduce the time-to-solution for a target accuracy.