Symplectic Modeling of Beam Loading in Electromagnetic Cavities

Simulating beam loading in radiofrequency accelerating structures is critical for understanding higher-order mode effects on beam dynamics, such as beam break-up instability in energy recovery linacs. Full wave simulations of beam loading in radiofrequency structures are computationally expensive, while reduced models can ignore essential physics and can be difficult to generalize. We present a self-consistent algorithm derived from the least-action principle which can model an arbitrary number of cavity eigenmodes and with a generic beam distribution.


I. INTRODUCTION
The interactions of charged-particle beams and electromagnetic fields drive a variety of applications and phenomena in accelerator systems. They range from klystrons to beam-loading in radiofrequency (RF) cavities to highorder-mode instabilities in energy-recovery linacs. Accurate, self-consistent modeling of these interactions is an important and challenging problem for many of these applications.
Self-consistent simulations using approaches such as finite-difference time-domain (FDTD) particle-in-cell (PIC) algorithms [1] are computationally demanding and suffer from a variety of numerical artifacts-for example, the numericalČerenkov instability [2], or dispersive errors induced by the mesh [3,4]-that can muddy the results. Detailed simulations of an electron beam passing through an RF cavity must resolve the current (i.e. the electron beam), and hence require cell dimensions that are small compared to the beam. Resolving the beam across the entire complex geometry can therefore require billions of cells, while the actual source terms are isolated to a mere few thousands of cells. As a consequence, load balancing becomes a serious concern for multi-core simulations.
The discrepancy of spatial scales is not the only challenge faced by FDTD-PIC algorithms. One must also initialize the fields on the numerical mesh in a manner that satisfies the appropriate numerical dispersion relation [3,5,6]. In addition, the presence of multiple bunches means that one must address the discrepancy of temporal scales that effectively prevents using such algorithms to study the effect of high-order modes on beam dynamics.
To bridge multiple scales, one often resorts to reduced models, as, for example, when studying beam loading (matbbu [7]) or beam instability (bi [8]). Such models typically treat the RF cavity as a thin lens, using reduced forms of the Shockley-Ramo theorem [9,10] to compute the energy transfer; and they can advance the field phases analytically once the beam passes. Reduced models are computationally much more efficient than full simulations. They also simplify the diagnostics, because the energy in each cavity mode is a dynamical variable and requires no * dabell@radiasoft.net; www.radiasoft.net additional computation to be extracted from the simulation. However, they can fail to demonstrate complicated phase-space structures that result from beam-beam disruption in a collider (as in the eRHIC energy recovery linac-ring design for an electron-ion collider [11,12]) or from a large energy spread (as results from advanced tapering schemes for free-electron lasers). The approximations made in reduced models frequently do not fail gracefully, and they are difficult to expand to next-leading-order, which can limit their versatility.
What we require is an algorithm that has the simplicity and physically intuitive feel of a reduced model, while being extensible for self-consistent simulations. In this paper, we present such an algorithm based on a symplectic map approach to electromagnetic particle-in-mode algorithms [13]. For this approach we use the cavity eigenmodes as our orthonormal basis for the electromagnetic field. We consider only the coupling of j z to A z , neglecting the transverse currents. And the modes themselves we compute using either an analytic model or interpolation of numerical data [14]. This approach leads to a fast, extensible model for beam loading with arbitrary numbers of modes and cavity geometries.

II. THE RATIONALE FOR A SPECTRAL TIME-BASED ALGORITHM
Self-consistent updates of the Maxwell equations require that the fields obey the boundary conditions. For FDTD-type algorithms-the canonical example being the staggered-Yee scheme [15]-this is not difficult: the basis used to decompose the fields is local in space, and it is easy to compute reflections off boundaries, transient effects, etc. FDTD algorithms are popular also because they can handle complex boundaries, and the electromagnetic part of the algorithm requires only local information to update the fields. However, these algorithms must resolve the smallest features in a simulation, typically the beam itself, and they have inherent dispersive errors that make accurate simulation of beam loading a challenge. Artificial instabilities such as numericalČerenkov add to the challenge, thus overall making the FDTD approach unsuitable for studying beam loading in most structures.
One possible solution is to use a spectral algorithm that does a global decomposition of the fields as a superpo-sition of the cavity eigenmodes. This is similar in spirit to the Condon method for computing wakefields [16][17][18]. Each mode frequency is known exactly 1 , meaning that in the absence of a source term it is possible to evolve the fields exactly. Adding the source term-the beam-is a convolution integral of the source currents and charges with each field eigenmode. The beam itself takes up a small volume of the simulation domain. Using the field eigenmodes removes the longest length scale from the simulation, dramatically reducing the computational resources required. It also greatly simplifies source deposition in an electromagnetic simulation, as compared to deposition for FDTD, which can encounter load balancing issues that impede parallel scaling 2 . This can make the spectral approach competitive for performance with FDTD algorithms. Because beam loading requires exact frequency information and field maps to compute properly, and the sources are located on a small fraction of the simulation domain which causes severe load balancing issues for FDTD particle-in-cell simulations, spectral algorithms are ideal for self-consistent beam loading simulations.
Conventional accelerator tracking codes almost universally use distance along the design orbit of a ring or linac as the independent variable (s-based tracking). When including collective effects, a code must address the issue of simultaneity in the Poisson equation: that leads to codes either adopting time as the independent variable (t-based tracking), or using various tricks to make the particle times simultaneous. Single particle integration through standing wave structures, such as radiofrequency cavities, can be accomplished so long as the self-consistent evolution of the fields is neglected.
In a waveguide type structure, which has translational symmetry in the longitudinal direction, s-or t-based tracking can be used. For systems with s-varying transverse geometries, such as RF cavities or traveling wave tubes, we must use t-based tracking. It is easiest to understand this by looking at how an ultra-relativistic beam (β ≈ 1) would affect an s-versus t-based spectral representation. The fields from the bunch radiate purely transversely, before encountering the cavity surface. If that surface has any tilt to it, the fields reflect off the cavity surface and back into the cavity. In an s-based approach, this would produce a local-in-s increase in the individual cavity eigenmode strengths on the surface, creating surface currents. These surface currents are not self-consistent with any single eigenmode, so the fields would re-radiate. An s-based eigenmode representation has a scattering term between multiple eigenmodes. This is not the case for a t-based approach, which changes the eigenmode amplitudes and their associated surface currents by a multiplicative constant, but which are still self-consistent for a single eigenmode.
Simply, an s-based algorithm with s-varying geometry has scattering terms between the eigenmodes, while a tbased algorithm would not. For the problem of a standing wave rf cavity, this means that a t-based algorithm is simpler. For a waveguide with no longitudinal variation in the boundaries, either s-or t-based algorithms can be used.

A. Beam Low Lagrangian
We begin with the Low Lagrangian [19], devised by Low in 1958 to provide a variational formulation for describing a non-relativistic ionized gas in an electromagnetic field. It has been used, for example, to study nonlinear waves in plasmas [20]. The extension to a relativistic form is straightforward, and has been used, in contexts related to the present work, to develop other algorithms for plasma simulations [21][22][23]. In Gaussian c.g.s. units, which we use throughout this work, it has the form Here the Boltzmann function ψ(x 0 , v 0 ) describes the phase-space distribution of particles at location (x 0 , v 0 ); 1 Or to within numerical tolerances, if the fields are solved numerically using a frequency domain solver. 2 It is worth noting that because spectral algorithms are global instead of local algorithms, they require an MPI AllReduce to consolidate the source terms. This prevents spectral algorithms from scaling indefinitely, and may require specialized algorithms which reduce the need for many-core simulations at the expense of generality.
φ and A denote the electromagnetic scalar and vector potentials; and τ = ct denotes our independent variable. The first integral describes the particles, including both their kinetic energy and their interaction with the electromagnetic field. The second integral describes oscillations of the electromagnetic field.
If we neglect the beam space charge, and consider only the cavity eigenmodes, then φ = 0. Furthermore, our beam has | dx ⊥ / dτ | dz/ dτ ∼ 1, while |A ⊥ | ∼ A z . We therefore choose to neglect the transverse coupling terms. We do not have to do this: all the algorithms and computational results described in this paper can be (and is some cases have been) done without making this simplification; but doing so leaves us with what we call the beam electromagnetic Low Lagrangian: Our approximation-in essence neglecting j ⊥ -is made in the interest of speed, as it dramatically reduces the number of required source depositions, which dominates the computation time for self-consistent simulations. For applications in which there is, for example, significant gyrotron motion, one may need to retain the transverse coupling terms. We assume the cavity eigenmodes are known and form a complete orthonormal basis, so that we may decompose the vector potential in these modes: where the f denote spatial eigenmodes, and a the corresponding mode amplitudes. Plugging this form into the Low Lagrangian, we obtain where overdot denotes differentiation with respect to the proper time τ . In going from (2) to (4), we have used the orthogonality of electric and magnetic field eigenmodes to eliminate the cross-terms in the second integral. As a convenience, we also define the mode inductance 1/L = 1 4π dx | ∇× f | 2 , and the mode capacitance To trace the particles, we introduce macroparticles by decomposing the phase-space density in discrete shapes in the usual manner [13,[22][23][24]: Here w j denote the macroparticle weights, δ the Dirac delta function, and Λ the normalized particle shape functions (so that dx Λ = 1). This decomposition transforms the Lagrangian (4) into a discrete set of coupled macroparticle-electromagnetic-mode Lagrangians: where represents the coupling of particle shape to field shape. Each of the individual x (j) and a are dynamical variables in this particle-in-mode discrete Lagrangian. One approach to numerically integrating the resulting equations of motion revolves around using time-discrete Lagrangians [25]. This approach, however, typically yields an implicit algorithm, which is substantially slower than an explicit version. We therefore opt to develop explicit symplectic maps using a canonical formalism.

B. Particle-Field Hamiltonian and a Split Operator Approach
We can compute the Hamiltonian corresponding to the Lagrangian from the canonical coördinates and conjugate momenta for the particles and fields: The particle conjugate momenta are the usual, but with the A ⊥ components neglected. The field conjugate momentum includes the mode capacitance, which acts like an effective mass for the mode. The Legendre transform to compute the coupled particle-field Hamiltonian is Carrying out this Legendre transform for these variables yields the particle-field Hamiltonian Here H pc denotes the particle-coupling Hamiltonian, which describes both the particle dynamics and the particle-field coupling; and H f denotes the field Hamiltonian, which describes the harmonic oscillation of the independent field eigenmodes.
The mode quantities C and L are determined only up to within a normalization constant. This means that the individual definitions of Q and P will depend on the normalization convention used for f . It will therefore prove useful to make a canonical transformation that combines the two into a single scale-invariant quantity, which would be an intrinsic quantity for a given cavity eigenmode.
By defining the canonically conjugate variables we transform the field and particle-coupling Hamiltonians into and All quantities in these Hamiltonians-including Q , P , and F / √ C -are independent of our choice of normal-ization for f . The invariant combination of the mode capacitance and inductance is which yields the eigenmode frequency, Ω .
C. Splitting the Hamiltonian The total particle-field Hamiltonian has no explicit time dependence, so the symplectic map for going from time τ to time τ + h is simply Here we use the colon notation introduced by Dragt [26][27][28] to denote that argument of the exponential is a Poisson-bracket Lie operator. The full map (14) is difficult to compute, but a symmetric splitting yields a map that is both straightforward to compute and second-order accurate in the step-size h: where and all maps are independent of the initial time τ . One can evaluate M f exactly and M pc to second order in h, leading to an overall second-order accurate symplectic integrator for both fields and particles.

Field Map
The field map M f , is simply the harmonic oscillator, and it acts only on the field phase-space coördinates: As is usually the case with second-order splitting, the two half-step field maps can, for simplicity, be combined into a single full-step map-so long as one retains the half-step maps at the ends of the simulation.

Particle-Coupling Map
The particle-coupling map M pc is not immediately integrable, but it can be split using a method described by Wu, Forest, and Robin [29] and exploited by Webb et al. [13] for a cylindrical electromagnetic algorithm. This method involves tracking each particle with respect to its own proper time, which allows us to re-write the Hamiltonian as an effectively non-relativistic Hamiltonian given by where γ (j) denotes the Lorentz factor of the j th macroparticle: Because the vector potential has no explicit time dependence, γ (j) is a constant of the motion for this Hamiltonian, and the proper time of each macroparticle can be mapped to the lab time simply by multiplying by γ (j) . Once written in the form (17), the Hamiltonian can be split and integrated using a pair of half-drift maps for the transverse coördinates, and another map for the longitudinal coördinate: This result is again second-order accurate in h and hence preserves the overall order of the integrator.
The leading and trailing half-drifts are particularly simple: We can now apply the technique in [29] to deal with the vector potential. Lie transformations obey an important similarity transformation property: e :f : e :g : e −:f : = e :e :f : g : .
This property allows us to rewrite the central map in (19) as the product of three maps: Moreover, each of these maps-the drift D z and the trans-formation A z -can be evaluated exactly: Because the algorithm is based on explicit symplectic maps, it is possible to implement each of these maps as a single function, and the update sequence is just a series of function calls. Once each of these maps is implemented, it is straightforward to begin simulations.
Before going on to discuss our numerical results, we should clarify a point that may seem mysterious: namely, how does the Lorentz γ factor change? The method given by Wu, Forest, and Robin is described in the context of magnetostatic systems. Moreover, the absence of (explicit) time dependence in the vector potential that appears in the Hamiltonian H pc of (17) is consistent with a magnetostatic system, which cannot change the Lorentz factor. And indeed the map M pc derived from the Hamiltonian H pc does not do so.
The resolve the apparent mystery, recall that H pc denotes just the particle-coupling term in the full Hamiltonian: the remaining part is the field Hamiltonian, given in (12a). When stepping through the full map (15) to update the system for a single time step, one must apply the field map M f of (16), which rotates the field coördinates, Q , and momenta, P ; and it is effectively this operation that updates-see (18)-the Lorentz γ factor.

IV. NUMERICAL RESULTS
To test the algorithm described in this paper, we considered a number of configurations. In all tests, we considered a rectangular pillbox cavity with a handful of test macroparticles. We list the specific parameters in table I. In all simulations, we considered four cavity modes: an accelerating mode, two transverse dipole modes, and a single transverse quadrupole mode.
An artificially high bunch charge was selected to demonstrate the field energy loss-otherwise the total energy of the bunch is so much smaller than the total energy in the fields that the variations are not easy to discern.
In the first test, we considered a configuration with a single macroparticle representing the entire bunch, accelerated from γ = 50 to γ = 51. The particle is on-axis, so no higher order modes should be excited. The result is shown in figure 1. As can be seen, the energy conservation is excellent, and none of the higher order modes are excited. Eigenmode Energies (one macro-particle on axis) TM110 TM210 TM120 TM220 FIG. 1. Energy gain and eigenmode energy decomposition for the case of one macro-particle on axis.
To test the excitation of higher order modes, we split the beam into four macroparticles-one on-axis; one offset in x and another offset in y to excite dipole modes; and one offset along the x-y diagonal to excite both dipole and quadrupole modes. As shown in figure 2, the total energy is conserved to high precision, and all the modes are excited, although the fundamental mode sees the most variation. Because this is a symplectic integrator, we expect the total energy to be well-behaved for arbitrarily long time. To test this, we considered a much longer rectangular pillbox, 125 cm long, with the same transverse eigenmodes. As shown in figure 3, the algorithm retains stable energy behavior for longer time scales. This captures oscillatory behavior as the beam accelerates and decelerates, and the particles slip in phase for the higher order modes.
When simulating considerably longer cavities, ∼ 180000 cm in length, we detected a slow spurious decay in the total energy. Over this length, the error was ∼ 10 −6 % total, as shown in figure 4. We have identified the origin of this spurious decay as numerical roundoff error in the field rotation matrix: for our implementation,the floating-point value of the determinant was very slightly less than 1. Over many time steps this can lead, depending on the roundoff error, to spurious cooling or heating of the system; but for a single bunch pass this is almost undetectable. Furthermore, when advancing the fields in between bunch passes, it suffices to perform a single matrix multiplication to advance the eigenmodes. Doing this reduces the roundoff-induced energy change in multi-bunch and multi-pass simulations performed when studying, for example, energy-recovery linacs.

V. ALGORITHM PERFORMANCE
To demonstrate the second-order behavior of the algorithm, we performed a simple test with a rectangular 10 cm × 10 cm × 20 cm cavity instantiated with three modes. We propagated four particles through the cavity for a single step and evaluated the change in the system Hamiltonian as the step size was varied. The resulting behavior is shown in figure 5. For decreasing step size, h, the error falls as h 3 , indicative of a second-order algorithm. For step sizes below 1 mm, the single-step error quickly approaches machine precision, on the order of 10 −14 relative error.
Initial benchmarks for the algorithm's performance demonstrate the scalability to many hundreds of modes and thousands of particles. We performed 100 steps of fixed step size on a single processor with varying macroparticle number and mode number to illustrate the single-step evaluation time per-macroparticle-per-mode; the results  are displayed in figure 6. Although Python-induced overhead associated with array creation hinders performance for small samples, the step times settle near ∆t ∼ 0.5 µs per-particle-per-mode. When particle and mode numbers cause array size to exceed cache size, a small drop in performance occurs, as can be seen in to the top right of the figure.

VI. CONCLUSIONS
We have presented a new approach to simulating beam loading in electromagnetic cavities. The approach is based on a particle-in-mode symplectic map approach to electromagnetic charged-particle simulations. The algorithm is suitable for modeling a variety of vacuum electronic devices, including radiofrequency cavities, traveling wave FIG. 6. Timings for a single step weighted per particle and per mode are shown for upwards of 256 modes and 12800 particles. For large systems, we find ∆t ∼ 0.5 µs per-particle-per-mode.
tubes, and klystrons. For applications in energy recovery linacs, the beam remains relativistic at all times, effectively suppressing space-charge forces. We therefore chose to ignore it in our current implementation. But it is possible to include space-charge within the framework of our algorithm, and we mention two possible approaches: One possibility is to stick with the Weyl gauge (φ = 0) and introduce additional high-order modes to represent the space-charge fields. When the bunch is very short, the required modes will correspond to very-high-frequency, essentially freespace, modes that are broadly separated from the cavity modes. This approach will work well if the space-charge variation is predominantly longitudinal, or if one includes the j ⊥ terms we neglected in transforming from the Lagrangian of (1) to that of (2). The other possibility is to use a spectral algorithm described in a paper by one of the present authors [24].
Another issue we did not address in our current implementation is the use of realistic rf cavities. To model realistic cavity eigenmodes, one may use the technology of generalized gradients described in reference [14]. There the fields are computed on the basis of generalized gradients, which one computes as integrals over a surface enclosing the relevant volume. This surface integration acts to damp any imprecision in the simulated fields, and it means that the field evaluations-including the derivatives-needed for the maps in (23) can be computed to high accuracy.
The algorithm presented here avoids the various pitfalls of both the finite-difference time-domain electromagnetic particle-in-cell and the reduced-model approaches that are typically applied to this problem.
Because the algorithm is directly spectral, there are no grid aliasing artifacts that can result in the numerical instabilities and dispersive errors that tend to infect the full electromagnetic PIC simulations. By using a modal decomposition, it is straightforward to identify, directly from the dynamical data, which cavity eigenmodes are excited and how they evolve-there is no need for post-processing of the simulation data to extract the eigenmodes. The approach also solves the load balancing problem, as the field data is global, and the particles can be evenly distributed across multiple processors. The approach is also much faster than the fully self-consistent approach: our early simulations of the single particle traversing the pillbox cavity, including the time to generate the figures, required less than a second using an implementation that performed symbolic mathematical computations and is far from optimized. See figure 6 for performance data based on a more sophisticated implementation.
The symplectic particle-in-mode (SymPIM) algorithm described in this paper in many ways resembles the reduced models, but it is distinct in a number of important ways. Because the cavity is not treated as a thin element, it is possible to model the complex phase-space of realistic beams, such as those with large energy spreads, or those disrupted by beam-beam collisions or various wakefield instabilities. The SymPIM algorithm thus lies between the fully self-consistent approach-which includes all resolvable cavity modes at great computational expense-and the reduced-model approach-which neglects many details of the beam distribution and how it evolves within the cavity. This makes the SymPIM algorithm suitable for determining the validity of reduced models, as well as for modeling those cases where the reduced models do not apply.

VII. ACKNOWLEDGEMENTS
This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award Number DE-SC0015211.