Nonlinear cooling of an annular beam distribution

In recent years, intense efforts have been devoted to studying how nonlinear effects can be used to shape the transverse beam distribution by means of an adiabatic crossing of nonlinear resonances. By this approach, it is possible to split the beams in the transverse plane, so that the initial single-Gaussian beam is divided into several distinct distributions. This is at the heart of the multiturn extraction process that is successfully in operation at the CERN Proton Synchrotron. Nonlinear effects can also be used to cool a beam by acting on its transverse beam distribution. In this paper, we present and discuss the special case of a beam with an annular distribution, showing how its emittance can be effectively reduced by means of properly devised manipulations based on nonlinear effects.


Introduction
Nonlinear effects introduce new beam dynamics phenomena that might open up the possibility of devising novel beam manipulation techniques. This is the case, for instance, when shaping the transverse beam distribution by means of adiabatic crossing of a stable nonlinear resonance. Such a process is at the heart of the so-called beam splitting that is used for the CERN Multiturn Extraction (MTE) [1,2,3] and has been successfully implemented as a routine part of operation of the CERN Proton Synchrotron since several years [4,5,6].
However, this is not the only nonlinear manipulation that can be devised. Indeed, under the inspiration of [7], it has been found that a controlled redistribution of the invariants can be achieved between the two transverse degrees of freedom [8], provided that an appropriate two-dimensional nonlinear resonance is crossed. This opens novel options in terms of manipulation of the transverse beam emittances.
It is therefore natural to study whether nonlinear effects can be used efficiently to reduce the linear invariants of a transverse beam distribution, thus generating a cooling of the transverse beam emittance. The basis of this approach to beam cooling is the observation that nonlinear effects do not preserve the linear invariant, i.e. the linear action, or the so-called Courant-Snyder invariant. In this sense, they can be used to reduce the value of the linear invariant without violating the symplectic character of the Hamiltonian dynamics. Therefore, the comparison of the value of the linear invariant before and after the action of the nonlinear forces, i.e. when the dynamics is linear and expressed as a rotation around the origin of phase space, is a correct indicator of the reduction of the invariant for each individual particle, and hence of the whole beam distribution and of the corresponding emittance.
In this paper, the initial step towards the development of a nonlinear cooling of a particle distribution is discussed. We present a framework to cool an annular beam distribution, i.e. a distribution with nonzero density in an interval of radii r 1 < r < r 2 , r 1 > 0 in the normalized phase space. It is well-known that annular beam distributions are generated as the result of applying a single transverse kick to a centered beam in the presence of decoherence. Hence, a potential application of annular beam cooling could be the restoration of the initial centered distribution after a transverse kick.
A general discussion of the systems that can be used to devise a cooling method for an annular beam distribution is presented in Section 2, while the considered models are presented in Section 3 together with some results of the theory of adiabatic trapping applied to the models. In the same section, several cooling protocols are defined, and their performance analyzed in detail by means of extensive numerical simulations, whose results are presented and discussed in Section 4. Finally, conclusions are drawn in Section 5, with some mathematical details reported in the Appendices.

General considerations on the model chosen
The general idea underlying the approach developed to achieve cooling of the emittance of an annular beam distribution is based on creating stable islands in phase space. This can be done by slowly modulating the parameters to vary the area of the islands to cause the particles to cross the separatrices. By then moving the resonance islands in phase space their action can be changed and eventually reduced.
To create stable phase-space islands, a resonance needs to be excited. The MTE experience suggests using a Hénonlike map as a model, close to stable low-order resonances, e.g. 1/4, 1/5. If the initial annulus lies outside the chain of islands, then by changing the linear frequency one can act on the area of the central region and of the islands to trap particles in the center. This reduces the action by a quantity equivalent to the area of the islands divided by 2π, according to the separatrix crossing theory.
A simple analysis of the scaling laws of the parameters of the islands, found in [9], suggests that this approach is feasible only for resonances of order n = 4. However, to get the best cooling results one needs two parameters to control the position and the area of the resonance islands. Acting on the sextupolar coefficient is not efficient since this acts as a global-scale parameter [9] and hence changes the dynamic aperture of the map. Therefore, an octupolar kick should be added to the sextupolar one to provide an additional free parameter. The estimates for the area of the islands and the central region can be derived using the results of [9] and [10]. However, the main drawback of this approach is the thick stochastic layer generated by the octupolar kick around the outer part of the separatrix of the four stable islands. This has the effect of inducing the loss of particles, which makes the method unreliable. These observations make the approach based on Hénon-like maps unsuitable for the application under consideration.
Ongoing studies suggest that trapping into islands and transport from within the islands can also be efficiently achieved using AC-modulated magnets [11]. The most straightforward option consists of creating one island using an AC dipole in a 1 : 1 resonance condition, i.e. with the oscillation frequency close to the linear tune of the system. It is worth recalling that AC dipoles have been widely studied in the field of accelerator physics, with essential applications to beam diagnostics (see e.g. [12,13,14,15,16,17,18,19], for an overview of AC dipole studies and applications). A cooling method for annular beams will therefore be devised based on Hamiltonian system modeling of the stable islands used to perform the adiabatic trapping, and subsequent transport, under the influence of an AC dipole.

The Hamiltonian model
Horizontal betatronic motion in the presence of an AC dipole can be described by the Hamiltonian of a generic oscillator with a sextupolar nonlinearity and a dipolar time-dependent excitation [14,12,13], namely where and B 0 ρ stands for the magnetic rigidity of the reference particle, B y is the transverse component of the magnetic field, and is the physical length of the magnetic element. We remark that the choice of the sextupolar nonlinearity is rather arbitrary, as other types of nonlinearity might be used, as long as they generate an amplitude-detuning term.
On the other hand, from the standpoint of applications, the use of a sextupolar nonlinearity is very convenient as it is present in all magnetic lattices of circular accelerators.
Using the action angle coordinates (φ, J) of the unperturbed (ε = 0) system and averaging on the fast Fourier components, the Hamiltonian reads where Ω 2 = g(ω 0 ) k 2 3 and g(ω 0 ) is a function of the linear frequency [9], representing an amplitude tuning term that can be derived using normal forms applied to the Hamiltonian (1). We recall that J(x, p x ) is an adiabatic invariant of the unperturbed system if the frequency ω 0 is slowly modulated.
If we change the coordinates to refer the system to a rotating reference frame with slow angle γ = φ − ωt, taking into account the generating function F = J(φ − ωt) and its time derivative ∂F ∂t = −ωJ, the transformation gives where ψ = ωt.
One can average the fast variable ψ, using yielding the new averaged Hamiltonian which, after a rescaling, can be written in the following form where Equation (7) represents a well-known Hamiltonian [20,21] that can be conveniently written in the form using the Cartesian coordinates X = √ 2J cos γ, Y = √ 2J sin γ. When λ > (3/2)µ 2/3 , a hyperbolic fixed point exists only for Y = 0 and where The phase space portrait of the Hamiltonian (9) is shown in Fig. 1, and it can be divided into three regions: the inner regions G 1 and G 2 (and G 3 = G 1 ∪ G 2 ) and the region outside them.
and J(γ) = 0 for γ = γ 0 with The area of G 1 in polar coordinates is thus given by while the area of G 3 is given by so that where Let us now consider a particle which lies in the outer region with an action J 0 > A 3 /(2π). The area enclosed by its orbit will be A 0 = 2πJ 0 . If we start a slow modulation of the parameters λ = λ(t), µ = µ(t), according to the theory of adiabatic separatrix crossing [20,22], at t = t * , when the condition A 3 = A 0 is met for λ = λ * , µ = µ * , the particle is captured into G 1 or G 2 as a random event. Defining the probability P i of trapping in G i , i = 1, 2 is given by After trapping, the resulting action J is given by A i /(2π), where A i is computed when trapping occurs, namely for λ = λ * and µ = µ * .
Given an initial distribution of particles, all of which have an initial action in the close neighborhood of J 0 , the expectation value of their final action is and we have J ≤ J 0 , since P 1 + P 2 = 1, Hence, the final expected action is smaller than the initial one, i.e. the Courant-Snyder invariant of the particle has been reduced. For a distribution of particles with action J 0 , this results in a cooling of the beam.
Furthermore, when trapping occurs at (λ * , µ * ) we have A 3 = 2πJ 0 and using A 3 = πλ * /2 + K 1 + K 2 = 2πJ 0 , we obtain the expression Substituting K 1 + K 2 into the expressions for A 1 and A 2 , one obtains We note that the values of A 1 and A 2 at the crossing time do not depend on µ * .

Cooling protocols
We envisage three possible protocols to achieve beam cooling, since we can trap particles by varying only λ(t), only µ(t), or both parameters. We will present the three possible processes in this order, referring to them as Protocol A, B and C, respectively.

Variation of λ (Protocol A)
If we keep µ constant, dA i /dt = ∂A i /∂λ · dλ/dt , and the probabilities are thus given by Their expressions have been computed in [20,21] and read where Figure 2 shows J /J 0 as a function of J 0 for different values of µ * . We find that the minimum value of J /J 0 is independent of µ * (the proof is given in Appendix A). A numerical computation of this minimum value gives J /J 0 = 0.3957. Given J 0 , we can always find a value µ that optimizes the cooling, with the final action reduced to ≈ 40% of the initial value.
Strictly speaking, when ε = 0, as in the final state of this protocol, the emittance is not equal to the average value of the adiabatic invariant. The reason for this is that the emittance is computed assuming that the dynamics induces a rotation around the origin, whereas the adiabatic invariant is computed with respect to the fixed point around which the initial conditions actually evolve. In fact, when ε = 0, and especially when particles are trapped both in G 1 and in G 2 , as in the final state of this protocol, they are not rotating around the origin. We observe also that if such a cooled beam were transferred to another accelerator, then its emittance would indeed be equal to the average action of the particle distribution. In this sense, the cooling ratio J /J 0 calculated from Eq. (24) is the lower bound to the actual ratio between the final and initial emittance values.
This situation could be solved or at least mitigated if it were possible to develop a protocol of adiabatic transport that, after the trapping phase, would preserve the actions of the particles while reducing µ to zero. However, one should consider that when trapping is achieved by means of a variation of λ only, the cooling is not particularly efficient, since at best the cooling ratio is ≈ 60%. The methods that we are going to present in the following sections are, in theory, capable of achieving total cooling.

Variation of µ and complete trapping in G 2 (Protocol B)
For the protocol based on the variation of µ, the area derivatives (i = 1, 2) are given by where Thus, we have ξ 1 = −1 and ξ 2 = 2, which means that P 1 = 0 and P 2 = 1. All particles are therefore trapped in G 2 , with an action value Cooling is possible in the interval λ * /4 ≤ J 0 ≤ λ * /2, i.e. 2J 0 ≤ λ * ≤ 4J 0 , which corresponds to the existence of the square roots λ * − 2x 2 c and 6x 2 c − λ * . On the other hand, for λ * > 4J 0 , the initial condition does not belong to the outer region but to the inner region, G 1 . In that case, the separatrix crossing occurs when A 1 = 2πJ 0 and the particle is trapped into G 2 at an action A 2 (λ * , µ * )/(2π). Using the expressions of A 1 and A 2 , we find that the resulting expected final action is which means that cooling is also possible for After being trapped in G 2 , the particle distribution has a smaller action than the initial one, but, as before, the definition of the adiabatic invariant, being µ = 0, is not related to (x 2 + p 2 x )/2. Therefore, a transport process must be designed to reduce µ to zero without losing particles from G 2 . Since the particles are trapped in region G 2 , we need to keep its area constant, i.e. dA 2 /dt = 0, or This can be used to derive a differential equation for µ(λ) Following this equation, as λ is reduced µ increases, and while A 2 remains constant A 1 is reduced to zero, which occurs when µ = (2λ/3) 3/2 . We can then safely reduce both µ and λ to zero, stopping the perturbation: in fact, as µ is kept below (2λ/3) 3/2 no island is present in the phase space.

Coupled variation of λ and µ and complete trapping in G 1 (Protocol C)
One could also devise a protocol in which both λ and µ are modulated. We can express µ as a function of λ, and the expression of the capture probabilities becomes where the prime symbol denotes the derivative w.r.t. λ.
The trapping probability is calculated at the jumping point (λ * , µ * ). Therefore, we can define an implicit function λ * (µ) that resolves the equation A 3 = A 0 (see Fig. 3, left). Then, we optimize the probability by imposing that: (a) all particles are trapped in region G 1 ; (b) the area A 1 is minimized at the trapping point. For the first condition, the equation P 1 = 1, P 2 = 0, gives the following condition on µ Note that the signs of the partial derivatives of A 2 w.r.t. λ and µ ensure that µ < 0.
When P 1 = 1, P 2 = 0 and 2π J = A 1 = λ * − 2J 0 , we can minimize J choosing the minimum λ * for which trapping is possible. This corresponds to A 1 = 0, from which λ * = 2J 0 , and the equation A 3 = 2πJ 0 becomes that can be solved by setting K 1 = πJ 0 and K 2 = 0. From K 1 = πJ 0 we have the equation which is solved when the argument of the arc-sine is 1, so It is straightforward to verify that this implies K 2 = 0. Additionally, this condition induces ∂A 2 /∂µ = 0, and µ diverges. Thus, a perfect cooling, i.e. in which the final value of the action is zero, would require to change µ infinitely fast, which contradicts the adiabatic condition we made to apply the theoretical results.
Although it is not possible to provide an analytical expression for the implicit solution λ * (µ * ) of equation A 3 = 2πJ 0 , we can prove that the graphs shown in Fig. 3 represent the unique solution after having properly scaled the axes. In particular, we find (the details are reported in Appendix A), that the graph of the implicit solution of equation A 3 = 2πJ 0 is independent of J 0 if we rescale λ * → λ * /J 0 and µ * → µ * /J 3/2 0 (see Fig. 3, left). Similar laws hold for the expected cooling J/J 0 , which is a function of the only variable µ * /J 3/2 0 (see Fig. 3, center), and for the required µ , which fulfills the functional relation µ / √ J 0 = f (µ * /J 3/2 0 ) (see Fig. 3, right).

Simulation results
We perform numerical simulations of the dynamics generated by the Hamiltonian of Eq. (1) varying λ and µ according to the protocols previously described. In these simulations, we set ω 0 /(2π) = 0.414, k 3 = 1, and invert the relations of Eq. (8) to obtain the values of ω and ε as a function of λ and µ at each time step. The amplitude-detuning parameter Ω 2 has been evaluated for the unperturbed Hamiltonian at ε = 0 by using the algorithm to evaluate the tune described in [23], to give Ω 2 = −0.3196.
The initial distributions used in the simulations are an infinitely thin annular distribution with initial action J 0 = (x 2 0 + p 2 x,0 )/2, while uniformly distributed according to the angle variable φ 0 = atan(p x,0 /x 0 ), i.e. with the p.d.f.

Protocol A: Cooling by varying λ
This protocol is divided in two phases. The first one is a matching phase, to slowly adapt the initial distribution to the phase space topology, as when µ = 0 the elliptic fixed point is shifted. We will increase µ until the chosen value µ * while keeping λ = 0.
In the first phase, for time t ∈ [0, t 1 ], we set    λ(t) = 0 The actual trapping occurs in the second phase. The parameter λ increases linearly from 0 to a value ∆λ. In order to trap particles at J 0 , one needs ∆λ > λ * , where λ * = λ * (µ * , J 0 ). We then set, for time t ∈ [t 1 , We remark that although the proposed protocol, for the sake of simplicity, envisages two phases of the same duration, it is certainly possible to remove this constraint to adapt the duration of each phase to make it as adiabatic as possible. Figure 4 shows the simulated J /J 0 for different annular distributions ρ J0 as a function of the initial action J 0 using three values of ∆λ (with µ * = 7.5 × 10 −3 ), and compares it with the theoretical estimate given by Eq. (24). We remark that the scale of λ and µ are related with that of J 0 and hence the selected values of µ * do not have any specific meaning, as any change would simply rescale the J 0 axis in Fig. 4.
We observe two effects that are the root of the difference between the theoretical reduction of J /J 0 and the observed behavior. For larger values of ∆λ, the cooling range is increased at the expense of the minimum cooling ratio. Given ∆λ and µ * , for large values of J 0 , λ is never big enough to achieve trapping, since the λ * value that solves A 3 (λ * , µ * ) is larger than ∆λ. Furthermore, increasing ∆λ to trap more particles moves the center of G 2 far from the origin of the phase space (all fixed points of Eq. (9), from the solution of the resulting cubic equation, are O( √ λ) for large values of λ), thus decreasing the effective cooling ratio.

Protocol B: Cooling by varying µ
This protocol consists of three phases: the first phase is used to perform particle trapping, with the second and the third needed to transport the particles back to the center of the phase space by progressively reducing the strength of the AC dipole. In the first phase, for times t ∈ [0, t 1 ], we have the following.
In the second phase, the differential equation (32) is solved. For t ∈ [t 1 , t 2 ], we set λ(t) = λ * −λ (t − t 1 ) and obtain µ(t) by numerically integrating the Cauchy problem where dµ/dλ is given by Eq. (32). The second phase is stopped at time t 2 once the condition µ(t 2 ) = µ 2 = (2λ(t 2 )/3) 3/2 is met. The third phase follows for times t ∈ [t 2 , t 2 + t 1 ], with The plots of the time evolution of λ and µ are shown in Fig. 5. We note that the theory presented earlier accurately describes the simulated cooling ratio unless it is in the vicinity of λ * = 4J 0 = 0.2, where the theory predicts total cooling, while in simulation, J /J 0 ≈ 10%. This is due to the angular dependence we averaged upon in our analysis, as can be inferred from Fig. 7. This figure shows the distributions at the end of each of the three phases of Protocol B for the same initial annular distribution for three values of λ * . We observe that at the end of each phase the action of the particles, which were all the same at the beginning, are spread according to their initial phase. For example red particles, which correspond to the initial phase π, result in the innermost position when λ * = 0.15 and in the outermost position when λ = 0.25. This behavior reverses for cyan particles, which have φ 0 = 0. This means that particles with different initial angles are trapped at slightly different values of J. Some particles are trapped earlier or later than expected, with a larger or smaller value of J than that given by theory. In the graphs, it is also visible that the inner and outer particles are reversed, depending on whether λ * < 4J 0 or λ * > 4J 0 . When λ ≈ 4J 0 , however, all particles are trapped at a higher value than expected no matter when they cross the separatrix, thus increasing J . In our simulations, we were able to reach J /J 0 = 0.078, for a cooling efficiency of 92%.
In Fig. 8 we show the dependence of the cooling ratio on the value of the initial action J 0 for three values of λ * . The range in which J /J 0 < 1 represents the possible interval of actions of a thick annular distribution that could be cooled using Protocol B. Note that according to the theoretical predictions cooling is possible in the range λ * /6 ≤ J 0 ≤ λ * /2 and the optimal cooling ratio is found at J 0 = λ * /4.
An animation of the trapping process for a thick annular distribution is available as Supplemental Material 2

Protocol C: Cooling by varying λ and µ
This protocol requires two phases: the first to adapt the phase space; the second for trapping and transport. Our goal, besides trapping the particles inside G 1 , is to ensure that both at the beginning and at the end of the process the adiabatic invariant is as close as possible to the linear action variable J = (x 2 + p 2 x )/2, which is true if the AC dipole is switched off, i.e. when µ = 0. Thus, in the first phase, µ is gradually increased, while keeping λ = 0 (i.e. ω = ω 0 ), until it reaches the value needed to initiate the trapping process. In the second phase, the derivative of µ(λ) is kept at a constant value µ while increasing λ, and taking advantage of the fact that as µ < 0, we can slowly reduce µ until it reaches zero to recover the equivalence between the adiabatic invariant and J.
In the first phase, for times t ∈ [0, t 1 ], we set where µ max = µ * + λ * |µ |. This ensures that during the second phase when λ = λ * , µ is exactly µ * and its derivative µ has the appropriate value. The values of µ * and λ * are obtained by choosing a solution of the implicit equation  In the second phase, where t ∈ [t 1 , When the process ends and µ = 0 is reached, G 2 disappears as the perturbation provided by the AC dipole has been switched off, and the particles trapped in G 1 have been transported to the center of the phase space. The values of λ and µ during the whole procedure are plotted in Fig. 9.
We remark that although the proposed protocol envisages two phases of the same duration, it is possible to remove this constraint to adapt the duration of each phase to make them as adiabatic as possible.
In Fig. 10 we show the simulated cooling ratio J /J 0 for an initial annular distribution ρ 0.05 (φ, J), as a function of µ * , and a comparison with the theoretically expected value J = (λ * (µ * ) − 2J 0 )/(2π). It can be seen that the agreement between theory and simulation is remarkable up to a certain breakdown value of µ * . This breakdown is due to the angular dynamics that has been neglected in the averaging process of the theory. In Fig. 11 we show the initial distribution, the situation at the end of the first phase and the final distribution of particles for two different values of µ * , using the hue to represent the initial angle φ 0 . For both values of µ * , we observe that the distribution after the first phase is no longer infinitely thin, and that the action of each particle depends on the initial angle. As a result each particle crosses the separatrix at a different time at the end of the second phase resulting in different values of the final action. For µ * smaller than the breakdown threshold, all particles are still trapped in G 1 , and this angular dependence is averaged out. On the other hand, for higher values of µ * , particles that at the end of the first phase are in the outer part of the distribution can also be trapped in G 2 at high amplitude, thus dramatically increasing the value of the final action. We again stress that we cannot expect to reach 100% cooling as |µ | and µ max would need to reach unlimited values. The best cooling that we could achieve in our numerical simulations is 92%, at J /J 0 = 0.08.
To study the applicability of the cooling protocol to a more realistic particle distribution, we have looked at an ensemble of infinitely thin annular distributions ρ J0 covering a certain interval in J 0 . The values of µ and µ max have been chosen to optimize the trapping for a particular value of J 0 ,Ĵ 0 = 0.05. The results are shown in Fig. 12. It is clearly visible that for different values of µ * , which translates into different cooling targets for particles atĴ 0 , a significant range of action values is actually cooled. The width of this cooling well, i.e. the range of J 0 where J /J 0 < 1, is the thickness of the annular distribution that the protocol can handle successfully.  Figure 9: Evolution of λ(t) and µ(t) during the two phases of Protocol C. The two values λ max and µ max have expressions in function of the computed λ * , µ * , and µ , i.e. µ max = µ * + λ|µ |, λ max = λ * + µ * /|µ |.
We note that, contrary to theoretical expectations, the minimum value of J /J 0 does not occur atĴ 0 , although this difference tends to decrease as µ * increases. This is due once more to the angular dynamics. Using the same parameters as the plots shown in Fig. 12, two final distributions are shown in Fig. 13 using the hue of the color to identify the initial phase. The right plot shows the case where the initial distribution is ρ 0.05 , i.e. the initial conditions are selected atĴ 0 , while the left plot shows the case where the initial distribution is ρ 0.045 , where the initial actions have a value J 0 <Ĵ 0 , but close to the minimum. In the left plot, a gap in the final distribution is clearly visible. This can be explained by the fact that in this case some particles are trapped earlier (the red dots in the plots) due to the spread of the action after the first phase. These can end up either in G 2 or in G 1 , according to the probability law, but when their areas are smaller. The average final action is therefore reduced more by this effect than by the increase induced by the particles in G 2 .
An animation of the trapping process for a thick annular distribution is available as Supplemental Material 3

Conclusions
In this paper, beam manipulations based on nonlinear beam dynamics have been devised with the goal of achieving cooling for annular transverse beam distributions. Such a beam distribution can be generated after a beam is kicked in the transverse direction. The possibility of achieving cooling by means of crossing stable resonances generated by static magnetic elements has been ruled out, however the use of an AC dipole for such manipulations has proven to be very successful. A Hamiltonian model describing the transverse dynamics in the presence of an AC dipole has been studied using concepts from the adiabatic theory for Hamiltonian systems. This has allowed the design of three cooling protocols, two of which proved to be extremely effective with a simulated best performance of ≈ 90% cooling. In physical terms this observed cooling is achieved by controlling the strength and frequency of the AC dipole according to the specifications of the proposed protocols.
Detailed numerical simulations carried out on the considered Hamiltonian systems have revealed a rich phenomenology that could be explained in detail by using adiabatic theory for Hamiltonian systems. Although an infinitely thin annular distribution was initially used, the two best protocols have been shown to have a significant cooling range. It therefore seems possible to be able to use them to cool a transverse annular beam distribution of finite thickness. Numerical studies on more realistic accelerator models will be considered in the future in view of experimental tests on a real machine.
Such annular beam distributions are also representative of the beam halo, which opens up the study of future applications to halo manipulation that could result in experimental tests at the LHC.  (1) has been used, with k 3 = 1, ω 0 /(2π) = 0.414, Ω 2 = −0.3196, t 1 = 1 × 10 5