On Nonintrusive Uncertainty Quantification and Surrogate Model Construction in Particle Accelerator Modeling | SIAM/ASA Journal on Uncertainty Quantification | Vol. 7, No. 2 | Society for Industrial and Applied Mathematics

Using a cyclotron-based model problem, we demonstrate for the first time the applicability and usefulness of an uncertainty quantification (UQ) approach in order to construct surrogate models. The surrogate model quantities, for example, emittance, energy spread, or the halo parameter, can be used to construct a global sensitivity model along with error propagation and error analysis. The model problem is chosen such that it represents a template for general high-intensity particle accelerator modeling tasks. The usefulness and applicability of the presented UQ approach is then demonstrated on an ongoing research project, aiming at the design of a compact high-intensity cyclotron. The proposed UQ approach is based on polynomial chaos expansions and relies on a well-defined number of high-fidelity particle accelerator simulations. Important uncertainty sources are identified using Sobol’ indices within the global sensitivity analysis.

1. Introduction.Uncertainty quantification (UQ) describes the origin, propagation, and interplay of different sources of uncertainties in the analysis and behavioral prediction of generally complex and high-dimensional systems, such as particle accelerators.With uncertainty, one might question how accurately a mathematical model can describe the true physics and what impact the model uncertainty (structural or parametric) has on the outputs from the model.Given a mathematical model, we need to estimate the error.How accurately is a specified output approximated by a given numerical method?Can the error in the numerical solutions and the specified outputs be reliably estimated and controlled by adapting resources?For example, in beam dynamics simulations with space charge, grid sizes would be such a resource.
UQ techniques allow one to quantify output variability in the presence of uncertainty.These techniques can generally tackle all sources of uncertainties, including structural ones.However, in this paper we focus on parametric uncertainty of input parameters.The moments of the output distributions are sampled using Monte Carlo [1] or quasi-Monte Carlo [2] methods, or newer approaches such as Multilevel Monte Carlo [3].Other approaches exist and are known as nonsampling-based methods.For an introduction to response surface methods, see [4,5].The most popular method these days, which is used in this paper, is the polynomial chaos (PC)-based method [6].Strictly speaking, PC also requires sampling, but it is not random sampling as in Monte-Carlo-type approaches.
In probabilistic UQ approaches, one represents uncertain model parameters as random variables or processes.Among these methods, stochastic spectral methods [16,17] based on PC expansions [6,18] have received special attention due to their advantages over traditional UQ techniques.For a more detailed discussion on that subject, consult the introduction of Hadigol, Maute, and Doostan [14] or, alternatively, the book of Smith [19].
In the field of particle accelerator science, nonintrusive methods are far more attractive than intrusive methods.The complexity of the physics model would most likely require a total rewrite of the existing simulation packages in order to facilitate intrusive methods.Because nonintrusive methods allow the use of existing beam dynamics codes as black boxes, they are the methods of choice.A nonintrusive method to solve an inverse was proposed in [20].A proton beam from a linear charged particle accelerator is focused through the use of successive quadrupoles.The goal of the inverse problem is to find the unknown initial state of the beam, in terms of particle position and momentum.Measurement data on the projection of the phase space was used where available beyond the focusing region.This setup is that of an inverse problem, in which a computer simulator is used to link an initial state configuration to observable values, and then inference is performed for the distribution of the initial state.The Bayesian approach used allows estimation of uncertainty in the initial distributions and beam predictions.
In this paper, we use OPAL [21,22] as the black-box solver.As we will see later, only independent solution realizations are needed, and hence embarrassingly parallel implementation is straightforward.
The proposed PC approach, first introduced in [16,23], computes the statistics for the Quantity of Interest (QoI) with a small number of accelerator simulations.However, in contrast to [16,23] we do not exploit the sparsity of expansion coefficients; this is subject to further research.Additionally, the presented UQ framework enables one to perform a global sensitivity analysis (SA) to identify the most important uncertain parameters affecting the variability of the output quantities.
To avoid confusion, we first point out a misnomer by mentioning that polynomial chaos [6] and chaos theory [24] are unrelated areas.Originally proposed by Wiener [6] in 1938 (prior to the development of chaos theory-hence the unfortunate usage of the term chaos), PC expansions are a popular method for propagating uncertainty through low-dimensional systems with smooth dynamics.
This work presents a sampling-based PC approach to study the effects of uncertainty in various model parameters of accelerators.As a model problem, we use the central region of a "PSI Injector 2 like" high-intensity cyclotron, where we only consider the first 10 turns of the cyclotron.While this paper's focus is mainly to introduce UQ to the field of particle accelerator science, we add a realistic example of an ongoing design effort.
1.1.Motivation in lieu of an actual research project.Searches for CP violation in the neutrino sector and "sterile" neutrinos, respectively, need a lot of statistics, i.e., events.This translates, in the Decay-At-rest Experiment for δ CP violation At a Laboratory for Underground Science (DAEδALUS) [40] and the Isotope Decay-At-Rest experiment (IsoDAR) [25], into high fluxes of protons and compact accelerators in our example cyclotrons.The detailed exposition on how UQ and PC expansion (PCE) is used in an ongoing research project is given in section 5; here we want to motivate this approach, and fix language and notation.In high-intensity accelerators, the working point is to a large extent defined with setting the flux of particles per time, i.e., the intensity I (cf.section 4).In the compact accelerator studied in section 5, the angle θ of the spiral inflector also defines a working point.The spiral inflector and part of the central region (magnets and cavities not shown) are depicted in Figure 1.The spiral inflector can be rotated around the z axes by an angle θ.
We consider these as design or controllable parameters.The machine is operated at only a few distinct different values, e.g., in high-intensity, i.e., production, mode, or for machine development in lower-intensity mode to prevent damage to or activation of the accelerator.Similar arguments can be made for θ; cf.section 5.
The other category of parameters are the model parameters.The model parameters are quantities that are either not measurable or measurable with an associated uncertainty.In the problem of section 5, the radius r of the injected particles and the associated radial momenta p r (cf. Figure 1) cannot be measured in-situ, and hence empirical values or values from simple models together with a meaningful PDF are used.Other quantities can be measured, for example, the phase of the cavity φ, but we want to find optimal values.In this paper, all design parameters are independent, bounded, and uniformly distributed.
Maximizing performance in high-intensity accelerators has two main dimensions: 1. maximize the number of transmitted particles throughout the accelerator and at the same time 2. minimize particle losses.In Figure 1, one can already see by eye that particle tracks are terminating at the not shown walls.A few are marked as red circles for illustration purposes.The tolerable particle losses have to be at levels of 3 to 4 standard deviations of the particle density.Particle losses are associated with "halo," i.e., particles that are sufficiently far away from the core of the distribution, such that they have a high probability of being lost.This all translates into the necessity to solve large N-body problems, taking into account the nonlinear particle-particle interaction, together with complicated boundary conditions.Furthermore, as hinted above, with the design of such complex scientific instruments, large-scale multiobjective optimization must be worked out and correlations and sensitivities identified.This motivated the search for inexpensive to evaluate surrogate models and is one of the main motivations behind this work.
In section 2, we present our stochastic modeling approach, which is based on nonintrusive PC expansions.After the derivation of the surrogate model, we then continue with reviewing a global sensitivity analysis approach using Sobol' indices.Section 3 introduces the simulation model and a model problem.Section 4 applies the UQ to the stated problem and shows the main features of this approach.The features are general in nature and not restricted to cyclotrons.Section 5 reports on an ongoing design effort using UQ.Conclusions are presented in section 6. [6] introduced PCE.In 1991, Ghanem and Spanos [16] reintroduced this technique to the field of engineering.They first studied problems with Gaussian input uncertainties and extended their method to non-Gaussian random inputs.In their studies, orthogonal polynomials of the Askey scheme were used.This is known as a generalized polynomial chaos (gPC) expansion [23].The method of gPC expansion provides a framework to approximate the solution of a stochastic system by projecting it onto a basis of polynomials of the random inputs.

UQ via PCE. Wiener in 1938
An overview and some details on the correspondence between distributions and polynomials can be found in [26].A framework to generate polynomials for arbitrary distributions has been developed in [27].The advantage of using PC is that it provides exponential convergence for smooth models.However, the approach suffers from the curse of dimensionality, making them challenging for problems with number of parameters in the range 10 to 50.To mitigate the curse of dimensionality, sparse grid techniques have traditionally been used [28,29].More recently, iterative methods to propagate uncertainty in complex networks have also been developed [30,31,32].
2.1.The surrogate model.Suppose one is designing or optimizing complex systems such as particle accelerators.As a particular example, consider the case of a high-intensity hadron machine.In such a machine, one needs to characterize and minimize some QoIs (for example, halo), and at the same time increase the beam quality.In order to accomplish this task, usually a large number of design and model parameters, in the search space D (cf. Figure 2), have to be considered.Let us furthermore assume that D is the admissible space, i.e., where the accelerator is working.The goal is to find a desired (optimal) working point ν, such that properties of the QoIs are met.The restriction to one point is arbitrary, but allows a more focused discussion.This endeavor is usually accompanied with large and extensive multiobjective optimizations.
In an ideal world, one would run a large number of high fidelity simulations (in some proportion to the size of D) to solve the problem.However, even with state-of-the-art tools, and in cases of practical interest, it is impossible to accomplish this task due to the prohibitive time to solution.
With the help of adequate surrogate models, there are at least two ways to tackle the . Admissible design parameter search space D and one of the many possible ideal configurations ν (working point) of the accelerator.The red circles depict the training points, from which the surrogate model will be constructed.The equidistance of these points is not necessary; however, it is sufficient to introduce the overall concept.We furthermore assume that subspace D * is much smaller than D.
problem.First, with a high fidelity simulator we build a surrogate model from a coarser, discrete search space, depicted by the red points in Figure 2.With this surrogate model we then predict ν, which eventually yields an optimal solution.
In the second option, we would first find the smaller domain D * , with the help of the surrogate model constructed from D. Because D * is much smaller than D, it is feasible to use the expensive high-fidelity model to obtain ν ∈ D * .
It is important to mention that the surrogate model does not really reduce the search space.Rather, it is an approximation to the full model over the area of the search space where one believes that the model matters the most.The goal of the surrogate model is to create a cheap-to-sample approximation of the full model.

Mathematical bases of UQ.
We briefly introduce the mathematical bases in the style and the notation of [19,16,23,17,14]; more details can be found in Appendix A.
All square integrable, second-order random variables with finite variance output, u(ξ) ∈ L 2 (Ω, F, P), can be written as Hence α i denotes the deterministic coefficients and Ψ i (ξ) are the multivariate PC basis functions [19, section 10.1.1],[16], and i is a multi-index.Note that the uncertain QoI, u, is represented by a vector of deterministic parameters α i .Input uncertainties of the system have been discretized and approximated by the random vector with Ψ i (ξ) certain orthogonal bases functions, and I d,p a set of multi-indices.
The number K of PC basis functions of total order p in dimension d can be calculated to Because of the orthogonality of Ψ i k (ξ k ) and the independence of ξ k , as p → ∞, the truncated PCE in (2.2) converges in the mean-square sense if and only if the following two conditions are fulfilled: (1) u(ξ) has finite variance, and (2) the coefficients α i are computed from the projection equation [23] (2.3) .

Nonintrusive PCE.
In PC-based methods, one obtains the coefficients of the solution expansion either intrusively [33] or nonintrusively [34].An intrusive approach requires significant modification of the deterministic solvers and increases the number of equations to solve.
Nonintrusive methods, on the other hand, can make use of existing deterministic solvers (M) as black boxes.First, one needs to generate a set of N deterministic or random samples of ξ, denoted by {ξ (i) } N i=1 .The second step is to generate N realizations of the output QoI, {u(ξ (i) )} N i=1 , with the available deterministic solver M and without any solver modifications.The third and final step is to solve for the PC coefficients using the obtained realizations.Methods such as least-squares regression [35], pseudospectral collocation [17], Monte Carlo sampling [36], and compressive sampling [37] are available.Along these lines, an in-depth discussion on least-squares regression and compressive sampling can be found in [14, sections 3.1.1and 3.1.2].
The mean, E[•], and variance, Var[•], of u(ξ) can be directly approximated from the PC coefficients because of the polynomial basis orthogonality given by (2.4) A more complete description will be shown in section 2.5.The expensive, deterministic high-fidelity particle accelerator model, M, is described by a function u = M( x), where the input x is a point inside D (cf. Figure 2) and u is a vector of QoIs.Finding correlations in these high-dimensional spaces is nontrivial; however, it is vital for a deep understanding of the underlying physics.For example, reducing the search space is of great interest in the modeling and optimization process.In the spirit of Sobol' [38], let u * = M( x * ) be the sought (true) solution.The local sensitivity of the solution u * with respect to x k is estimated by (∂ u/∂x k ) x= x * .On the contrary, the global sensitivity approach does not specify the input x = u * ; it only considers the model M( x).Therefore, global sensitivity analysis should be regarded as a tool for studying the mathematical model rather than a specific solution ( x = x * ).For details, we refer the reader to Appendix A.1.

The UQTk-based framework.
In this section, a detailed description is provided on how the particle accelerator UQ framework is constructed.The framework is based on the uncertainty quantification toolkit (UQTk) [39], a lightweight C++/Python library that helps perform basic UQ tasks, including intrusive and nonintrusive forward propagation.The UQTk can also be used for inverse modeling via Bayesian or optimization techniques.The corresponding tools used from the UQTk are indicated in typewriter style in the following algorithm.
Let's denote M as the black-box solver, λ as the model parameters, and x as the design or controllable parameter, with l distinct values. 1 The nonintrusive propagation of uncertainty from the d-dimensional model parameter λ to the output u i = M( λ, x i ) follows a collocation procedure, given a K-dimensional basis Ψ = (Ψ 1 , . . ., Ψ K ) and K = (d+p)!d!p! multivariate basis terms, with p being the maximal polynomial order.Algorithm 1. generate for each x i (design or controllable) a PC surrogate model 1. generate N = (p + 1) d quadrature point-weight pairs ( ξ n , w n ) (generate quad).2. for each of quadrature point ξ n , compute corresponding model input λ n by 3. create the training points with high-fidelity simulations ( OPAL) 4. calculate the expectation via orthogonal projection (pce resp) using quadrature

5.
Given the computed α ki values for each i and k, one assembles the PCE Remark 1.The input PC in (2.6) is assumed to be given by an expert.For example, often only bounds for the inputs are known, in which case (2.6) is simply a linear PC or just scaling from Remark 2. If samples ξ n are randomly selected from the distribution of ξ, then the projection formula (2.8) still holds, as long as one sets w n = 1/N for all n, and it becomes an importance sampling Monte Carlo.
Remark 3. In Figure 3, a design parameter x is introduced.In case of p + 1 < l, i.e., if one only has a few discrete values for the design parameter, a reduced number model evaluation is obtained.Instead of sampling this parameter, you create l different response surfaces.) the utility pce eval can be used to evaluate ûi (2.9).

Sensitivity analysis.
As shown in section 2.4, the same information used in the surrogate model construction can be used in the sensitivity analysis.In the UQTk, pce sens will compute the total and joint sensitivities along with the variance fraction of each PC term individually.
3. A general model problem.Charged particle accelerators are among the largest and most complex scientific instruments.The application of charged particle accelerators ranges from material science to biology to fundamental physics questions, currently addressed, for example, with the LHC or in the future maybe with experiments like DAEδALUS/IsoDAR [40,25] (cf.section 5).There exist a wide range of different accelerator types and a commonly used classification in linear and circular types based on the geometrical nature.Given this taxonomy, a circular accelerator with a nonconstant radius of curvature is the most general accelerator and is used as a template for all other conventional types.Even in the simplest incarnation of a cyclotron [41], a rich dynamic is present, periodic or near periodic orbits or, in general, twist maps.
The Hamiltonian that describes the motion of a classical relativistic charged particle in a general magnetic field [42] is given by We are neglecting in this discussion the spin and radiation for the sake of simplicity.All external electromagnetic fields (magnets, etc.) are absorbed in the vector potential A, and with P we denote the generalized momenta.Charge and mass are denoted by m and q, respectively, and c is the speed of light.In the case of a cyclotron, all quantities are expressed in a Frenet-Serret coordinate system with nonconstant curvature h (cf. Figure 4).The scalar potential φ represents the nonlinear particle-particle interaction.The computation of φ and the resulting nonlinear force is computationally very expensive, and the effect of these forces is a limiting factor of high-intensity particle accelerators.The limiting aspect is based on the fact that these repulsive forces create a halo around the core of the particles.This halo has a tendency to separate from the core and contribute to particle losses, which in the end activates the machine to a level where maintenance is difficult or even impossible.The case of a cyclotron represents a large class of accelerator topologies.For example, in case of vanishing curvature h in (3.1) the case of a linear accelerator is recovered.We select the cyclotron in order to demonstrate the applicability of this framework in a very general context.
Solving such a problem under relevant circumstances is equivalent to solving a large Nbody problem with nontrivial boundary conditions.This, together with the multiscale nature of the problem-in time and phase space-is calling for a hierarchy of models.On the extreme end we are using a 1:1 ratio between simulation and macroparticles and computational times of days on high end parallel computers.Low-dimensional models, on the other hand, are important to narrow a potential high-dimensional search space and to make the problem more accessible.
Surrogate models as introduced in the previous section are in between the two extremes.They cover some nonlinearities but are much faster to evaluate compared to the high-fidelity model.With the sensitivity analysis we will get insight into a correlated space of QoIs and model parameters.
3.1.The accelerator simulation model.For this discussion, we briefly introduce OPALcycl [43], one of the three flavors of OPAL.OPAL will be used as the back-box solver denoted by M in (2.7).
3.1.1.Governing equation.In the cyclotron under consideration, the collision between particles can be neglected because the typical bunch density is low.In time domain, the general equations of motion for charged particles in electromagnetic fields can be expressed by We denote p = mcγβ as the momentum of a particle, β = (β x , β y , β z ) as the normalized velocity vector, and γ as the relativistic factor.In general, the time (t) and the position (x) dependent electric and magnetic vector fields are written in an abbreviated form as B and E.
If p is normalized by m 0 c, then (3.2) can be written in Cartesian coordinates as The evolution of the beam's distribution function, f (x, cβ, t) : (R M × R M × R) → R, can be expressed by a collisionless Vlasov equation: Here it is assumed that M particles are within the beam.In this particular case, E and B include both externally applied fields and space charge fields: all other fields are neglected.
3.1.2.The self-fields.The space charge fields can be obtained by a quasi-static approximation.In this approach, the relative motion of the particles is nonrelativistic in the beam rest frame, and thus the self-induced magnetic field is practically absent and the electric field can be computed by solving Poisson's equation Because of the large vertical gap in our cyclotron, the contributions from image charges and currents are minor compared to space charge effects [44], and hence it is a good approximation to use open boundary conditions.Details on the space charge calculation methods utilized in OPAL can be found in [43,45,46].

External fields.
With respect to the external magnetic field, two possible situations can be considered.In the first situation, the real field map is available on the median plane of the existing cyclotron machine using measurement equipment.
In most cases concerning cyclotrons, the vertical field, B z , is measured on the median plane (z = 0) only.Since the magnetic field outside the median plane is required to compute trajectories with z = 0, the field needs to be expanded in the z-direction.
According to the approach given by Gordon and Taivassalo [47], by using a magnetic potential and measured B z on the median plane at the point (r, θ, z) in cylindrical polar coordinates, the third-order field can be written as where B z ≡ B z (r, θ, 0) and All the partial differential coefficients are computed on the median plane data by interpolation, using Lagrange's 5-point formula.
In the second situation, a 3D field map for the region of interest is calculated numerically from a 3D model of the cyclotron.This is generally performed during the design phase of the cyclotron and utilizes commercial software.In this case, the calculated field will be more accurate, especially at large distances from the median plane; i.e., a full 3D field map can be calculated.For all calculations in this paper, we use the method of Gordon and Taivassalo [47].
For the radio-frequency cavities, a radial voltage profile V (r) along the radius of the cavity is used.The gap-width, g, is included in order to correct for the transit time.In addition, a voltage profile varying along the radius will give a phase compression of the bunch, which is induced by an additional magnetic field component B z in the gap,  Finally, in this paper, both the external fields and space charge fields are used to track particles for one time step using a fourth-order Runge-Kutta (RK) integrator.This means the fields are evaluated four times in each time step.Space charge fields are assumed to be constant during one time step because their variation is typically much slower than that of external fields.

Application of the UQ framework to a model problem.
To demonstrate the usefulness and strength of UQ, consider a simplified model of the PSI Injector 2 cyclotron, which is sketched in Figure 4.The simplifications are as follows: (1) only energies up to 8.5 MeV (turn 10) are considered to reduce the computational burden; (2) a Gaussian distribution, linearly matched to the injection energy of 870 keV, is used for the initial conditions; (3) the magnetic field and RF structures are the same as in our full production simulation; (4) P r and R are obtained from equilibrium orbit simulations; and (5) one collimator is introduced in order to mimic bunch shaping.Full scale high-fidelity simulations of this kind can be found in [48,22], where similar physics goals were pursued.

Model parameters.
In typical design studies of high-power cyclotrons, the high number of model parameters is such that one cannot fully scan their entire range.For this feasibility study, one model parameter out of a family of three important categories (cf. Figure 4) was chosen: 1. initial conditions: model parameter xp x , correlation between the initial x and p x phase-space variables; 2. collimator settings: model parameter ∆C 1 , position of the collimator; 3. RF phase settings: model parameter φ 1 defines the phase of the acceleration cavity.From previous experience, these three categories have the most influence when designing and optimizing high-precision models of high-power cyclotrons.The relationship of the parameters with uncertainties, λ 1 , λ 2 , λ 3 , is shown in Figure 3.

Quantities of interest (QoIs).
The phase space spanned by M macroparticles, in the high-fidelity OPAL model (simulation), is given by ( q i (t), p i (t)) ∈ Γ ⊂ R (2M +1) and i = x, y, z.We identify a subset of interesting QoIs such as 1. εx = q 2 x p 2 x − q x p x 2 the rms projected emittance and x the rms beam size; 2. the kinetic energy E and rms energy spread ∆E; 3. h t = q 4 x q 2 x 2 − k, the halo parameter in the x-direction at the end of turn t with k ∈ R, and a distribution dependent normalization constant.The rms beam size x is one of the better quantities that can be directly measured and hence among the first candidates for characterization of the particle beam.A measure of the projected phase-space volume is the emittance εx .This quantity is often used for the estimation of the beam quality.The two energy related parameters E and ∆E are target values to achieve.The first one, E, is closely related to the experiment, where the particle beam is designed for.The energy spread, ∆E, is directly related to the beam quality in the case of the presented model problem.Minimizing the halo of the particle beam is equal to minimizing losses, the most important quantity to optimize in high-power hadron accelerators.In the formulation of h t , this parameter is deviating from 1 if and only if the initial chosen distribution is changing.If the initial distribution is a stationary distribution, this measure can be attributed to the mechanism of halo generation, in the case of a deviation from the value 1.
In the case of a high-intensity cyclotron model, we choose the controllable parameter y as the average current.

UQ model setup.
The controllable parameters are not modeled with polynomials, but rather given by 10 equidistant values from 1 to 10 mA.As a next step, the polynomial type for the model parameter is chosen according to the Wiener-Askey scheme (cf.Appendix B).The distribution of the three model parameters xp x , ∆C 1 , and the phase φ 1 are modeled according to a uniform distribution using polynomials of the Legendre type.The bounds of the distribution are given in Table 1.
Other parameters for the UQ model are listed in Table 2.  Overall, the expected convergence is observed when increasing p as shown in Figures 5-10, and furthermore this is supported by the L 2 error shown in section 4.6.

4.4.1.
Projected emittance and beam size.Given the fact that the emittance is a very sensitive quantity, measuring phase space volume, it is surprising, but also promising, that such a good agreement between the surrogate model and the high-fidelity model can be achieved.This is graphically illustrated in Figures 5 and 6.The maximum training error in % is given in Tables 3 and 4 and is below 7% for all considered cases.

Table 3
Maximum training error in % between the high-fidelity and surrogate models for the projected emittance εx of the beam.High Fidelity OPAL Simulation High Fidelity OPAL Simulation  High Fidelity OPAL Simulation  2.
third harmonic cavity for acceleration).For the given experiment, only the last two turns are contributing.This fact is even better illustrated, when looking at the maximum training error, which is ≤ 0.07%, as seen in Table 5.

rms energy spread.
Despite the fact that the rms energy spread is influenced by space charge, the collimation, and the change in phase, a very good agreement with absolute deviations ≤ 5% was obtained.Table 6 and Figure 8 show details.
4.4.4.Halo parameters.The halo parameter was evaluated at turn 5 (Figure 9) and at turn 10 (Figure 10).As anticipated, the halo grows and the surrogate model has a maximum absolute error of ≤ 5%, again a very good accuracy.2.

Table 5
Maximum training error in % between the high-fidelity and surrogate models for the final energy of the beam.4.5.Sensitivity analysis.S k in (A.8) can be interpreted as the fraction of the variance in model M that can be attributed to the ith input parameter only.S T k in (A.9) measures the fractional contribution to the total variance due to the ith parameter and its interactions with all other model parameters.In what follows, an analysis based on S T k is shown for the model problem.
Figure 11 shows, for a subset of the controllable parameter I, sensitivities of the QoIs with respect to the model parameters.The polynomial order is p = 4; similar correlations for other   2.
orders are not shown.Correlations, for example the insensitivity of the energy, and x, p x or the significant energy phase correlation, are consistent with what is anticipated.A very mild dependence on x, p x is observed and expected.There is a phase correlation appearing in the case of I = 5 mA, which seems to be suppressed at other intensities, and the initial correlation of the distribution seems to become insignificant.A closer inspection of the phase space, beyond the scope of this article, hints that the halo at this intensity has a minimum.This could explain the observed behavior and is subject to a deeper investigation.
These are very interesting findings that can guide new designs but also improve existing accelerators, and show the quintessential merit and power of such a sensitivity analysis.2.
model.In addition, the accuracy was checked using a hold out model of N rs = 100 uniform random samples over the model parameter domain λ.

4.7.
Predictions.The surrogate model is constructed by selecting an appropriate number of training points in order to sample the input uncertainties of the design parameter space.These finite number of training points are depicted as yellow points in Figure 13.However, with the surrogate model we can choose any point within the lower and upper bounds specified (a i , b i in (2.10)) in order to obtain λ in (2.6).In Figure 13, the red points are arbitrarily chosen within the specified bounds and they are very well within the bounds of the surrogate model and the 95% confidence level (CL) obtained by evaluating the Student-t test.The data presented in Figure 13 are only from experiment 3 in the case of 1 mA.
4.8.Performance.The presented surrogate model is the most simple but gives, for the nontrivial model problem, statistically sound results.This fact and the remark that the evaluation of the surrogate model is ∼ 800× faster than the high-fidelity model (0.5 seconds vs. 400 seconds) opens up unprecedented possibilities in research areas such as on-line modeling and multiobjective [49,50] optimization of charged particle accelerators.4.9.Conclusions for the model problem.For a representative and at the same time nontrivial model problem, an accurate and fast to evaluate surrogate model is presented.From the sensitivity analysis, a phase correlation, in the case of I = 5 mA, could be observed.A surrogate model for the halo parameter with high fidelity is constructed.Due to the low computational cost of the surrogate mode, future optimization, minimization of the halo, is conceivable.This model problem should be understood as a "showcase" demonstrating the applicability of this approach in a generalized accelerator setting.
5. Contribution to the DAEDδALUS/IsoDAR accelerator design effort.The Decay-Atrest Experiment for δ CP violation At the Laboratory for Underground Science (DAEδALUS) [40] and the Isotope Decay-At-Rest experiment (IsoDAR) [25]   search for CP violation in the neutrino sector and "sterile" neutrinos, respectively.In order to be decisive within 5 years, the neutrino flux and, consequently, the driver beam current, produced by a chain of cyclotrons (cf. Figure 14), must be high, higher than achieved today.
With CP violation we recognize two broken symmetries: the broken C-symmetry (charge conjugation symmetry) and P-symmetry (parity symmetry).If CP-symmetry would not be broken, the law of physics is the same if a particle is interchanged with its antiparticle (C symmetry), while its spatial coordinates are inverted ("mirror" or P-symmetry).
The hypothetical sterile neutrinos are particles that interact only via gravity.To distinguish them from the known active neutrinos, the name sterile is stipulated.
5.1.Physics motivation.The standard model of particle physics includes three so-called "flavors" of neutrinos, ν e , ν µ , and ν τ , and their respective antiparticles.These particles can change flavor (neutrino oscillations), a process that can be described using a mixing matrix.This means that neutrinos must have a small mass [51].In addition, some experiments aimed at measuring these oscillations in more detail have shown anomalies that led to the postulation of "sterile" neutrinos which would take part in the oscillation, but, contrary to the three known flavors, do not interact through the weak force [52].Another important question is whether the three-neutrino model can give rise to a CP-violating phase δ CP [53], which might explain the matter-antimatter asymmetry in the universe today.The main challenge, from the accelerator point of view, is the handling of the highintensity beams.Of utmost importance is the minimization of particle losses and hence the understanding and mitigation of the particle halo.A second and related task is the optimization of the exit path out of the cyclotron.Here the separation of the last two turns  The following surrogate model construction and sensitivity analysis of the IsoDAR cyclotron is research in progress, i.e., far from complete, but should illustrate the potential of the introduced methods on an ongoing design effort.
5.2.Initial conditions for maximal energy and turn separation.In order to run the physics experiment with the highest efficiency, a target energy of 60 MeV/amu and lowest particle losses has to be reached.
A large turn separation between the extracted turn n and the turn n − 1 allows the insertion of a septum to change the sign of curvature of the nth orbit and hence facilitate clean (lossless) extraction of the beam.Detailed initial conditions for the cyclotron simulation are obtained from a 3D spiral inflector model, as shown schematically in Figure 15 [54].From the exit of the spiral inflector, we need to find optimal initial conditions for the full cyclotron favorable.We restrict the number of parameters to three model parameters describing the beam initial conditions, which are injection radius r, radial momenta p r (cf. Figure 1), and the phase φ of the radio frequency of the acceleration cavities (not shown in Figure 1).The controllable parameter is the angle θ of the spiral inflector, which brings the beam from the vertical direction into the midplane.We varied the azimuthal angle, θ, of the spiral inflector over a range of 5 deg as shown in Figure 15.With the guidance of the sensitivity analysis (cf. Figure 16), we selected the most favorable case of 140 degrees for the inflector angle to minimize the impact of p r on the turn separation (ts).The radial position r of the beam and the phase φ(r) of the cavities are directly accessible to control, while the radial momenta p r are not directly accessible for control, and hence a low sensitivity of the QoI with respect to p r is desired.
The influence of the spiral inflector position was known to have an impact on the extraction efficiency; a direct quantification of this fact was, to our knowledge, never described.
Having fixed the spiral inflector position, a surrogate model was constructed to estimate the final energy E and ts.We concentrate on the model for the energy and remark that the performance of the ts's is very similar.A detailed discussion about the ts as well as the influence of the spiral inflector position will be given in a forthcoming physics paper.
In Figure 17, a random sample N rs = 100 is used to compare the high-fidelity model to the surrogate model with orders 2 and 5.The second-order model is behaving very well until, at high energies, nonlinearities from the curvature of the RF sine wave are present.In Figure 17(b), the performance of the fifth-order model, in the high energy sector, is visible; in Table 7, the L 2 error is given.at higher energies and activate and/or damage the machine.For this purpose, collimators are inserted, just after injection (not shown in the figures).These collimators can be spatially adjusted and will deliberately remove particles with wrong dynamical properties, such as large vertical momenta.
The following multiobjective optimization problem needs to be solved: given a range of target emittance at extraction (indicate the quality of the beam), maximize the transmission (minimize losses).We remark that in order to solve this problem, many time consuming particle-in-cell simulations have to be conducted.Hence, an accurate surrogate model could have a substantial impact on the time to solution.
In order to construct such a model, we consider four collimators as model parameter λ and search for a surrogate model for the transmission Q = N inj /N ext × 100%, with N inj the injected number of particles and N ext the surviving (to be extracted) number of particles.
The second QoIs are the emittances defined in section 4.2.As in the previous section, we use a uniform random sample N rs = 100 to evaluate the quality of the surrogate model.
In Figure 18, we recognize a very good agreement between the random samples (from highfidelity OPAL simulations) and the surrogate model.The same is true for the beam quality shown in Figure 19.Hence we can conclude that for the four QoIs, energy E, transmission Q, and emittances ε x and ε y , we can construct high-fidelity surrogate models.
One PIC simulation, used to train the model, runs for approximately 13090 seconds on 8 cores (Intel KNL).The evaluation of the fourth-order multivariate polynomial takes less than 0.002 seconds using the UQTk software.This represents a speedup of ≈ 6.5 × 10 6 and allows us to do large-scale multiobjective optimization, using surrogate models as a forward solver.These optimizations are part of ongoing research for the IsoDAR compact cyclotron design.

Conclusions.
A sampling-based UQ approach is presented to study, for the first time, the effects of input uncertainties on the performance of particle accelerators.A particular, but complex, example in the form of a high-intensity cyclotron was used to demonstrate the usefulness of the surrogate model as well as the global sensitivity analysis via computing the total Sobol' indices.The presented physics problem is a model problem, with the aim of demonstrating the usefulness and applicability of the presented UQ approach.However, we claim to present a problem that can be recognized as a template for many high-intensity modeling attempts and beyond.
The proposed UQ approach is based on PCE using the UQTk software.The goal is to achieve an accurate estimation of solution statistics using a minimal number of high-fidelity simulations.For several QoIs, a surrogate model was constructed; the validity is proved by comparing to a high-fidelity model.L 2 error norms show the expected convergence with regard to the degree of the PCE.For the rms beam size (x), holdout points, i.e., points that are not used in the training set, were evaluated and compared to the statistical expectations from the model.We found that the values are consistent with the surrogate model and clearly within the 95% CL.

Figure 1 .
Figure 1.Spiral inflector with particle trajectories, from the IsoDAR example presented in section 5.

Figure 3 .
Figure 3.The UQ framework, with the discretized input uncertainties of the system denoted by ξ; cf.(2.1).In case of different design (or controllable) parameters x, we would build l separate response surfaces.Details can be found in Algorithm 1 and Appendix A.
and ASA.Published by SIAM and ASA under the terms of the Creative Commons 4.0 license Downloaded 01/16/20 to 129.129.214.179.Redistribution subject to CCBY license where φ and ρ are the electrostatic potential and the spatial charge density in the beam rest frame.The electric field can then be calculated by (3.6) E sc = −∇φ and transformed back to yield both the electric and the magnetic fields, in the lab frame, as required in (3.4) by means of a Lorentz transformation.
cos(ω rf t − φ), c 2019 SIAM and ASA.Published by SIAM and ASA under the terms of the Creative Commons 4.0 license Downloaded 01/16/20 to 129.129.214.179.Redistribution subject to CCBY license with F denoting the transit time factor (F = 1 2 ω rf ∆t) and ∆t the transit time defined by (3.10) ∆t = g βc .

Figure 4 .
Figure 4.The cyclotron model problem setup.The two red lines indicate the two double gap flat-top resonators, the blue line represents a collimator, and the yellow circle stands for the initial conditions.

4. 4 .
High-fidelity simulations vs. surrogate model.As a first method to determine the validity of the surrogate model, the values of the high-fidelity OPAL simulations on the xaxis and the values of the surrogate model on the y-axis were compared.The distance of the corresponding point to the line x = y is a measure of the surrogate model's quality.The QoIs, c 2019 SIAM and ASA.Published by SIAM and ASA under the terms of the Creative Commons 4.0 license Downloaded 01/16/20 to 129.129.214.179.Redistribution subject to CCBY license

P = 4 P = 3 P = 2 I 2 . 4 2 I 3 3
Final energy.The energy dependence shown in Figure 7 for 10 mA serves as an illustration of the expected behavior for all other intensities.This is because of the small gain the third harmonic cavity is supposed to deliver (in the PSI Injector 2, we use the c 2019 SIAM and ASA.Published by SIAM and ASA under the terms of the Creative Commons 4.0 license Downloaded 01/16/20 to 129.129.214.179.Redistribution subject to CCBY licenseTable Maximum training error in % between the high-fidelity and surrogate models for the rms beam size x of the beam.P = 4 P = 3 P = ˜ x (mm−mr) Experiment 2 ˜ x (mm−mr) Experiment 1 ˜ x (mm−mr) ˜ x (mm−mr) Experiment 2 ˜ x (mm−mr) Experiment 1 ˜ x (mm−mr)

Figure 6 .
Figure 6.The rms beam size x (mm) for all three experiments described in Table2.

Figure 7 .
Figure 7. Final energy E (MeV) for I = 10 mA and all experiments described in Table2.

Figure 8 .
Figure 8. Energy spread ∆E (keV) for all three experiments described in Table2.

Figure 12 .
Figure 12.Medium values and variances are shown as dots and error bars on the left y-axis for the extraction energy E. The global L2 error (lines) between the high-fidelity and the surrogate models, for the final energy of the particle beam, is shown on the right y-axis.

Figure 13 .
Figure 13.The surrogate mode for x, together with training and prediction points.The 95% CL of the model is also shown.

Figure 14 .
Figure 14.Cartoon picture of one single DAEδALUS module with the injector part that can be used for IsoDAR highlighted on the right.The spiral inflector (cf.Figure15) would be located in the middle of the DIC, at the end of the line marked with LEBT (low energy beam transfer).
in the cyclotron has to be maximized.The conducted research by the DAEδALUS/IsoDAR collaboration over the last couple of years suggests that it is feasible, albeit challenging, to accelerate 5 mA of H + 2 to 60 MeV/amu in a compact cyclotron and boost it to 800 MeV/amu in the DSRC (DAEδALUS Superconducting Ring Cyclotron) with clean extraction in both cases.

Figure 16 .Figure 17 .
Figure 16.Global sensitivity analysis for finding the most favored spiral inflector position.On the left side, the sensitivities for a spiral inflector position of θ = 135 deg is shown, while on the right side a more favorable position of θ = 140 deg, reducing the influence of pr on the turn separation ts.

5. 3 .
Maximal transmission.Inspection of Figure 1 reveals the fact that particles will terminate at some location very early in the machine.Particles with wrong dynamical properties need to be removed from the ensemble at low energies; otherwise, we would lose them c 2019 SIAM and ASA.Published by SIAM and ASA under the terms of the Creative Commons 4.0 license Downloaded 01/16/20 to 129.129.214.179.Redistribution subject to CCBY license

Figure 18 .Figure 19 .
Figure 18.Surrogate model for the transmission Q, compared with random sampling.

Table 1
Upper and lower bounds of the design parameters.

Table 2
Summary of UQ related parameters for the presented results.The dimension for all the experiments is d = 3.The one controllable parameter x has length l = 4.
as defined in section 4.2, are compared for a subset of controllable parameters, 1, 5, 8, and 10 mA, and for three different orders of the surrogate model, as described in Table2.All data from the surrogate model and the high-fidelity model are taken at the end of turn 10 in our model problem.The maximum training error is calculated from the dataset used to create the surrogate model.

Table 2 .
c 2019 SIAM and ASA.Published by SIAM and ASA under the terms of the Creative Commons 4.0 license Downloaded 01/16/20 to 129.129.214.179.Redistribution subject to CCBY license

Table 6
Maximum training error in % between the high-fidelity and surrogate models for the energy spred ∆E of the beam.
c 2019 SIAM and ASA.Published by SIAM and ASA under the terms of the Creative Commons 4.0 license Downloaded 01/16/20 to 129.129.214.179.Redistribution subject to CCBY license are proposed experiments to c 2019 SIAM and ASA.Published by SIAM and ASA under the terms of the Creative Commons 4.0 license Downloaded 01/16/20 to 129.129.214.179.Redistribution subject to CCBY license

Table 7
L2 error in energy as a function of the order of the surrogate model.