Transverse-Longitudinal Coupling by Space Charge in Cyclotrons

A method is presented that enables to compute the parameters of matched beams with space charge in cyclotrons with emphasis on the effect of the transverse-longitudinal coupling. Equations describing the transverse-longitudinal coupling and corresponding tune-shifts in first order are derived for the model of an azimuthally symmetric cyclotron. The eigenellipsoid of the beam is calculated and the transfer matrix is transformed into block-diagonal form. The influence of the slope of the phase curve on the transverse-longitudinal coupling is accounted for. The results are generalized and numerical procedures for the case of an AVF cyclotron are presented. The algorithm is applied to the PSI Injector II and Ring cyclotron and the results are compared to TRANSPORT.


INTRODUCTION
There is continuous interest in the understanding of space charge effects in isochronous cyclotrons [1][2][3][4][5][6][7][8][9] . The area is of special relevance for the conceptual designs of cyclotrons as energy-efficient drivers for accelerator driven systems (ADS) 10,11 . Nevertheless there is no selfconsistent first order theory of matched bunched beams with space charge in cyclotrons up to date known to the author. Bertrand and Ricaud showed that it is possible to set up and solve linearized equations of motion for a simple cyclotron model 9 , but they did not determine the parameters of matched beams in cyclotrons and they did not generalize their results in order to find a numerical method to compute matched beams for sector-focused cyclotrons. Information on the matched ellipsoid is of special interest for modern high-performance codes like OPAL, that enable to simulate high intensity beams in the space charge dominated regime 12,13 .
First we describe the usual simplified azimuthally symmetric cyclotron model (see for instance 14 or 15 ), but include the linearized space charge forces. The model neglects possible effects of rf-acceleration and is restricted to the description of a coasting beam. A modified version of Teng's method to decouple longitudinal and horizontal motion will be presented and applied to the simplified analytical model 16,17 . Later we describe how this method can be applied in an iterative numerical procedure to determine the parameters of matched beams with space charge in sector-focused cyclotrons. Finally we present some results as computed for the cyclotrons of the PSI high intensity accelerator facility. a) Electronic mail: christian.baumgarten@psi.ch

I. THE SIMPLIFIED CYCLOTRON MODEL
The simplest relativistic cyclotron model is based on a magnetic field with cylindrical and mid-plane symmetrie. The nominal orbital frequency of an ion beam coasting in an isochronous cyclotron is where q is the charge and m the rest mass of the ions. The axial magnetic field is B = B(r) = B 0 γ. We define a cyclotron length unit a = c ωo using the orbital frequency ω o . The particle velocity is v = ω o r so that β = v/c = r/a. Hence the relativistic factor γ can be written as The axial magnetic field in dependence of the radius r is then given by where ε is a small distortion of the isochronism (ε ≪ 1). The radial field increase compensates the relativistic mass increase and the cyclotron is strictly isochronous for ε(r) = 0. The field index n = r B dB dr is given by so that The phase shift per turn ∆φ ≡ dφ dn is (see for instance 14 : where θ is the azimuthal angle, ω c =θ is the real orbital frequency and ω hf = N h ω o the frequency of the accelerating high frequency system operated at the harmonic number N h . The bending radius r of a particle with charge q in a magnetic field B is given by r = p q B and hence ω c can be written as so that Hence dε dr is proportional to the change of the phase shift per turn with radius. The reference orbit corresponding to the kinetic energy E = m c 2 (γ − 1) = E 0 (γ − 1) is a circle with the radius We will use x = r − r(E) as the horizontal coordinate of a specific orbit relative to the reference orbit, y as the axial deviation from the median plane and z as the longitudinal coordinate.

II. SOME GENERAL REMARKS
From the assumption of midplane symmetrie it follows that the axial coordinate y is not coupled neither to the horizontal nor to the longitudinal coordinate. Therefore it is sufficient to describe the coupling of horizontal and longitudinal motion. Since we assume azimuthal symmetrie, the matched beam ellipse must also be azimuthally symmetric, i.e. constant along the orbital length coordinate s. If the matched beam is constant, then the forces are constant -even if space charge is taken into account.
The vector describing the (horizontal and longitudinal) displacement of an orbit at position s relative to the reference orbit is x = (x, x ′ , z, δ) T where x ′ = dx ds is the horizontal direction angle with the reference orbit and δ = p−p0 p0 is the momentum deviation. We use the formalism of linear Hamiltonian systems as for instance described by R. Talman 18 including the method of A. Wolski to compute the matched beam matrix 19 . The equations of motion (to first order) can be written as where F is a constant 4 × 4-matrix. The general solution is where M(s) is the transfer matrix. The matrix exponential M(s) is Let S be the skew-symmetric matrix with S 2 = −1 so that -since M is sympletic -the following relation holds: If there exists an invertible transformation matrix E and a diagonal matrix λ such that then its is easy to show that the following relation holds: E is the matrix of columnwise eigenvectors and λ is the (diagonal) matrix of the corresponding eigenvalues of F. The (imaginary) eigenvalues of F are the betatron frequencies of the orbital motion. F and M have the same eigenvectors -and the eigenvalues of M are the exponentials of the eigenvalues of F. The eigenellipsoid σ E defined by can then be written as follows 19 : where D is the diagonal matrix with the eigenvalues of σ E S -which are the emittances (apart from factors ±i).

III. THE EQUATIONS OF MOTION (EQOM)
The Hamiltonian is given by where the canonical coordinates are x and l and the momenta are given by x ′ and δ. The EQOM (in first order) are: r is the inverse bending radius and k x = h 2 (1 + n) is the horizontally restoring force. K x and K z represent the strength of the horizontal and longitudinal space charge forces 20 : The eigenvalues of F and M are: Note that b must be positive to give a real-valued frequency ω and hence a longitudinally focused orbit. This is especially important, if we allow for a (small) field error ε(r). In this case the change of the orbital frequency as given by dε dr can play a significant role: If K x < h dε dr , then the longitudinal focusing frequency ω is imaginary and the longitudinal beam size increases exponentially with s. This in fact is a surprising feature of this type of coupling in combination with isochronism: Longitudinal focusing requires a strong enough horizontally defocusing space charge force. On the other hand, if dε dr < 0, i.e. if the radial increase of the magnetic field is below the isochronous field increase, then the longitudinal focusing is strengthened.
The matrix of eigenvectors E of the force matrix F is The inverse matrix E −1 is given by: where The transfer matrix M can be computed by using eq. 16 and one obtains: where C = cos (Ω s) S = sin (Ω s) c = cos (ω s)s = sin (ω s)

IV. THE EIGENELLIPSOID
The matrix D is explicitely given by D = 19 . The matched eigenellipsoid can then be calculated using eq. 18: The beam dimensions are given by the diagonal elements of the matrix representing the eigenellipsoid: where ε 1 = ε x and ε 2 = ε z are identified with the horizontal and longitudinal emittance, respectively. The ax-ial motion can be treated separately and one finds If one considers EQ. 21 together with EQ. 30 and EQ. 31, then it is obvious, that this result does not enable to start a straightforward calculation, since the beam sizes depend on the space charge forces and vice versa in an algebraically complicated way. But with the additional assumption of a spherical beam, it is possible to derive a 4th order equation for the beam size as will be shown in Sec. VII.

V. DECOUPLING LONGITUDINAL AND TRANSVERSE MOTION
If sectors are considered, then the azimuthal symmetrie is broken and hence the force terms in the EQOM and consequently the beam ellipsoid depend on the position s. The vertical beam dynamics can still be treated separately, but is has to be taken into account that the periodic change of the vertical beam size influences the space charge factors K x and K z . The design orbit usually has to be computed numerically as described by Gordon 21 . If this has been done, it is possible to compute the transfer matrix for known starting conditions, since the equations of motion are known and can be integrated. The problem is to find the correct beam dimensions for a given beam current and given emittances such that the beam is matched, i.e. such that EQ. 17 is fulfilled.
Teng and Edwards described a parametrization for coupled motion in two and more dimensions that allows to find the decoupling matrix 16,17 . They called their method symplectic rotation. We will give a brief summary of the method and apply it to the problem of decoupling longitudinal and transverse motion.
Given two 2×2 symplectic matrices A and B that form a block-diagonal (i.e. decoupled) transfer matrix T 0 : where α i , β i and γ i are the familiar twiss-parameters. Then a general symplectic transfer matrix M of coupled motion can be written as where M , N , m and n are 2 × 2 matrices and R is the symplectic rotation matrix. Teng suggested to write R in the form where D is a symplectic 2 × 2 transfer matrix that describes the structure of the coupling. Then one ob- tains 16,17 : The rotation angle φ can be computed using If we apply this method to the transfer matrix computed according to EQ. 16, we obtain: From EQ. 26 we find that A > 0 and B > 0 and B > A since Ω 2 > ω 2 , so that the method fails in the case under study. The method of symplectic rotation is not generally applicable and has to be extended. A workaround solution was found by replacing trigonometic by hyperbolic functions 1 .
with Det(D) = −1, i.e. D is not symplectic -but R is symplectic. Tab. I compares the formulas of the symplectic rotation and of symplectic "Lorentz boost". Applying this method to the case under study then gives: The signs in B suggest a negative β 2 -but this can be compensated by using either a negative emittance or a negative frequency ω, since the transfer matrix is invariant against a change of the sign of ω, while σ E is not. The transformation matrices are explicitely given by: Note: R and R −1 are symplectic. To make the method complete, we give the formula to split 2 × 2 matrices of the form of A or B: The matched beam ellipsoid can be computed for given emittances according to EQ. 18 as soon as the diagonalization of the transfer matrix M is known. The fact that we have to extend the decoupling method of Teng and Edwards raises the question, whether the description of decoupling is now complete or if there are other cases that require further modifications or extensions. The answer to this question requires a proper twodimensional extension of the Courant-Snyder theory and a complete survey of all possible symplectic transformations as presented in an accompanying paper 22 .

VI. THE ITERATION PROCESS AND EXAMPLES
We have shown for the symmetric analytical example that a modified version of the method of Teng and Edwards enables to bring the transfer matrix in blockdiagonal form according to EQ. 32. The diagonalization of the block-diagonal matrices is given by EQ. 43 and hence the matrix of eigenvectors E and the matched beam ellipsoid can be constructed.
In order to take advantage of this method for the case of sectored cyclotrons, the method has to be applied iteratively. The goal is to compute the properties of a matched beam. Assumed that the beam emittances and the current are given as boundary conditions, one proceeds as follows: 1. Provide an initial guess of the beam dimensions σ x (s), σ y (s) and σ z (s).
2. Compute the driving terms of the space charge forces K x (s), K y (s) and K z (s).
3. Compute the one-turn transfer matrix M(s) for all azimuthal angles.

Compute the eigenvectors of M(s) by diagonalization.
5. Use the beam emittances to obtain the eigenellipsoid σ E according to EQ. 18.
6. Take the beam sizes from σ E and go back to step 2, if beam sizes changed significantly compared to the previous iteration.
Optionally one can start the iteration with a reduced beam current and increase it either during the iteration process or use the result of the reduced beam current as an initial guess for higher currents. The (speed of) convergence strongly depends on the assumptions about beam current and emittance. In the following examples, the process converged typically after less than 20 iterations.

A. Example: PSI Injector II
The PSI high intensity proton accelerator facility 23,24 (HIPA) consists of a Cockcroft-Walton pre-accelerator, a four-sector injector cyclotron (Injector II, 72 MeV) and the eight sector Ring cyclotron (590 MeV). In routine operation the beam current is 2.2 mA.
The Injector II cyclotron can be modelled reasonably well by a hard edge approximation of the sector magnets. This has the advantage that the computed matched beam parameters can be compared to TRANSPORT [25][26][27] . We have computed the matched beam parameters for Injector II. The results have then be used as input parameters for TRANSPORT (with space charge). The results are shown in Fig. 1.

B. Example: PSI Ring Cyclotron
The PSI ring cyclotron is a separated sector isochronous cyclotron with 8 sectors. But since the field of the magnets does not fall off as sharp as it does in Injector II, a hard edge approximation does not work equally well. Furthermore the sectors are much more spiralled. For the sake of precision we used the measured field map B(r, θ) to compute the equilibrium orbits (EO) on a radial grid 28 . The radius R eo (θ) of the equilibrium orbit (EO) is then given as a function of the angle and of the starting radius R eo (θ 0 ). Besides the pure geometry of the orbit one also obtains a precise value for ε = 1 − ωo ω for each EO 21 , which then allows to compute dε dr . In order to obtain the force matrix according to EQ. 20, the following additional calculations have been done: 3. Path length element ∆s = √ r ′2 + r 2 ∆θ.
5. h(E, θ), n(E, θ), γ(E) and the space charge forces K x , K y , K z are then used to compose the force matrix F. 9. The average beam sizes are obtained from the eigenellipsoid and used to recompute the force coefficients.
10. The convergence is then checked and either the next iteration starts with step 5) or is stopped. Fig. 2 shows some results of the matched beam with space charge for the PSI ring cyclotron, the phase shift per turn vs. energy in the upper, the transversal tunes in the second and the axial tune in the last graph. The region where the slope of the phase shift is negative has been marked by vertical lines. In this region the longitudinal tune due to space charge forces apparently becomes imaginary.

VII. SPHERICAL SYMMETRY
The first item in the described algorithm is to provide an initial guess of the beam dimensions. In the following we show how to obtain a possible choice of the initial guess for "nearly" spherical beams. If the horizontallongitudinal coupling is strong enough, the beam will usually be approximately circular in this plane 29 . The axial size depends mainly on the vertical emittance and tune.
From EQ. 22 and EQ. 26 one quickly derives: With the assumption of an isochronous cyclotron (k x = h 2 γ 2 ) and spherical symmetrie, i.e.
one finds The beam size σ then yields using EQ. 30: so that one obtains If one defines σ 0 = 2 r ε γ , then EQ. 49 can be written as where A polynom of 4th order has 4 solutions, namely the solutions of EQ. 50 are where the abbreviations used are defined by The minus sign in the square root of EQ. 52 belongs to a pair of complex conjugate solutions and may not be used. Starting from low values of α, Y monotonically increases starting from zero so that we have to take the plus in front of the square root as well in order to obtain positive solutions also for small values of α, so that Analyzing EQ. 50 one finds that so that the function x(α) = (1 + α) 1/3 is a reasonable approximation. The maximal relative deviation of about −0.035 appears at α ≈ 3/2. Therefore we suggest to use the approximations Using the normalized emittanceε = β γ ε, the wavelength λ = c ν rf = 2 π c ωo N h , the cyclotron radius a = c ω0 = r β and the relation ε 0 c 2 = 1 µ0 , the values of K 3 , α and σ 0 can be written as follows: EQ. 56 and 57 describe the matched beam sizes for a spherical beam in a perfectly isochronous cyclotron. The real beam sizes will usually differ from these values since the spherical symmetrie requires that the vertical tune roughly equals the horizontal tune, while in most cyclotrons the vertical tune is below the horizontal tune. Nevertheless the equations derived for the special case of a spherical beam can be used as starting conditions for the described iterative matching procedure.

VIII. SUMMARY
A method has been developed that allows to compute the parameters of a matched beam with space charge in cyclotrons in linear approximation for given beam current and known emittances. As an example, a matched beam in the PSI INJECTOR II cyclotron has been computed. The result has been used as input to TRANSPORT and it has been shown that the beam is matched.
Furthermore it has been shown that the longitudinal focusing as induced by the space charge forces depends on the isochronism of the cyclotron. In case of the PSI Ring cyclotron, the negative slope of the phase shift per turn may reduce or even destroy the longitudinal focusing effect. The results are important for the design of high intensity cyclotrons. If beam emittance and current are chosen accordingly, then the space charge induced longitudinal focusing allows to replace flat-top cavities by accelerating cavities and thus increase the energy gain per turn and the turn separation significantly. This is important to minimize beam losses and activation of components. The flat-top cavities of the PSI Injector II cyclotron are already operating as accelerating cavities and a complete replacement is planned 29,30 .