Differentiable self-consistent space-charge simulation for accelerator design

The nonlinear space-charge effects in a high intensity or high brightness accelerator can have a significant impact on the beam properties through the accelerator. These effects are included in the accelerator design via self-consistent multi-particle tracking simulations. In order to study the sensitivity of the final beam's properties with respect to the accelerator design parameters, one has to carry out the time-consuming space-charge simulation multiple times. In this paper, we propose a differentiable self-consistent space-charge simulation model that enables the study of the beam sensitivity through only one simulation. Such a model can also be used with gradient-based numerical optimizers for accelerator design optimizations including the self-consistent space-charge effects.


I. INTRODUCTION
The nonlinear space-charge effects from the Coulomb interaction inside a charged particle beam can have a significant impact on the beam transport through an accelerator by causing beam emittance growth, halo formation, and even particle losses. These effects are normally studied in the accelerator design via self-consistent simulations. To simulate the space-charge effects self-consistently, multi-particle tracking with a particle-in-cell method has been widely employed in the accelerator community [1][2][3][4][5][6][7][8][9][10][11]. However, none of these codes directly calculates the derivatives of the beam property with respect to the accelerator machine parameters.
The derivatives of the beam property with respect to the accelerator machine parameters are important in the accelerator design. These derivatives provide quantitative measurements of the sensitivity of the final beam property with respect to the machine parameters. Such sensitivities can be used to set the tolerance limits of the corresponding machine parameters. In a typical machine design, the sensitivity can be obtained by running the space-charge simulation multiple times, each time with a small change of a single machine parameter. A numerical scheme such as the finite difference method is used to calculate the derivative of the beam property with respect to that parameter. In an accelerator, there can be hundreds and thousands machine parameters. To compute the derivatives with respect to all machine parameters using the finite difference method will involve hundreds and thousands self-consistent spacecharge simulations and can be very time consuming. In this paper, we propose a differentiable self-consistent spacecharge simulation model. The derivatives of the final beam property with respect to the entire machine parameters can be obtained through a single simulation. Such a differentiable space-charge simulation can also be integrated into a gradient-based optimizer for accelerator design optimization.
The differentiable simulator is a simulator that can automatically compute derivatives of the simulation result with respect to its input parameters. This can be done through automatic differentiation [12] that has been widely used in the artificial intelligent/machine learning (AI/ML) community to train a neural network through gradientbased optimization method. A number of AI/ML frameworks such as PyTorch [13] and TensorFlow [14] include this capability. Recently, a differentiable model was developed to optimize magnet design by implementing the nonparametric version of the Preisach model output using the PyTorch framework [15]. In this paper, instead of using the automatic differentiation framework of the AI/ML community, we proposed a differentiable self-consistent spacecharge simulation model using the truncated power series algebra that was developed in the accelerator community. This method can also be implemented in many other available self-consistent space-charge codes for high intensity, high brightness accelerator design study.
The truncated power series algebra (TPSA) was first introduced by Berz in 1989 for accelerator applications [16]. Since then, it has been used to calculate transfer maps of beam line elements in several optics codes [17][18][19][20]. It was also used to calculate the space-charge potential/fields in the fast multipole method [21][22][23] and to solve the Poisson's equation [24]. Recently, it was also used to extract transfer maps in the presence of space-charge effects [25][26][27]. However, so far, these applications of the TPSA have no direct connections to the control parameters in the accelerator design. In this study, we present a new application of the TPSA to the accelerator design in the self-consistent spacecharge simulation. In this differentiable space-charge model, the final beam property such as emittance is connected to the accelerator control parameters. This enables to study the sensitivity of the final beam property with respect to these control parameters through only one simulation, which is different from the other uncertainty quantification methods such as the surrogate model method in which a number of simulations have to be used to train the model [28].
The organization of this paper is as follows: after the introduction, we introduce the truncated power series algebra in Section II, present a differentiable space-charge model in Section III, two application examples in Section IV, and draw conclusions in Section VI.

II. TRUNCATED POWER SERIES ALGEBRA
The truncated power series algebra (TPSA) is an effective tool to calculate derivatives of a function with respect to its variables using an algebraic method. Consider the Taylor series approximation of a one-dimensional function f (x) at a point x 0 , the derivatives in the above equation can be calculated using a numerical finite difference method, for example, Such a way to calculate the derivative introduces numerical errors and requires multiple function evaluations for a multi-variable function. In the truncated power series algebra method, the derivative up to N th order can be regarded as a point in a function space spanned by the bases: These derivatives can be represented as a vector: Such a vector is also called a TPSA variable. For example, for a constant c, its derivative representation is Dc = [c, 0, 0, · · · , 0], and for a variable x, Dx = [x, 1, 0, · · · , 0]. By using the above derivative vector representation, the derivatives of function with respect to its variable can be written as the function of that derivative vector, i.e.
The computing of the derivatives of a function becomes the algebraic function evaluations. The evaluation of the derivative vector inside a function can be broken down as the operations of addition and multiplication. Given two derivative vectors , the sum of two vectors will be: The multiplication of these two vectors will be, where Using the above rules of addition and multiplication, the operation of a derivative vector inside a function can be calculated algebraically. For example, the reciprocal of a derivative vector 1/Df x0 can be calculated as: For a concrete example, f (x) = 1 1+x+x 2 , one can use the above derivative vector operations to obtain the first and the second derivative of this function at x = 1. That is, x = 1, Dx = [1, 1, 0], and This yields the first derivative f (1) = − 1 3 and the second derivative f (1) = 4 9 . The truncated power series algebra changes the calculation of the derivatives of a function with respect to its individual variable into the evaluation of a function of derivative vector, i.e. a function of TPSA variable.
The above single variable function example can be extended to a multiple variable function with the more complicated multiplication rule for the derivative vector [29]. The truncated power series algebra libraries have been developed to handle some general special functions such as the exponential function and the trigonometry function [30,31].

III. DIFFERENTIABLE SELF-CONSISTENT SPACE-CHARGE MODEL
The above truncated power series algebra is used to develop a differentiable self-consistent space-charge model using TPSA variables. The self-consistent space-charge simulation is done by solving the following Hamilton equations: where H(r 1 , p 1 , r 2 , p 2 , · · · , s) denotes the Hamiltonian of the system, and r i and p i denote canonical coordinates and momenta of particle i, respectively. Let ζ denote a 6N or 4N-vector of coordinates, the above Hamilton's equation can be rewritten as: where [ , ] is the Poisson bracket. A formal solution for the above equation after a single step τ can be written as: Here, we have defined a differential operator : H : as : where A z denotes the longitudinal vector potential associated with the external focusing fields and where K = qI/(2π 0 p 0 v 2 0 γ 2 0 ) is the generalized perveance, q is the charge of particle, I is the beam current, 0 is the dielectric constant in vacuum, p 0 is the momentum of the reference particle, v 0 is the speed of the reference particle, γ 0 is the relativistic factor of the reference particle, and ϕ is the space charge Coulomb interaction potential. In this Hamiltonian, the effects of the direct electric potential and the longitudinal vector potential are combined together. For a Hamiltonian that can be written as a sum of two terms H = H 1 + H 2 , an approximate solution to the above formal solution can be written as: Let exp(− 1 2 τ : H 1 :) define a transfer map M 1 and exp(−τ : H 2 :) a transfer map M 2 , for a single step, the above splitting results in a second order numerical integrator for the original Hamilton's equation as: For the external focusing with quadrupole magnets, the single step transfer map M 1 in the focusing plane can be written as: and in the defocusing plane as: where k is the normalized focusing strength k = qg/p 0 and g is the magnetic field gradient. For the space-charge Hamiltonian H 2 (r), the single step transfer map M 2 can be written as: The electric Coulomb potential in the Hamiltonian H 2 can be obtained from the solution of the Poisson equation. In the following, we assume that the coasting beam is inside a rectangular perfectly conducting pipe. In this case, the two-dimensional Poisson's equation can be written as: where φ is the electric potential, and ρ is the particle density distribution of the beam. The boundary conditions for the electric potential inside the rectangular perfectly conducting pipe are: where a is the horizontal width of the pipe and b is the vertical width of the pipe. Given the boundary conditions in Eqs. 23-26, using a Galerkin spectral approximation method, one obtains the space-charge Hamiltonian H 2 as [32]: The resultant one-step symplectic transfer map M 2 of the particle i with this Hamiltonian is given as: where the space-charge potential in the spectral domain is given as: Here, both p xi and p yi are normalized by the reference particle momentum p 0 .
In the differentiable space-charge simulation, the above space-charge model will be rewritten using TPSA variables. The phase space coordinates r i and p i of the particle i will be replaced by the corresponding TPSA variables Dr i and Dp i defined in the last section. The potential φ lm is replaced by Dφ lm of the TPSA variable. The momentum updates after a single step due to the space-charge effects are given by: Dφ lm α l cos(α l Dx i ) sin(β m Dy i ) Dφ lm β m sin(α l Dx i ) cos(β m Dy i ) (30) where the space-charge potential in TPSA variable is: The map corresponding to the external quadrupole field can also be written in TPSA variable as: where Dk is the TPSA variable of the quadrupole focusing strength, and Dτ is the TPSA variable of the step size that is the quadrupole length divided by the number of steps. A similar expression can be written in the defocusing plane of the quadrupole. The charged particle beam initial distribution parameters, beam energy, and current can also be written using TPSA variables if needed. The final beam properties such as emittances are defined using the TPSA variables. The horizontal emittance D x is given as where The vertical emittance has a similar expression with x replaced by y.

IV. APPLICATION EXAMPLES
In the following, we will use two application examples to illustrate the above differentiable self-consistent spacecharge simulation model: one for parameter sensitivity study, and the other one for design optimization study.

FIG. 2:
Derivatives of final horizontal and vertical emittances with respect to seven lattice parameters from the single differentiable self-consistent space-charge simulation and from the finite difference approximation to the first derivative.
In the first illustrative application, we studied sensitivity of the final emittances of a 1 GeV coasting proton beam transporting through a transverse focusing FODO lattice inside a rectangular perfectly conducting pipe with respect to the lattice parameters and the initial beam parameters. A schematic plot of the FODO lattice is shown in Fig. 1. It consists of a drift (D1) of 0.2 meters, a quadrupole (Q1) of 0.1 meters and focusing strength 29.6/m 2 for transverse focusing, another drift (D2) of 0.4 meters, another quadrupole (Q2) of the same length as the first quadrupole but opposite sign of focusing, and another drift (D3) of length 0.2 meters. The rectangular pipe has an aperture size of 13 by 13 millimeters. The zero current phase advance of the FODO lattice is 87.0 degrees. The current of the proton beam is 200 Amperes with 1 mm mrad normalized emittance, which results in a depressed phase advance of 63.1 degrees. The initial distribution is assumed to be a four-dimensional Gaussian distribution given by: The parameters in the above distribution are chosen to be RMS matched through the FODO lattice including the space-charge effects. We first checked the sensitivity of the final emittances with respect to seven FODO lattice parameters, i.e. lengths of the drift and quadrupole elements, and focusing strengths of two quadrupole elements. In the this example, we used 5000 macroparticles and 12×12 spectral modes in the space-charge simulation. The sensitivity is measured by the first derivatives of the final emittances (normalized by initial emittance) with respect to these lattice parameters. Figure 2 shows the sensitivities of the final horizontal and vertical emittances with respect to the seven lattice parameters from the single differentiable self-consistent space-charge simulation and from the numerical finite difference approximation to the first derivative using eight space-charge simulations. Here, the seven lattice parameters are drift one length (D1 L), quadrupole one length (Q1 L), quadrupole one focusing strength (Q1 k), drift two length (D2 L), quadrupole two length (Q2 L), quadrupole two strength (Q2 k), drift three length (D3 L). It is seen that the first derivatives calculated from the differentiable space-charge model and from the finite difference approximation agree with each other quite well. This provides a verification of the differentiable self-consistent space-charge simulation model. From this figure, one can see that the final emittance is much more sensitive to the length of the first quadrupole than the other lattice parameters such as drift lengths and quadrupole strengths. The change of the first quadrupole length causes significant beam envelope variation and results in large final emittance change.
Next, we checked the sensitivity of the final proton beam emittance with respect to the initial beam parameters. These beam parameters are proton beam energy (Eng.), beam current (Cur.), and six beam distribution parameters (σ x , σ px , µ xpx , σ y , σ py , and µ ypy ). These parameters were represented using TPSA variables. A number of macroparticle coordinates (also in TPSA variables) were generated from the regular sampling method using the beam distribution parameters. The quadrupole magnetic field gradient was used instead of the focusing strength since the latter includes the proton beam energy. The change of the beam energy affects both the transverse focusing and the space-charge effects. Figure 3 shows the first derivatives of the final emittances with respect to the above eight parameters. It is seen that the final emittances are more sensitive to the initial distribution parameter σs. The perturbation of these parameters affects the matching of the initial distribution to the FODO lattice and can cause significant emittance growth with a mismatched beam [33,34]. The differentiable self-consistent space-charge simulation through the accelerator lattice produces the final beam properties and their derivatives with respect to the lattice control parameters at the exit of the accelerator as shown in the first example. These derivatives can be used in a gradient-based parameter optimizer for accelerator lattice control parameter optimization. In the following example, we integrated the differentiable self-consistent space-charge simulation model into a conjugate gradient optimizer to attain the quadrupole strengths inside a matching section in front of a periodic FODO lattice.
A schematic plot of the matching section lattice and the periodic FODO lattice is shown in Fig. 4. The first four quadrupoles in the figure were used to match an initial distribution to the given Twiss parameters at the entrance to the periodic FODO lattice. The Polak-Ribiere conjugate gradient optimization method [35] was used to minimize the objective function that is defined as follows: where k is a set of control variables, α xt , β xt , α yt , and β yt are the target Twiss parameters at the entrance to the periodic lattice, and the α x , β x , α y , and β y are the beam Twiss parameters calculated from the self-consistent spacecharge simulation. These calculated beam Twiss parameters depending on the focusing strengths of the quadrupoles inside the matching section. These strengths are control variables in the above objective function. Using the above differentiable space-charge simulation model, the first derivatives of the objective function with respect to the four control variables were obtained in addition to the objective function value. These derivatives were used to construct a conjugate direction of the gradient direction to guide the search for the minimum solution. Figure 5 shows the proton beam transverse RMS size evolution through the above FODO lattice without the quadrupole matching and with the quadrupole matching including the space-charge effects. Here the quadrupole strengths inside the matching section without the matching were set based on the zero current matched solution. It is seen that with a 200 Ampere beam, the space-charge effects are significant so that the initial zero-current matched quadrupole strengths no longer produce a matched RMS evolution inside the periodic lattice. After reoptimizing the four quadrupole strengths including the space-charge effects through the self-consistent simulations, the RMS evolution inside the periodic lattice becomes well-matched and results in much less emittance growth (less than 10%) than that from the mismatched quadrupole setting (greater than 50%).

V. CONCLUSIONS
The self-consistent space-charge simulation is an important part in the high intensity, high brightness accelerator design. In this paper, we proposed a differentiable self-consistent space-charge model that can be used to efficiently study the sensitivity of the final beam properties with respect to the accelerator design parameters. Using the differentiable self-consistent space-charge model, only one simulation is needed to attain the derivatives of the final beam properties with respect to all accelerator design parameters instead of multiple simulations of the conventional space-charge model. The resultant first derivatives measure the sensitivity of the final beam properties with respect to those design parameters. Some highly sensitive machine parameters can be quickly identified after one differentiable self-consistent space-charge simulation.
As an illustration, we presented two application examples. One example computed the sensitivities of the final proton beam emittances with respect to seven lattice parameters or eight beam parameters through a single differentiable self-consistent simulation. The second example showed that the derivatives of the objective function with respect to the quadrupole strengths inside a matching section from the differentiable self-consistent space-charge simulation were used in a conjugate gradient optimizer to attain the space-charge matched solution to a periodic FODO lattice. The success of these two examples shows that the differentiable self-consistent space-charge simulation model can be a useful tool in the accelerator design.
In this study, we used a spectral solver as an illustration of a differentiable space-charge model. In general, some other space-charge solvers such as fast multipole solver or Green function solver can also be used to as the differentiable spacecharge model as long as the space-charge fields from these solvers are represented using TPSA variables. Furthermore, this can be generalized beyond the space-charge effects. The other collective effects such as beam-beam effects can also be represented using TPSA variables. Together with the particle coordinates and accelerator machine parameters that are represented using the TPSA variables, one can develop a general differentiable simulation tool that includes a variety of collective effects for accelerator design. Such a tool will be useful to study the beam quality sensitivity to the accelerator machine parameters subject the collective effects through only one simulation.
In the present study, the computational speed of the above differentiable space-charge solver is slow compared with the conventional space-charge solver due to the overhead associated with the computation using TPSA variables in the available TPSA library. From the discussion with the author of the TPSA library used in this study, the computational speed can be substantially improved if only the first-order derivative is needed [36]. A number of performance optimization strategies can be employed to improve the computational efficiency of the library. This will be pursued in the future study by working with the TPSA library developers.