Applications of Differentiable Physics Simulations in
Particle Accelerator Modeling

Ryan Roussel
SLAC National Accelerator Laboratory
Menlo Park, CA 94025
rroussel@slac.stanford.edu

Auralee Edelen
SLAC National Accelerator Laboratory
Menlo Park, CA 94025
edelen@slac.stanford.edu

Abstract
Current physics models used to interpret experimental measurements of particle
beams require either simplifying assumptions to be made in order to ensure analytical tractability, or black box optimization methods to perform model based
inference. This reduces the quantity and quality of information gained from experimental measurements, in a system where measurements have a limited availability.
However differentiable physics modeling, combined with machine learning techniques, can overcome these analysis limitations, enabling accurate, detailed model
creation of physical accelerators. Here we examine two applications of differentiable modeling, first to characterize beam responses to accelerator elements
exhibiting hysteretic behavior, and second to characterize beam distributions inhigh dimensional phase spaces.

1

Introduction

Particle accelerators are indispensable tools for making discoveries in the physical, material, chemicaland biological sciences. Detailed physics models have been developed to describe the generation andmanipulation of particle beams for these applications. Despite this wealth of knowledge, accurately
replicating realistic conditions in the physical accelerator is challenging due to the limited availabilityof measurements and the complex nature of real accelerator components and beam distributions. In
this work, we examine the application of differentiable computing techniques combined with machinelearning towards enabling detailed, accurate modeling of realistic accelerator elements and beam
properties.

2

Hysteresis Modeling

Hysteresis is a well-known physical phenomenon where the state of a given system is dependent
on its historical path through state-space. Hysteresis effects in magnetic [1], mechanical [2] andmaterial [3] elements of particle accelerators makes optimizing the performance of current accelerator
facilities used for scientific discovery challenging using standard black box optimization algorithms,
such as Bayesian optimization (BO) [4].
The Preisach model [5, 6] is a numerical representation of hysteresis, comprised of a discrete setof hysterons, which when added together, model the output of a hysteretic system f (t) for a timedependent input u(t). Given a set of discrete time ordered inputs ui = u(ti ), the hysteron state isrepresented by the hysteron operator γ̂αβ shown in Fig. 1a, which has an output of ±1, where α andβ describe the input required to switch the hysteron between its two possible states. The number ofhysterons with values (α, β) is given by the hysteron density function µ(α, β), plotted on the Preisach(α-β) plane (Fig. 1b).
Machine Learning and the Physical Sciences workshop, NeurIPS 2022.

Figure 1: Elements of the differentiable non-parametric Preisach hysteresis model. (a) Output of
the hysteron operator γ̂αβ and the approximate differentiable hysteron operator ζ̂αβ acting on the
input u(t). (b) Discretization of the density on the Preisach (α-β) plane. Note that µ(α, β) = 0 if
α < β, α > αm or β < βm where αm , βm are equal to the maximum and minimum inputs of themodel respectively. (c) S + and S − sub-domains after three time steps, where u1 > u3 > u2 > βm ,assuming that all hysterons are in the negative state initially. Inset: Corresponding model output.

The Preisach model output is represented byZZ
f (t) = Γ̂u(t) =µ(α, β)γ̂αβ u(t)dαdβ

(1)

α≥β

where α ≥ β results from physical conditions of the hysteron operator. This integral is evaluatedthrough a geometric interpretation, shown in Fig. 1(c). Given the sequence of input values ui , we candetermine sub-regions of the Preisach plane, S + and S − , where hysteron operators output positive
and negative states respectively using a geometric interpretation.
Fitting a Preisach model to experimental data requires the determination of the hysteron density
function µ(α, β), often referred to as the identification problem. We discretize the hysterion density
function onto a fine mesh µi = µ(αi , βi ) and treat the density at each mesh point as a free tuning
parameter to fit experimental measurements. Due to the large number of mesh points and corresponding free parameters, gradient enabled optimization must be used to fit the model to experimental
measurements. This is achieved by implementing the Preisach hysteresis calculation in PyTorch [7],which uses backwards auto-differentiation to efficiently calculate gradients.
Our goal is to determine hysteresis properties from beam-based measurements, so we combine thedifferentiable hysteresis model with a Gaussian process (GP) model [8], implemented in GPyTorch [9],to describe beam propagation as a function of magnetic fields due to hysteresis. We simultaneously
infer hysteresis model parameters and GP hyperparameters from training data by maximizing a
differentiable calculation of the log likelihood with respect to the parameters of the joint modelthrough gradient descent.
We demonstrate the effectiveness of our joint hysteresis-GP (H-GP) model by fitting the beam
response with respect to the current applied to a focusing magnet located in the Advanced Photon
Source (APS) injector [10]. Measurements from this experiment, shown in Fig. 2, have two sourcesof uncertainty, one from random noise inherent in the accelerator (aleatoric uncertainty) and one dueto the unknown properties of magnetic hysteresis (epistemic uncertainty). We trained two models
on the entire data set using L-BFGS with gradient information, the results of which are shown in
Fig. 2. First, we trained a normal GP model (Fig. 2a) which does not take into account the existenceof hysteresis. As a result it interprets epistemic errors due to hysteresis as aleatoric uncertainty,
overestimating uncertainties in portions of the input domain. However, the joint hysteresis-GP model(Fig. 2b), is able to resolve hysteresis cycles inside the data, removing epistemic uncertainties in
the model prediction, thus improving model accuracy and reducing uncertainty. The increase inaccuracy from joint hysteresis-GP models has ramifications for model-based, online optimization ofaccelerators using BO [11].
2Figure 2: Comparison between GP modeling and joint hysteresis-GP modeling of beam transmission
as a function of quadrupole current at the APS injector. (Left) GP model prediction with training dataover three cycles (see inset). Shading denotes 2σ confidence region. (Right) Hysteresis-GP model
prediction, colored by cycle index, compared to confidence interval of normal GP model.

3

Phase Space Reconstruction

Tomographic measurement techniques are used in accelerators to determine the density distribution ofbeam particles in phase space ρ(x, px , y, py , z, pz ) from limited measurements [12, 13, 14, 15, 16, 17].
While these methods have been shown to effectively reconstruct 2D phase spaces from image
projections using algebraic methods, application to higher-dimensional spaces requires independence
assumptions between the phase spaces of principal coordinate axes (x, y, z), complicated phase spacerotation procedures [18], or measurement of multiple 2D sub-spaces with specialized diagnostic
hardware [19]. It is also possible to directly fit a beam distribution to experimental data through theuse of black box optimization algorithms [20, 21, 22] and particle tracking simulations. This however
increases the computational cost of inferring the beam distribution from simulation, resulting in
simplified beam representations to keep costs within a fixed computational budget.
We propose reducing the computational cost of determining beam phase space distributions by
combining methods for parameterizing arbitrary beam probability distributions in 6D phase space
using neural networks with differentiable particle tracking simulations. Particle tracking simulations
propagate samples from the initial beam distribution down a model of the beamline to simulated
diagnostics. This allows us to learn the beam distribution from arbitrary downstream accelerator
measurements [23]. Our beam distribution parameterization is heavily inspired by normalizing
flows [24], however we are only concerned with forward transformations of drawn samples from abase distribution, and not the transformed probability distribution itself (which requires invertabletransformations). Thus we use simple densely connected neural networks instead of a series of
invertable transformations.
Fitting neural network parameters to experimental measurements is done by minimizing a loss
function to determine the most likely initial beam distribution, subject to the constraint that it
reproduces experimental measurements, similar to the maximum entropy tomography algorithm[25]. The likelihood of an initial beam distribution in phase space is maximized by maximizing
the distribution entropy, which is proportional to the log of the 6D beam emittance ε6D [26]. Thus,
we specify a loss function that minimizes the negative entropy of the proposal beam distribution,
penalized by the degree to which the proposal distribution reproduces measurements of the transverse
beam distribution at the screen location. To evaluate the penalty for a given proposal distribution, we
track the proposal distribution through a batch of accelerator simulations that mimic experimental(i,j)
conditions to generate a set of simulated images Qn to compare with experimental measurements.
As in the previous example, a loss function that is differentiable is needed to fit the large number ofneural network parameters. Unfortunately, no particle tracking codes which are differentiable withrespect to particle coordinates currently exist. To satisfy this requirement, we implement particle
tracking and image creation using PyTorch [7]. We demonstrated reconstruction of a beam in phase3

Figure 3: Description of our approach for reconstructing 6D phase space beam distributions. First,
a base distribution is transformed via a neural network, parameterized by θ, into a proposed initial
distribution. This distribution is then transported through a differentiable accelerator simulationof the tomographic beamline. The quadrupole is scanned to produce a series of images on thescreen, both in simulation and on the operating accelerator. The images produced both from the
(i,j)
(i,j)
simulation Qn and the accelerator Rn are then compared with a custom loss function, which
when optimized, maximizes the entropy of the proposal distribution constrained on accurately
reproducing experimental measurements.

space using a synthetic test case, shown in Fig. 4. We see excellent agreement between the averagereconstructed and synthetic projections in both transverse correlated and uncorrelated phase spaces,
as well as transverse beam images.

Figure 4: (Left) Comparisons between the synthetic and reconstructed beam probability distributions
using our method. A synthetic beam with a tailored phase space distribution is transported througha simulated diagnostic beamline with varying quadrupole strengths and reconstructed using ouralgorithm. (Right, top) Learning rate schedule for snapshot ensembling. (Right, bottom) Evolution of
the proposal distribution beam size in each phase space coordinate during training in synthetic case.
Dotted lines denote ground truth values. Vertical lines denote snapshot locations.
We characterize the confidence of our reconstruction using snapshot ensembling [27]. This encourages
the optimizer to find multiple possible solutions (if they exist). It is instructive to examine the evolutionof the proposal distribution during model training (shown in Fig. 4). The phase space componentsthat have the strongest impact on beam transport through the beamline as a function of quadrupole
strength converge quickly, whereas the ones that have little-to-no impact (e.g. the longitudinal
distribution characteristics) do not. In particular, we observe that each snapshot converges to beam
distributions which have slightly varying energy spreads, signifying uncertainty in that aspect of thebeam’s phase space, likely due to weak coupling between energy spread and transverse beam size
through chromatic aberrations in the quadrupole. The convergence of the proposal distribution thus
provides a useful indicator of which components of the phase space can be reliably reconstructed forarbitrary (possibly unknown) tomographic beamlines.
44Conclusion

Here we have shown the advances in accelerator modeling that can be achieved with differentiable
computations. When coupled with machine learning techniques these models can surpass limitations
faced by historically used methods for analyzing experimental data. Further investments in creating
fully differentiable accelerator physics simulations will continue to improve model detail and accuracy.

5

Impact Statement

This work is expected to have significant impact on the ability of accelerator diagnostics in the future.
As the work is limited within the scope of accelerator science we expect there to be no ethical aspects
or future societal impacts beyond the improvement of accelerator operations.

References[1] N. J. Sammut, H. Burkhardt, S. M. White, C. Giloux, and W. Venturini-Delsolaro, “Measurementand Effects of the Magnetic Hysteresis on the LHC Crossing Angle and Separation Bumps,”
Aug. 2008. Number: LHC-PROJECT-Report-1107.
[2] N. Huque, M. Abdelwhab, E. Daly, and Y. Pischalnikov, “Accelerated Life Testing of LCLS-II
Cavity Tuner Motor,” in 17th International Conference on RF Superconductivity, p. THPB062,
2015.
[3] D. A. Turner, G. Burt, and T. Junginger, “No DC superheating and increased surface pinning in
low temperature baked niobium,” tech. rep., Jan. 2022. ISSN: 2693-5015 Type: article.
[4] J. Snoek, H. Larochelle, and R. P. Adams, “Practical Bayesian Optimization of Machine
Learning Algorithms,” in Advances in Neural Information Processing Systems 25 (F. Pereira,C. J. C. Burges, L. Bottou, and K. Q. Weinberger, eds.), pp. 2951–2959, Curran Associates, Inc.,
2012.
[5] I. D. Mayergoyz and G. Friedman, “Generalized Preisach model of hysteresis,” IEEE Transactions on Magnetics, vol. 24, pp. 212–217, Jan. 1988. Conference Name: IEEE Transactions on
Magnetics.
[6] G. Bertotti and I. D. Mayergoyz, eds., The science of hysteresis. Amsterdam ; Boston: Academic,
1st ed ed., 2006.
[7] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin,N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani,
S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “PyTorch: An Imperative Style,High-Performance Deep Learning Library,” in Advances in Neural Information Processing
Systems 32 (H. Wallach, H. Larochelle, A. Beygelzimer, F. d. Alché-Buc, E. Fox, and R. Garnett,eds.), pp. 8024–8035, Curran Associates, Inc., 2019.
[8] C. E. Rasmussen and C. K. I. Williams, “Bayesian Regression and Gaussian processes,” Gaussian Processes for Machine Learning, p. Chapter 2, 2006. ISBN: 026218253X.
[9] J. R. Gardner, G. Pleiss, D. Bindel, K. Q. Weinberger, and A. G. Wilson, “GPyTorch: Blackbox
Matrix-Matrix Gaussian Process Inference with GPU Acceleration,” in Advances in Neural
Information Processing Systems, 2018.
[10] Y. Sun, M. Borland, G. I. Fystro, X. Huang, and H. Shang, “Recent Operational Experience
with Thermionic RF Guns at the APS,” in Proc. 12th Int. Particle Accelerator Conf. (IPAC’21),
pp. 3959–3962, JACoW Publishing, May 2021.
[11] R. Roussel, A. Edelen, D. Ratner, K. Dubey, J. Gonzalez-Aguilera, Y. Kim, and N. Kuklev,“Differentiable Preisach Modeling for Characterization and Optimization of Particle Accelerator
Systems with Hysteresis,” Physical Review Letters, vol. 128, p. 204801, May 2022. Publisher:
American Physical Society.
[12] C. McKee, P. O’Shea, and J. Madey, “Phase space tomography of relativistic electron beams,”
Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers,
Detectors and Associated Equipment, vol. 358, no. 1, pp. 264–267, 1995.
5

[13] S. Hancock, M. Lindroos, E. McIntosh, and M. Metcalf, “Tomographic measurements oflongitudinal phase space density,” Computer Physics Communications, vol. 118, no. 1, pp. 61–
70, 1999.
[14] D. Stratakis, R. A. Kishek, I. Haber, M. Walter, R. B. Fiorito, S. Bernal, J. Thangaraj, K. Tian,
C. Papadopoulos, M. Reiser, and P. G. O’Shea, “Phase space tomography of beams with extreme
space charge,” in 2007 IEEE Particle Accelerator Conference (PAC), pp. 2025–2029, 2007.
[15] V. Yakimenko, M. Babzien, I. Ben-Zvi, R. Malone, and X.-J. Wang, “Electron beam phase-spacemeasurement using a high-precision tomography technique,” Physical Review Special Topics Accelerators and Beams, vol. 6, p. 122801, Dec. 2003.
[16] M. Röhrs, C. Gerth, H. Schlarb, B. Schmidt, and P. Schmüser, “Time-resolved electron beamphase space tomography at a soft x-ray free-electron laser,” Physical Review Special Topics Accelerators and Beams, vol. 12, p. 050704, May 2009.
[17] M. Gordon, W. Li, M. Andorf, A. Bartnik, C. Duncan, M. Kaemingk, C. Pennington, I. Bazarov,
Y.-K. Kim, and J. Maxson, “Four-dimensional emittance measurements of ultrafast electron
diffraction optics corrected up to sextupole order,” Physical Review Accelerators and Beams,
vol. 25, no. 8, p. 084001, 2022.
[18] K. Hock and A. Wolski, “Tomographic reconstruction of the full 4d transverse phase space,”
Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers,
Detectors and Associated Equipment, vol. 726, pp. 8–16, 2013.
[19] J. C. Wong, A. Shishlo, A. Aleksandrov, Y. Liu, and C. Long, “4D transverse phase space
tomography of an operational hydrogen ion beam via noninvasive 2D measurements using laser
wires,” Physical Review Accelerators and Beams, vol. 25, p. 042801, Apr. 2022.
[20] M. Wang, Z. Wang, D. Wang, W. Liu, B. Wang, M. Wang, M. Qiu, X. Guan, X. Wang, W. Huang,and S. Zheng, “Four-dimensional phase space measurement using multiple two-dimensional
profiles,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators,
Spectrometers, Detectors and Associated Equipment, vol. 943, p. 162438, Nov. 2019.
[21] B. Hermann, V. A. Guzenko, O. R. Hürzeler, A. Kirchner, G. L. Orlandi, E. Prat, and R. Ischebeck, “Electron beam transverse phase space tomography using nanofabricated wire scanners
with submicrometer resolution,” Physical Review Accelerators and Beams, vol. 24, p. 022802,
Feb. 2021.
[22] A. Scheinker, F. Cropp, S. Paiagua, and D. Filippetto, “An adaptive approach to machinelearning for compact particle accelerators,” Scientific reports, vol. 11, no. 1, pp. 1–11, 2021.
[23] R. Roussel, A. Edelen, C. Mayes, D. Ratner, J. P. Gonzalez-Aguilera, S. Kim, E. Wisniewski,
and J. Power, “Phase space reconstruction from accelerator beam measurements using neuralnetworks and differentiable simulations,” arXiv preprint arXiv:2209.04505, 2022.
[24] I. Kobyzev, S. Prince, and M. Brubaker, “Normalizing flows: An introduction and review ofcurrent methods, arxiv e-prints,” arXiv preprint arXiv:1908.09257, 2019.
[25] K. M. Hock and M. G. Ibison, “A study of the maximum entropy technique for phase space
tomography,” Journal of Instrumentation, vol. 8, pp. P02003–P02003, feb 2013.
[26] J. Lawson, R. Gluckstern, and P. M. Lapostolle, “Emittance, entropy and information,” Part.
Accel., vol. 5, pp. 61–65, 1973.
[27] G. Huang, Y. Li, G. Pleiss, Z. Liu, J. E. Hopcroft, and K. Q. Weinberger, “Snapshot ensembles:
Train 1, get m for free,” arXiv preprint arXiv:1704.00109, 2017.

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [No](c) Did you discuss any potential negative societal impacts of your work? [N/A](d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]6

2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A](b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [No] Implementations
included in cited work(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they
were chosen)? [N/A](c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [N/A](d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [N/A]
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [N/A](b) Did you mention the license of the assets? [N/A](c) Did you include any new assets either in the supplemental material or as a URL? [N/A](d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]7