Equilibrium ion distribution in the presence of clearing electrodes and its influence on electron dynamics

Here we compute the ion distribution produced by an electron beam when ion-clearing electrodes are installed. This ion density is established as an equilibrium between gas ionization and ion clearing. The transverse ion distributions are shown to strongly peak in the beam's center, producing very nonlinear forces on the electron beam. We will analyze perturbations to the beam properties by these nonlinear fields. To obtain reasonable simulation speeds, we develop fast algorithms that take advantage of adiabatic invariants and scaling properties of Maxwell's equations and the Lorentz force. Our results are very relevant for high current Energy Recovery Linacs, where ions are produced relatively quickly, and where clearing gaps in the electron beam cannot easily be used for ion elimination. The examples in this paper therefore use parameters of the Cornell Energy Recovery Linac project. For simplicity we only consider the case of a circular electron beam of changing diameter. However, we parameterize this model to approximate non-round beams well. We find suitable places for clearing electrodes and compute the equilibrium ion density and its effect on electron-emittance growth and halo development. We find that it is not sufficient to place clearing electrodes only at the minimum of the electron beam potential where ions are accumulated.


I. INTRODUCTION
Several processes can contribute to the production of ions in the vicinity of accelerated electron beams [1]. These positively charged ions are attracted to the negative beam. Anti-proton accumulation systems have suffered under the field from accumulated ions [2,3], and so have electron beams [4,5].
Ions can damage the beam in various ways. Firstly, they produce a focusing field that changes in a very nonlinear way with an electron's distance from the beam's center. The resulting nonlinear motion can increase the beam emittance, can lead to particle loss, and can produce a halo around the beam. Secondly, the ion distribution oscillates within the electron beam, while the electrons are attracted to the ion distribution. This system of coupled oscillators can become unstable [6,7,8,9] and lead to large transverse beam oscillation amplitudes or to an increase of the apparent transverse beam size.
It is therefore important to reduce the density of ions in the vicinity of the beam to a tolerable amount. Different accelerators achieve this differently. Storage rings typically use ion-clearing gaps. These are short gaps in the filling pattern that lead to an absence of focusing forces for the ions every time this gap has traveled around the ring. When the gap length and frequency are chosen suitably, the ions get over-focused which lets the ions oscillate with increasing amplitude until they have moved outside the beam region. The length of the beam-filled region is typically a few microseconds long for a large accelerator, whereas the gap is typically shorter than one microsecond. In pulsed linacs, the gaps are often much longer and allow ions to drift out of the beam region.
In rings with coasting beams, there is no ion-clearing gap, and obviously the beam cannot be turned off regularly. Similarly, in Energy Recovery Linacs (ERLs) [1] where the beam's energy is dumped in RF cavities and is immediately used to accelerate new electrons, one cannot easily turn off the beam (because this would interrupt the ERL process) and one can also not easily introduce short gaps in the beam (because this would disrupt the gun or the linac that injects large currents into the ERL). In both of these cases, ion-clearing electrodes may have to be used [10,11].
The electron beam diameter varies along the accelerator, and this variation produces longitudinal forces, guiding ions to a location where the electron density is relatively large, typically close to the waists of the electron beam. Clearing electrodes are placed along the beam-line at such places of ion accumulation.
Because these longitudinal forces are relatively weak, it typically takes a few milliseconds for an ion to move from the place where it is produced to an ion-clearing electrode a few meters down the beam-line. During this time, new ions are created, leading to a remnant equilibrium ion density that establishes itself in the presence of clearing electrodes. Creating a number of ions per length as large as the number of electrons per length typically takes only a few seconds, even for very good vacuum in the nTorr range [1]. Because the motion to clearing electrodes that are spaced many meters apart takes several milliseconds, such electrodes tend to produce a linear ion density that is in the order of about one part in a thousand of the linear charge density in the beam.
Computing an equilibrium density is typically very time consuming. Here we show how to use scaling properties of the Maxwell equations and the Lorentz force, as well as the adiabatic invariance of the action integral, to compute the equilibrium ion density very efficiently.
We subsequently apply this technique to analyze the Cornell ERL project [12]. We investigate whether the remnant ion density after placement of clearing electrodes can lead to any of the discussed damages to the electron beam, so that additional clearing electrodes would have to be installed. We assume the parameters listed in Tab. I, where the ionization cross-section is taken from [13] for 5GeV. To be conservative we used the high energy cross section for all energies, even though for lower energies the cross-section is up to 40% smaller. The gas density of 3·10 13 m −3 corresponds to a pressure of 1nTorr for room temperature sections of the accelerator. The unnormalized emittances are obtained with the relativistic γ by ǫ x = ε nx /γ.

II. COMPUTATION OF THE ION EQUILIBRIUM
A. The 3D beam-beam force For a rotationally symmetric electron beam, the density is given by where λ is the linear particle density. With Gauss' law, the transverse velocity kick on a singly charged, non relativistic ion becomes where σ is the rms width of the electron-beam profile, r is the distance from the beam centerline, N the number of electrons in the bunch, r p the classical proton radius and A the atomic weight of the ion. For non-round Gaussian electron beams, the transverse kick of the ion can be evaluated analytically in terms of complex error functions [14]. If the size of the electron beam changes as a function of s, the ion will also experience a longitudinal force from the passing bunch as shown in Fig. 1. Linearizing in the slopes of the electron trajectories, this longitudinal force can be expressed as a function of transverse force components and beam parameters [15]: The velocity changes that are produced by one bunch passage can be seen in Fig. 2. Here we used the values of Tab. I for top energy and Twiss parameters of β x = β y = 2m, α x = α y = 1.
In our restriction to rotationally symmetric beams, we have a dispersion of η = 0, α = α x = α y , ε = ε x = ε y , and the longitudinal kick becomes For the beam parameters we studied, the longitudinal kick is typically five orders of magnitude weaker than the transverse kick, as shown in Fig. 2. In the special case of rotational symmetry, the longitudinal kick on the ion will be in the same direction for all positions. This is not true in the general case; if α x and α y have different signs, the kick direction can depend on the distance from the centerline [15].

B. Adiabatic invariants of ion motion
In a brute force simulation, the state of each ion is characterized by its x and s position, as well as by the velocities in both directions. Every kick has to be calculated separately from the beam-beam force, which is assumed to be a function of the x and s position. Such a simulation would be very inefficient because there are two different timescales: A large number of kicks is required to resolve the sharply focused transverse ion motion, but the longitudinal motion hardly changes in a single oscillation. It is therefore not possible to simulate the motion of a realistic number of ions in this way.
The period of transverse oscillations in the beam's potential [16,17,18] depends on the average beam density, and for typical high-brightness beam parameters as those of Tab. I it is of the order of several to tens of bunch crossings. During one transverse oscillation the longitudinal position s of the ion changes typically only by a fraction of a millimeter even at the largest possible ion speeds.
Ions thus oscillate in a potential that slowly changes over many periods. Because of this, an ion's action integral over one oscillation period is an adiabatic invariant of the motion. This can be used to drastically speed up the simulation. In this improved simulation, the state of the ion is now described by its longitudinal position s, its longitudinal speed v s and its action integral J. For each position s, J is a measure of the transverse oscillation amplitude a, because the action increases with a. This can most easily be seen in (x, p x ) phase space where the graph (x(t), p x (t)) for any oscillation is a closed curve, and J is the enclosed area. A closed curve that starts with a larger oscillation amplitude has to completely enclose one that starts with a smaller amplitude because different phase-space trajectories cannot cross. Solving the motion in (J, s) rather than in (x, s) coordinates has two advantages: (1) The degrees of freedom are reduced from 2 to 1. (2) Because ∆v s changes much slower than ∆v x , the integration steps can be vastly increased, typically by about a factor of 10000. If the density at s is to be computed, it is not enough to know the action at s, but one oscillation is sufficient to compute the density contribution of particles with action J.
The longitudinal acceleration averaged over one oscillation is given by the time average over individual kickṡ Calculating this time average again and again for each ion is computationally demanding and wasteful. Instead, we pre-compute a table of possible values of this time average for many beam sizes σ and amplitudes a of the oscillation.
The program was further accelerated by noting that the time averaged kick as a function of J can be calculated for a typical standard beam size and then rescaled for regions with other beam sizes.
To derive this scaling behavior, we use a general scaling theorem for Maxwell's equations and the Lorentz force. If E( x, t), B( x, t), ρ( x, t), and j( x, t) satisfy Maxwell's equation, and a charged particle with mass m and trajectory r(t) satisfies the Lorentz force equation, then the following scaled quantities also satisfy these equations, The effect of this scaling transformation is to reduce spatial distances and time intervals by a factor α, so that velocities remain unchanged. Also, the charge and current densities are increased by a factor α 3 , so that the total charge in a scaled volume remains constant.
We now want to find out how the spatial derivative of the transverse kicks changes when the transverse scale of the beam is reduced by a factor α. If the rotationally symmetric beam produces a field E(x, s), the 3-D scaled field would be α 2 E(αx, αs). However, only the transverse scale changes from one beam width to another, and the longitudinal scale stays the same, signified by the unchanged line density. The charge within a 3-D scaled volume is therefore smaller by 1/α, and the field after transverse scaling is therefore E 1/α (x, s) = α E(αx, s).
The transverse scaling transformation therefore gives ¿From this, we see that a rescaling of the two relevant length scales σ and a has the effect This is true for both derivatives, so that the time average of their sum multiplied by σ 2 is only a function of a/σ, Our goal is to express this time average as a function of the beam size σ and the action integral J, so we also have to find the functional dependence of J on a and σ.
As noted above, the electric field at a location where the beam size is reduced by a factor α, is scaled in the transverse to αE(αx, s). To obtain scaled transverse ion motion, we do not change the mass of the ions, i.e. m s = m.
The transverse velocity is then not changed by transverse scaling, and according to Eq. (5), the action integral is reduced by a factor α. From this we conclude Combining this result with the above analysis, we finally find To speed up the program, a lookup table for the function f is therefore computed for a standard beam size σ and a range of action integrals J. The longitudinal acceleration of the ions is then given bẏ C. Fast ion propagation to equilibrium The toy model for our ion simulations consists of a field-free beam region of variable length L, so that the β-function is given by with waist at s = 0. We then simulate the creation of ions in the beampipe through ionization processes. The ionization rate in a given volume V is given by where ρ e is the local particle density of the electrons as given by the Gaussian beam profile ρ g . By integrating over the transverse directions, one obtains an equation for the number of ions created in a length L of the beam ¿From this it is clear that ions are created with equal rate at all locations along the beam. Under the influence of the longitudinal beam force, the ions will then slowly propagate to the point with minimal potential.
In our simulation, we assume that clearing electrodes have been placed at the waist, so that ions crossing the s = 0 line are removed. An equilibrium situation is reached when an equal number of ions are produced and removed during a given time interval. The time needed to As a first example, we show a L=20m long part of the beam with a waist of β * = 1m at its center. The longitudinal equilibrium ion distribution for 10m up to the waist is shown in Fig. 3. There is a sharp increase in the ion density near the electrode at s = 0m that is caused by two effects. The first is that ions which have been created at all points of the beam pass through this region and contribute to the density. The second effect is that the α-function becomes smaller near the waist, so that ions which are created there experience a weaker longitudinal force and move more slowly to the electrode.
The transverse ion distribution can be characterized by its standard deviation, which is plotted in Fig. 4. Two features of this graph can be explained easily. (a) The rms width of the ions at the start of the section (at s = −10m) is smaller than that of the electron beam, even though ions are created in proportion to the Gaussian electron distribution. The reason is that ions start to oscillate within the electron beam after their creation and therefore contribute to the density around zero, no matter where they were created. (b) The ion distribution contracts less than the electron beam. This can be explained by the adiabatic invariance and the scaling of the action integral. Ions that are created with an action J in regions where the electron beam is wide will still have that action by the time they reach a part of the beam that is narrower by a factor α. At that position, their oscillation amplitude will have increased relative to the beam size, because the action for a proportionally reduced oscillation amplitude would be reduced by a factor α according to Eq. (9).
The density profile of the transverse ion distribution at s = ±1 near the clearing electrode is shown in Fig. 5. Near the beam axis, the density diverges approximately as 1/r, so that it is not appropriate to fit the ion distribution to a Gaussian. For a rotationally symmetric particle density, the number of ions created in a length ds and in a radius element dr per time iṡ n(r, s) drds = cσ col ρ g ρ e (r, s) 2πr drds. (15) Subsequent to the creation of an ion at (r 0 , s 0 ), it travels along the beamline to s = S(t; r 0 , s 0 ) and oscillates through the electron beam to R(s; r 0 , s 0 ). Ions created between s m and s can contribute to the distribution at s, and a ring element drds with radius r therefore contains dn ions, dn(r, s) = The number of particles per radius interval dn drds (r) is clearly not zero for r = 0, and the density close to zero therefore has a singularity of the form Because of repulsive forces between the ions, the density at the origin cannot diverge, and we therefore estimate the radius r m below which these forces become larger than the forces from the electron beam. It will turn out that this radius is a very small fraction of the electron beam's width, so that deviation from the 1/r density can be neglected.
We assume that the region with radius r i in which the ion density can be approximated as ρ ion = A r contains most of the charge. At some radius, the density has to fall off faster than 1/r to integrate to a finite value. Our assumption therefore leads to an overestimate of the ion force. Integrating ρ ion up to r i leads to the ion line density 2πAr i = λ ion .
The region that contains most of the charge has a radius r i that is usually smaller than the rms width of the electron beam. Only if ions were created in a very wide region of the beam, and it were to be focused down to an extremely small radius, could the ion width be much larger than the electron width. However, such extreme focusing does not occur in most beams, especially not in the Cornell ERL.
Gauss' law within r i leads to the constant radial electric field The radial electric field of the electron beam has already been used in Eq. (2) and for radii much smaller than σ simplifies to E e = λr 4πǫ0σ 2 . We can now find the radius R for which the ion field is as large as the electron field, R = σ λion λ 2σ ri . As long as the number of ions per length is much smaller than the number of electrons per length, the 1/r scaling of the ion density is therefore correct and not altered by the ion forces in by far the largest part of the beam.

A. Electron propagation
Because of the 1/r ion density characteristics, the electric field of the ion distribution approaches a constant for small distances, as specified in Eq. (18). This is shown in Fig. 6. To investigate the strength of the resulting beam distortion, we simulated the propagation of a representative set of beam electrons through the ion field. Our simulation was set up to include a piece of the beampipe of length L, with the clearing electrode and the waist of the beam placed in the center.
The initial electron distribution is determined by random betatron amplitudes J x and phases φ x with the density distribution We also investigated the effect of the beam transversing the region of length L several times, which could lead to the build-up of damaging resonances. For this study, the linear optics has to be made periodic after crossing the field-free region. This can be accomplished by including a thin half-quad with the transfer matrix before and after each run through the ion field. The α-function is α − = 0 before and α + = β kl/2 after the thin quadrupole. Simultaneously, this α-function must correspond to that of a drift, i.e. α + = L 2β * . Therefore, The simulation of the electron motion is implemented with a dynamical time step proportional to the distance from the beam axis, so that the effects of the strong field near the centerline are accounted for accurately.
As a worst case scenario, we simulated two drift regions with 200m and 100m length and a beam waist of β * = 100m at the center in both cases. The final phase-space distributions with and without the ion field are shown in Figs. 7 and 8.
While passing through the ion field once, the emittance increases from 30.6pm to 53.1pm in the 100m field and to 50.7pm in the 200m field, as can be seen in Fig. 9. The The previous two examples show that sections of only a few times 10 meters between clearing electrodes can produce intolerable emittance growth if the beta function is large, because in that case the beam is not very divergent and the longitudinal force that clears ions is consequently small. It therefore has to be tested whether the optics in the Cornell ERL provides fast enough ion motion to clearing electrodes, so that emittance growth is limited.
Because ions travel to the minima of the electrostatic potential, a clearing electrode has to be located at every such minimum. We therefore first calculate the potential in the beam's center with the approximate equation [17] for a round beam pipe with radius R and a beam with the relativistic factor γ. This equation assumes a uniform transverse density distribution. Fig. 10 shows the relevant beam optics for the Cornell ERL and Fig. 11 shows the resulting potential in the center of the beam.
To find a round beam approximation suitable to our simulations, we first look at the minima of this potential and calculate the β-function that would give the same potential at those locations in the case of a round beam, We can then define a round beam approximation for the entire ERL lattice by using a free drift from one potential maximum to the next with β(s) from Eq. (12). The longitudinal ion velocities near the center of the beam (where most of the ions are located) depend only on this potential gradient. As will be demonstrated by a comparison with particle tracking in Figs approximation. This clearing speed determines the linear ion density also for non-round beams, and our round beam model should therefore be a good approximation. The number of ions created per unit length at the longitudinal position s 0 during a time interval t − t 0 is given by where the number density of the residual gas is ρ n , the ionization cross section is σ and e is the elementary charge. In the equilibrium state this is the same as the rate of ions per unit length originating from s 0 passing through every surface perpendicular to the beam axis, provided that the surface is located downstream from s 0 . During a time interval t − t 0 , the total number of ions per unit interval from s 0 crossing the surface is given by where ρ s0 (s) is the local linear density of ions from s 0 per unit length at the point s. Solving for this quantity, we obtain where the velocity of ions from s 0 at s is given by Integrating over the relevant upstream points from the nearest potential maximum to s yields the total local ion density at s Using this method, and assuming that the residual gas density in the linac sections is reduced by a factor 100 as specified in Tab. I, we find the linear ion densities for the Cornell ERL shown in Fig. 12.
The two sections of the ERL with the largest linear ion densities are located at s ≈ 960m (section A) and at s ≈ 1760m (section B). We therefore studied these sections in detail by performing numerical ion and beam simulations. The resulting equilibrium ion densities can be found in Figs. 13 and 14. The figures demonstrate that our analytical approximation agrees well with the simulations even though the approximation assumes that all ions are in the center of the electron beam. We checked that the agreement becomes perfect in this case.
To reduce the impact of the two sections, additional clearing electrodes can be placed between the maximum and minimum of the potential. We find that by including one extra electrode, the emittance increase can be reduced from 1.04pm to 0.46pm for section A, and from 1.55pm to 0.71pm for section B.
To get a better estimate of the ion effects in the full ERL, we also simulated the ion distribution in a 34m long region around s = 1530m, which corresponds to one of the medium high peaks in Fig. 12. Sending the beam through this section 1000 times results in the phase-space distribution in Fig. 15. We find that after 1000 passes through the ion field, 6.2% of the beam electrons have left the main bunch and have migrated to the four separated islands in phase space. The electronic version of this paper contains an animation of this development of the beam halo.

C. Outlook for non-round systems
While we have argued that our circular beam model should be quite accurate, it would be good to extend our simulations to non-round beams. However, in that case the ion propagation cannot be simulated in the same way using adiabatic invariants. Two independent invariant quantities would be needed to account for oscillations in the x and y directions. The problem is that the action integrals cannot be defined as the area in a phase space plane, because the ion motion is generally not a simple superposition of two independent periodic oscillations in phase-space planes. An example of such phase-space motion can be seen in Fig. 16.

IV. CONCLUSION
We computed the ion distribution produced by an electron beam when ion-clearing electrodes are installed. The transverse ion distributions are shown to strongly peak in the beam's center, producing very nonlinear forces on the electron beam. These can produce strong perturbations of beam properties leading to emittance growth and halo development. These simulations rely on fast algorithms that take advantage of adiabatic invariants and scaling properties of Maxwell's equations and the Lorentz force.
Our results are very relevant for high current Energy Recovery Linacs, where ions are produced relatively quickly, and where clearing gaps in the electron beam cannot easily be used for ion elimination. As an example we used the Cornell Energy Recovery Linac project. For simplicity we only consider the case of a circular electron beam of changing diameter. However, we parameterize this model to approximate the non-round beams of the ERL well. We found suitable places for clearing electrodes and computed the equilibrium ion density and the electron-emittance growth it would produce even after installing extra clearing electrodes at the two most dangerous places. We also showed that the nonlinear-ion forces would lead to a significant halo development.