Dynamical Stability of Slip-stacking Particles

We study the stability of particles in slip-stacking configuration, used to nearly double proton beam intensity at Fermilab. We introduce universal area factors to calculate the available phase space area for any set of beam parameters without individual simulation. We find perturbative solutions for stable particle trajectories. We establish Booster beam quality requirements to achieve 97\% slip-stacking efficiency. We show that slip-stacking dynamics directly correspond to the driven pendulum and to the system of two standing-wave traps moving with respect to each other.


Introduction
Slip-stacking is integral to high-intensity operation at Fermilab and will likely play a central role in upgrades to the accelerator complex [1][2] [3]. Particle loss in the slip-stacking process is a limiting factor on ultimate performance [1] [4]. Single-particle dynamics associated with slip-stacking contribute directly to the particle losses. This paper analyzes these dynamics at depth, both analytically and numerically. Our numerical results completely characterize the stable phase-space boundary. We use these results to recommend an upgrade to the Fermilab Booster that would substantially reduce slip-stacking losses.
Our analytical results provide insight into slip-stacking by presenting a perturbative general solution and new parameteric resonances. These results should also be of interest to the greater field of dynamical mathematics because, as we demonstrate, the dynamics of slip-stacking are isomorphic to the well-studied dynamics of the driven pendulum. The analysis in this paper is also intended to facilitate application of slip-stacking to other accelerators and non-accelerator systems with analogous dynamics.

Background
Slip-stacking is a particle accelerator configuration that permits two high-energy particle beams of different momenta to use the same transverse space in a cyclic accelerator. The two beams are longitudinally focused by two sets of rf cavities with a small frequency difference between them. Each frequency is tuned to the momentum of one of the beams.
The two azimuthal beam distributions are manipulated as a consequence of their difference in rf frequency. The two beams injected on separated portions of azimuth with a small frequency difference will overlap gradually, allowing injection [4]. When the cyclic accelerator is filled, the azimuthal distribution of the beams will coincide at a certain tune and can then be accelerated simultaneously. The accelerating rf cavities operate at the average frequency, capturing both beams as one. The potential beam intensity of a synchrotron is doubled through the application of this technique.
A preliminary study explored the beam dynamics in a 2-rf system [5]. The slipping of bunched beams was first demonstrated at the CERN SPS [6] but the emittance growth led to unacceptable particle losses. Fermilab has subsequently implemented slip-stacking operationally since 2004 [7][8] [4]. Initially, the higher beam intensity was used to increase antiproton production for proton-antiproton collider experiments [9]. Subsequently, slip-stacking was applied to neutrino production for Neutrinos at Main Injector (NuMI) experiments [10] [11][12] [13].
Beam-loading effects can impact the effectveness of slip-stacking and were addressed in the Main Injector by the development of a beam-loading compensation system with -14dB feedback and -20dB feedforward [14][15] [16]. The beam-loading effects on slip-stacking in the Recycler will be an order of magnitude weaker than in the Main Injector and can be compensated if necessary. The typical beam-loading voltage is ∼2kV [15] compared to a typical rf cavity voltage of 90kV [16]. In the Main Injector the R sh /Q of the rf cavities is 100Ω [15], while in the Recycler the R sh /Q is 13Ω [17]. This paper focuses on the constraints on the stable phase-space area from the single-particle dynamics of the two-rf system; direct space-charge effects are an order of magnitude weaker.
For the particular case of slip-stacking at Fermilab, the difference between the two rf frequencies must be equal to the product of the harmonic number of the Booster rf and the cycle rate of the Fermilab Booster. The cycle-rate of the Booster is 15-Hz and therefore ∆f = h B f B = 1260 Hz. A possible upgrade to a 20-Hz Booster is also analyzed, for which ∆f = 1680 Hz. A 20-Hz Booster would enable slip-stacking buckets with substantially greater phase-space area. However, the slip-stacking cyclic accelerator (either the Main Injector or the Recycler) must be able to simultaneously accommodate beams in a range of momentum corresponding to their frequency difference.

Single-rf Longitudinal Stability
The motion of a single particle under the influence of a single stationary rf cavity is described in terms of its phase space coordinates φ and δ. The phase φ is the phase of the particle relative to the resonating electromagnetic field in the rf cavity. δ is the fractional deviation from the reference momentum; δ = 0 corresponds to a particle whose revolution frequency f rev is a subharmonic of the frequency of the rf cavity f rf = hf rev . The phase-slip factor η is used to describe how the revolution period T changes with δ and is given by ηδ = ∆T /T (see [18]).
The equations of motion associated with the trajectory of a single particle under the influence of a single stationary rf cavity [18] are given by: (1) V is the effective voltage of the rf cavity, e is the charge of the particle, β = v/c is the velocity fraction of the speed of light, E is the total energy of the particle. Let V δ be equal to eV β 2 E , the maximum change in δ during a single revolution.
For small φ, Eq. 1 has a stable solution known as a synchrotron oscillation where the synchrotron frequency ω s is given by ω s = 2πf rev V δ hη 2π . The amplitude ρ and the initial phase ψ are set by initial conditions. More generally, stable oscillatory motion is bound in phase and momentum by the separatrix The equations of motion given in Eq. 1 are isomorphic to that of a simple pendulum. The stable region of phase space within the separatrix is referred to as the rf bucket and the region outside of the separatrix is referred to as slipping with respect to the bucket.
In contrast, the dynamics of slip-stacking are explicitly time-dependent and there is no simple separatrix delineating the bucket boundary. The lack of a clearly defined bucket confounds beam operation as the incoming particles cannot be conventionally inserted into the bucket. We broaden the term bucket to include cases without a separatrix: a particle trajectory is in a particular rf bucket if the particle phase with respect to the rf cavity is bounded and averages to zero.

Slip-stacking and the Driven Pendulum
The equations of motion for a single particle under the influence of two rf cavities with identical voltage and different frequencies are: ω φ is what we refer to as the phase-slipping frequency, the angular frequency separation between the two rf cavities ω φ = 2π∆f . φ A is the average of the phases and φ D is the difference between the phases for the two rf cavities at t = 0. Without loss of generality, we eliminate φ D with t → t − φ D /ω φ . Applying a substitution to Eq. 4 yields the corresponding second-order equation of motion: The corresponding Hamiltonian is then given by: This Hamiltonian leads to nonlinear, resonant, and chaotic phase-space trajectories. We find that this Hamiltonian is isomorphic to that of a pendulum under a sinusoidal driving force (in the absence of gravity). The driven pendulum is a type of nonlinear Mathieu equation that is a subject of ongoing research in computational mathematics [19] [20]. A canonical form for the driven pendulum may be parameterized as [19]: Eq. 7 is obtained from Eq. 5 under the substitution The parameter a corresponds to the force of gravity and a = 0 in this analogue. The accelerator literature [5][6] [7] has identified the importance of the slip-stacking parameter as the criterion for effective slip-stacking. The parameter b in Eq. 8 is a function of α s , indicating that all nontrivial dynamics of slip-stacking depend only on α s . For example, if one slip-stacking configuration has phaseslipping frequency ω φ and another configuration with the same α s has phase-slipping frequency ω φ then the second phase space diagram is isomorphic to the first where the δ axis must be scaled by ω φ /ω φ .

General Perturbative Solution
We analyze the dynamics within a slip-stacking bucket by shifting the origin from one synchronized to the average frequency and phase of the two rf cavities to one synchronized to the rf cavity with higher frequency: Substituting one equation into the other we have: Eq. 11 can be expanded into powers of φ to consider the small φ perturbation: We use the Poincare-Lindstedt method (see Ch. 2 of [21]) to find the perturbative solution to Eq. 12 as a linear combination of oscillatory terms. We substitute the case φ = 0 into Eq. 12 and we solve to generate φ = −α 2 s sin(ω φ t). Next we use φ = −α 2 s sin(ω φ t) to generate A n sin(nω φ t) terms. The coefficients A n are of the order α −2n s and do not depend on the initial coordinates of the particle. These terms form the particular solution: There is no stable equilibrium point inside of the bucket. The particular solution is analogous to a moving bucket center; we term the trajectory of a particle at the moving bucket center to be quasi-synchronous because the frequency spectrum will depend only on harmonics of ω φ . We continue the perturbation to the general case by using φ = φ p +ρ sin(ω s t + ψ), the sum of the particular solution (Eq. 13) and the small-oscillation single-rf solution (Eq. 2). Using φ = φ p + ρ sin(ω s t + ψ) generates terms of the form B m,n sin[m(1 + σ)ω s t + nω φ t + mψ]. The shift in the synchrotron oscillation frequency σ is a necessary contribution to the coefficient of sin[(1 + σ)ω s t + ψ] in Eq. 12 to counterbalance the contribution made by the cross-multiplication of higher order B m,n sin[m(1 + σ)ω s t + nω φ t + mψ] terms. For any integer m > 0 and any integer n, the coefficients B m,n are of the order ρ m α −2|n| s ; except when m is even and n = 0, in which case the coefficients B m,0 are of the order ρ m α −2 s . Writing out the full perturbative solution, we have: The trajectory of a particle in a slip-stacking rf bucket is referred to as a rotating solution in the driven-pendulum literature. The particular solution was previously obtained by Zhang and Ma [22]. An alternate perturbative approach for the general solution in implicit form is given in [23]. We are first to find a general and explicit solution.
The perturbative solution for the small oscillations around the moving bucket center can be expressed in coefficients up to order α −4 s and ρα −2 s . The derivation shown in the appendix leads to the equations of motion: The parameters ρ and ψ are determined by initial conditions, shown explicitly in the appendix. We are the first to discover and calculate σ, the synchrotron frequency shift from a slip-stacking perturbation. Generally, particles within a slip-stacking bucket will undergo synchrotron oscillations at a higher frequency than the corresponding single-rf bucket. Substituting the perturbative terms from Eq. 14 into Eq. 11 indicates that a new parametric resonance will occur wherever mω s (1 + σ) = nω φ . For example, the ρ m α −2(n−1) s sin[m(1 + σ)ωt − (n − 1)ω φ t + mψ] term will be multiplied by cos(ω φ t) in Eq. 11 and lead to a growth term proportional to ρ m α −2(n−1) s sin(mψ). The case where mω s = ω φ was previously investigated by Mills [5]. An analytical description of the stable phasespace boundary would require a complete determination of the cases in which parametric resonances lead to particle loss; this may be the subject of future work.

Stability Maps & Area Factors
The size and shape of slip-stacking buckets determine which portion of an injected beam distribution is lost. In application, lost particles migrate to an incorrect azimuthal location and consequently collide into the beampipe during injection, extraction, or acceleration [4]. We map the stability of initial particle positions by integrating the equations of motion for each position. The integration is iterated for a sufficiently large number of revolutions (at least 30 synchrotron periods). The stability of the particle is tested after every phase-slipping period. A particle is considered lost if its phase with respect to each of the first rf cavity, the second rf cavity, and the average of the two rf cavities, is larger than a certain cut-off (we used 3π/2). The remaining particles therefore belong to one of four stable regions shown in Fig. 1: one for the higher frequency, one for the lower frequency, one for the average frequency and average phase, and one with average frequency but π offset from the average phase. These two stable regions at the average frequency are the original examples of dynamic stabilization [24]. Fig. 2 shows the stability of initial coordinates in the higher bucket for α s = 3.6 and α s = 4.1, in which the effects from slip-stacking resonances are evident. The supplemental material [25] shows the stability maps of the higher slip-stacking bucket for values of α s from 2 to 8 in descending 0.1 increments.
We find some trajectories are "metastable" because they lead to particle loss only after thousands of revolutions. The stable phase-space area as a function of time is shown in Fig. 3 for several values of α s . The color corresponds to the number of synchrotron periods a particle with the corresponding initial coordinates survives before it is lost. The two large stable regions correspond to the higher and lower rf buckets where beam is injected and maintained. The two stable regions along the δA = 0 axis are created by the interaction between the two rf cavities.
The bucket area is computed as the product of the total number of ultimately surviving points and the cell area. We define the slip-stacking area factor F (α s ) = A s /A 0 as the ratio of the slip-stacking bucket area to that of a single-rf bucket with the same rf voltage and frequency. The area factor follows the notation of Lee (Ch. 3.II of [18]) for accelerating beams, in which the ratio of running bucket area to stationary bucket area is used. Particles in the bucket are described by Eq. 14 with finite coefficients, therefore the bucket area is conserved. Consequently F (α s ) does not depend on the initial rf phase difference used to generate the stability map. We write the phase space area (φ · δ units) using F (α s ): The color corresponds to the number of synchrotron periods a particle with the corresponding initial coordinates survives before it is lost. On the left, αs = 3.6 and four resonance islands can be seen due to the ωs(1 + σ) = 4ω φ resonance. On the right, αs = 4.1 and five resonance islands can be seen due to the ωs(1 + σ) = 5ω φ resonance. Fig. 4(a) plots the numerically derived slip-stacking area factor F (α s ). Using Fig. 4(a) with Eq. 21 provides the first method for calculating the slip-stacking stable phase-space area without requiring each case to be simulated individually. F (α s ) increases rapidly above α s ≈ 3 and asymptotically approaches 1. F (α s ) has several local minimum where resonances are crossed; this loss of area occurs when large amplitude trajectories have a parametric resonance and therefore does not occur at precise integer values of α s .
In application, slip-stacking is tuned to maximize stable phase-space area while holding ω φ constant. The value of ω φ is generally constrained by gross features of the accelerators, for example the harmonic number and cycle time. The slip-stacking parameter α s is tuned through changing ω s which is proportional to the square root of the applied rf voltage. Furthermore ω s changes the bucket area by both the slip-stacking area factor F (α s ) and the single-rf bucket area, so there is an optimal voltage in which phase space area is maximized. We rewrite Eq. 21 to separate the parameters that are The rapid losses at the beginning corresponds to regions of phase space in which particles rapidly slip by both slip-stacking buckets. In the next phase, the metastable particle loss occurs asymptotically. As αs increases the distance between the rf buckets becomes greater, the buckets become more independent, and the slipstacking bucket area approaches the single-rf bucket area. (b) The modified slip-stacking area factor as a function of αs. The modified slip-stacking area factor is maximized near αs = 6.2.
held constant from those dependent on α s : This modified area factor Z(α s ) is graphed in Fig. 4(b). Z(α s ) is maximal near α s = 6.2 and when considering other optimization criteria 5.5 to 7 is is a practical tuning range for α s .

Injection Efficiency, Emittance, and Aspect Ratio
The stability maps can also be used to analyze injection scenarios, by weighting the (scaled) stability maps according to a distribution that represents the number of incoming particles injected into that region of phasespace. We used this technique to identify the greatest longitudinal emittance an incoming Gaussian-distributed beam could have and still achieve 97% injection efficiency at its optimal value of α s . The longitudinal beam emittance is given in Eq. 23 below: The current accelerator upgrade proposal, Proton Improvement Plan II (PIP-II) [1], defines a minimum 97% slip-stacking efficiency required to maintain current loss levels while increasing intensity. Fig. 5 shows the 97% longitudinal emittance as a function of aspect ratio and demonstrates the consequences of a mismatched injection into a slip-stacking bucket. Fig. 6 shows the optimal value of α s as a function of aspect ratio. The optimal value of α s determines the optimal rf cavity voltage, shown in Fig. 7. These results were obtaining using parameter values specific to slip-stacking in the Fermilab Recycler (see Table. I).  Table I. A nominal value for the Booster emittance is 0.12 eV·s [26]. The Fermilab Booster uses bunch rotation via quadrupole excitation [27] [28], with parameters that are actively tuned to minimize losses. With bunch rotation,  Table I.  Table I. the aspect ratio of at least 2.6 MeV/ns is achievable at extraction from the Booster [26]. At Recycler rf cavity voltage V 0 = 100kV, the slip-stacking parameter for the Recycler is α s (V 0 ) ≈ 4.39 for a 15-Hz Booster cycle-rate and α s (V 0 ) ≈ 5.86 for a 20-Hz Booster cycle-rate. For other voltages, the Recycler slip-stacking parameter is given by α s (V ) = α s (V 0 ) V /V 0 .
We examine the 97% efficiency benchmark not only for a 15-Hz Booster cycle-rate but also for a proposed 20-Hz  Booster cycle-rate. A 20-Hz Booster cycle-rate would necessarily increase the phase-slipping frequency (rf frequency separation) by a factor of 4/3 and therefore the bucket height would increase by a factor of 4/3 (and have the same α s ). Slip-stacking losses could be reduced by a factor 4-10 for the same emittance [29]. Alternatively, the 97% efficiency benchmark is achievable at an emittance up to a factor of ∼1.7 greater. Other implications of a 20-Hz Booster cycle-rate are discussed in a Fermilab technical memo [29]. A 20-Hz Booster cycle-rate is clearly superior for high-intensity operation.
The scaling symmetry used to analyze the 20-Hz Booster cycle-rate can generalized. An optimization at phase-slipping frequency ω φ and aspect ratio r is equivalent to an optimization at phase-slipping frequency ω φ and aspect ratio (ω φ /ω φ )r. The same optimal slipstacking parameter would be obtained at a higher synchrotron frequency (ω φ /ω φ )ω s , increasing the rf voltage at (ω φ /ω φ )r to (ω φ /ω φ ) 2 V .

Application to Other Physical Systems
In general, the dynamics discussed in this paper apply to any system governed by (nearly) identical sinusoidal potentials moving with respect to each other (or equivalently, a sinusoidal potential oscillating in amplitude). This is relevant to standing wave traps which are used in optical and acoustic physics and are instances of a controllable sinusoidal potential. Optical lattices are a type of standing wave trap used in ultracold atomic physics. Two optical lattices moving with respect to each other could occupy the same transverse space yet focus two groups of atoms with independent momenta. We know of no such experiment that utilizes this optical slip-stacking laser configuration, but we believe it may be relevant to at least two optical applications: trap accumulation [30] and atomic collisions [31]. If M is the mass of the atom, V H is the potential barrier height, and v is the relative velocity between the two standing wave traps, then the optical slip-stacking parameter is given by: The relative velocity v can be calculated from the frequency difference in the standing waves v = ∆f [32]. Our approach suggests a value of α s equal to at least 5 for efficient stacking. Acoustic standing waves can be used to trap small spheres or droplets in a sinusoidal potential originally described by Gor'kov [33]. This technique has grown in sophistication and application [34] [35]. Eq. 24 would also determine the stability of objects (with mass M ) in a possible acoustic slip-stacking configuration.

Conclusion
In summary, we have provided a framework for addressing both the trajectory and stability of particles in a slip-stacking potential. We introduce the slip-stacking area factor F (α s ) and the modified area factor Z(α s ) as tools to calculate the stable slip-stacking bucket area for any combination of accelerator parameters. We introduce the quasi-synchronous particle trajectory and provide a perturbative solution near it. We describe a series of new parametric resonances in slip-stacking. We provide a general method for analyzing slip-stacking injection scenarios and describe the implications for the operation and upgrades of the Fermilab Booster. We identify for the first time how the dynamics of slip-stacking correspond to the driven pendulum and moving standing wave traps.

Acknowledgments
This work is supported in part by grants from the US Department of Energy under contract DE-FG02-12ER41800 and the National Science Foundation NSF PHY-1205431. Special thanks to SY Lee for providing a crucial mentoring role immediately prior to the beginning of this research.

Appendix A: Derivation of Perturbative Solution
In this paper we show that the second-order equation of motion for a single particle in a slip-stacking bucket is given by Eq. 11 and the perturbative solution is given by Eq. 14. Recall also, that the coefficients A n are of order α −2n s and B m,n are of order ρ m α −2|n| s except the B 2k,0 coefficients which are of the order ρ 2k α −2 s . The coefficient B 1,0 is defined to be equal to ρ. The parameters ρ and ψ are set by initial conditions. For clarity, we adopt a short-hand notation for the oscillatory terms as follows: In this appendix, we explicitly obtain the perturbative solution up to order α −4 s and ρα −2 s (or equivalently, all coefficients up to order α −5 s with ρ ∼ α −3 s ). Therefore, we use the coefficients A 1 , A 2 , ρ, B 1,1 and B 1,−1 ; all other coefficients are neglected at this precision. We start by assuming a solution form and will demonstrate it to be self-consistent: φ = A 1 s 0,1 +A 2 s 0,2 +ρs 1,0 +B 1,1 s 1,1 +B 1,−1 s 1,−1 . (A2) It will be sufficient to substitute this expression into the form of Eq. 11 expanded up to second order in φ: Here we have divided both sides by ω 2 φ to make the order of the perturbation terms more explicit.

(A5)
where the crossed-out terms are higher order than the precision of this analysis. We rewrite the RHS without these immediately negligible terms, using the trigonometric product-to-sum rules, and grouping by the oscillatory term: We then equate Eq. A4 with Eq. A6 for all time. Each oscillatory term corresponds to its own equation: Solving Eq. A8 for A 2 we obtain: Solving Eq. A7 for the linear A 1 terms we obtain: Since A 1 of the order α −2 s (as expected) then the α −2 s A 2 term and the α −2 s A 2 1 are of order α −6 s and are neglected.
Rewriting Eq. A7 to reflect this, we have: All expressions on the RHS side of Eq. A12 are of order α −6 s , therefore we are self-consistent to exclude the A 3 term.
We solve Eq. A10 and Eq. A11 for B 1,±1 and obtain: where σ makes a negligible contribution to B 1,±1 . All expressions on the RHS of Eq. A13 and Eq. A14 are of order α −4 s ρ, therefore we are self-consistent to exclude the B 1,2 and B 1,−2 terms.
For Eq. A9 we call the σ 2 terms negligible, subtract the bare ρ term from each side and solve for σ to obtain the shift in synchrotron frequency: To calculate Eq. A19 first we must calculate: (A20) We substitute Eq. A20 into Eq. A19 to obtain: To summarize, we write these coefficients together: Derivation of ρ and ψ from Initial Conditions At this the synchrotron amplitude ρ and initial synchrotron phase ψ are still undetermined, but we can express them in terms of the initial coordinates. At time t = 0 (which is fixed to a time in which the relative phase between the rf cavities is zero), let δ = δ 0 and φ = φ 0 .
We have shown that φ takes the form give in Eq. A2. We evaluate this expression for φ at t = 0, leave the short-hand notation (Eq. A1), and find φ 0 in terms of ρ and ψ: We can calculate δ from our solution for φ by taking the derivative: We evaluate this expression for δ at t = 0, leave the short-hand notation (Eq. A1), and find δ 0 in terms of ρ and ψ: Next we solve Eq. A26 and Eq. A29 for Φ 0 = ρ sin(ψ) and ∆ 0 = ρ cos(ψ): φ 0 and δ 0 have been translated by the initial position of the bucket center and rescaled to obtain the expressions for Φ 0 and ∆ 0 . Using Φ 0 = ρ sin(ψ) and ∆ 0 = ρ cos(ψ), the solution for ρ and ψ be written as: Eq. A30 and Eq. A31 can be further simplified by writing the B 1,1 and B 1,−1 terms explicitly in terms of α s . We calculate: We apply Eq. A20 and Eq. A34 to Eq. A30 and Eq. A31 to obtain: Φ 0 and ∆ 0 are fully expanded as follows: Eq. A38 and Eq. A39 can then be substituted into Eq. A32 and Eq. A33 to obtain ρ and ψ respectively. Then ρ can be substituted into Eq. A24 to obtain B 1,1 and B 1,−1 .
[2] C. Mariani (LBNE/DUSEL Collaborations), in Proceedings of the Neutrino Oscillation Workshop, Otranto, This supplemental material shows the stability maps of the higher slip-stacking bucket for values of α s from 2 to 8 in descending 0.1 increments. Recall that stability map is a visualization of a numerical integration of the following (1) Where φ is the phase of the particle relative to the rf cavity and δ is the fractional deviation from a reference momentum. Each pixel in a stability map represents a different set of initial coordinates and the color corresponds to the number of synchrotron periods a particle survives before it is lost. The stability maps are scaled to accurately reflect the bucket shapes obtained in the Fermilab Recycler with a 15-Hz Booster cycle-rate. The corresponding Recycler rf voltage is listed for each plot. Some stability maps are particularly noteworthy. Recall that α s = 6.2 maximizes the phase-space area and α s = 5.5 maximizes the Booster injection efficiency. In these figures one can clearly see the 7th-order resonance (α s = 5.0 to 5.1), the 6th-order resonance (α s = 4.5 to 4.8), the 5th-order resonance (α s = 4.0 to 4.3), the 4th-order resonance (α s = 3.5 to 3.7), and the 3.5-order resonance (α s = 3.2), the 3rd-order resonance (α s = 2.9 to 3.0), the 2.5-order resonance (α s = 2.5), and the 2nd-order resonance (α s = 2.2 to 2.3). For α s = 2.1 and lower, there is no stable phase-space area.