Analytic fluid theory of beam spiraling in high-intensity cyclotrons

Using a two-dimensional fluid description, we investigate the nonlinear radial-longitudinal dynamics of intense beams in storage rings and cyclotrons. With a multiscale analysis separating the time scale associated with the betatron motion and the slower time scale associated with space-charge effects, we show that the longitudinal-radial vortex motion can be understood in the frame moving with the charged beam as the nonlinear advection of the beam by the $\mathbf{E}\times\mathbf{B}$ velocity field, where $\mathbf{E}$ is the electric field due to the space charge and $\mathbf{B}$ is the external magnetic field. This interpretation provides simple explanations for the stability of round beams and for the development of spiral halos in elongated beams. By numerically solving the nonlinear advection equation for the beam density, we find that it is also in quantitative agreement with results obtained in PIC simulations.


I. INTRODUCTION
In recent years, novel uses of particle accelerators for nuclear energy, nuclear security, and medicine and material science applications have triggered the design, development, and construction of moderate size, high-intensity cyclotrons. The higher beam intensities expected in these new machines provide more powerful tools for these new endeavors, but also come with very stringent uncontrolled beam loss requirements. Satisfying these low beam loss criteria requires detailed knowledge of the beam dynamics in such machines. In particular, collective modes and instabilities associated with space-charge effects have to be taken into account and need to be studied and understood. Theory and simulation tools are being developed to that purpose. The computational effort is crucial because of the complexity of the problem at hand. One can hardly imagine a purely analytic study that would accurately account for the complex magnetic geometry of modern cyclotrons and at the same time properly model the nonlinear beam dynamics due to space-charge effects. However, analytical theory work is equally as important in order to provide an interpretation for the results from complicated, allinclusive simulations as well as to identify the basic contributing mechanisms.
Interestingly, the computational effort in high-intensity cyclotrons is almost entirely focused on particle-in-cell (PIC) methods [1][2][3][4]. PIC methods have a number of advantages. They are intuitive, naturally lead to parallel solvers, and can easily take exact magnetic geometries as inputs. Nowadays, PIC methods are commonly used to design, interpret, and corroborate experiments in which space-charge effects play a significant role in the beam dynamics [5][6][7]. However, precisely because of their conceptual simplicity, PIC codes tend not to yield as much insight into the physical mechanisms at work in a given observed phenomena as do continuum descriptions. Theoretical work based on continuum descriptions, whether kinetic or fluid, is often better suited to identify the key nondimensional parameters and to derive the typical scalings for the dependence of the quantities of interest on these nondimensional parameters. Analytic and numerical work based on continuum descriptions is therefore a necessary and useful complement to PIC simulations.
In this article, we illustrate the constructive interplay of experiments, PIC simulations, and theory in the study of a topic of high interest in beam dynamics in isochronous cyclotrons, namely the coherent longitudinal-radial vortex motion associated with the nonlinear interplay of radial and longitudinal space-charge forces. It has been observed in PIC simulations [3,[8][9][10]] that a coasting beam which is elongated in the longitudinal direction (the direction of propagation of the beam) develops a spiral galaxy shape due to space-charge effects. This phenomenon, which has sometimes been called the ''spiraling instability'' [8], can play a positive role when the spiral arms are lost after a few turns, at low beam energy, and the center of the spiral becomes a very compact stable beam [3,[10][11][12]. The spiraling of the beam can also play a negative role. It is thought to lead to beam breakup, as confirmed by experiment and PIC simulations [13,14]. At lower beam currents, the spiral arms can also grow over longer time scales, and therefore form a beam halo only after a large number of turns, at which point it has been accelerated to high energies. It is important here to note that beam spiraling in this context occurs in configuration space. This is in contrast to beam spiraling and halo formation more commonly studied as a phenomenon happening in phase space, often in the context of linear accelerators [15,16].
Until now, aside from PIC simulations, the theoretical explanations of beam spiraling have been based on singleparticle pictures [10,13,14,17]. In these single-particle descriptions, the evolution of the beam in the presence of space-charge forces is extrapolated from the modifications of the motion of single particles due to this force. In contrast, using a fluid model of the charged beam we present a new self-consistent nonlinear description of the coherent radial-longitudinal vortex motion and give an intuitive explanation for the formation of the spiral galaxy shape in an elongated charged particle beam. With a multiscale analysis that separates the time scale associated with the betatron motion and the time scale associated with the space-charge force, we show that in the frame rotating with the center of the beam, the evolution of the beam on the space-charge time scale is determined by the nonlinear advection of the beam density by the E Â B velocity. Here, E is the self-electric field, and B is the external cyclotron magnetic field. With this new fluid picture of the radial-longitudinal vortex motion, the observed behavior in experiments and PIC simulations is well understood, and beam spiraling has a simple explanation.
Imagine looking at a beam from a top view, aligned with the externally applied magnetic field, and consider first a round beam, whose charge density is only a function of its radius. In the frame rotating with the beam, the self-electric field of such a charge distribution is only along the direction of the beam radius, and the E Â B velocity vector in the radial-longitudinal plane is therefore purely perpendicular to the beam radius at every point. This implies that the advection in the E Â B field will have no effect on the radial distribution of the beam density: a round beam thus has a stationary distribution. In an elongated beam however, the cylindrical symmetry of the problem is broken since the density distribution depends both on radius and angle. The beam density distribution will thus be distorted by the E Â B velocity field. Since the E Â B velocity is largest at the extremities of the beam, the distortion will be largest away from the beam center, which leads to the formation of spiral galactic arms. This intuitive explanation is confirmed by simulations in which we numerically solve the nonlinear fluid equations, and we demonstrate quantitative agreement between our fluid description and PIC simulations.
The structure of the paper is as follows. In Sec. II, we describe the geometry of the problem and introduce the fluid description of the beam which we use for our analysis. In Sec. III, we carry out a multiscale expansion of the fluid equations, and derive the advection equation for the evolution of the beam density on the time scale associated with the space-charge force. We solve this equation numerically in Sec. IV for different cases of interest, and demonstrate both the formation of spiral galactic arms for elongated beams and the stability of round beams. We end with conclusions in Sec. V.

II. FLUID DESCRIPTION OF THE BEAM
One of the main motivations for this work is to offer an analytic derivation of the self-consistent beam vortex motion under the effect of space charge. Such a goal implies that we will have to make simplifications in order to make our analysis tractable. In making these simplifications, our guiding principle is to simplify the geometric aspects of the problem enough so that only a minimal amount of beam physics assumptions need to be made.

A. Geometry of the problem
To make our analysis more tractable, we restrict our study to the two-dimensional problem of a nonrelativistic coasting beam in the radial-longitudinal plane, immersed in a homogeneous magnetic field along the z axis, which coincides with the vertical direction. Throughout our analysis we will thus have B¼Be z with B¼constant, and @=@z ¼ 0 will hold. One implication of these assumptions is that the formation of spiral arms does not depend on relativistic effects or field gradient along the vertical direction.
Our choice to focus on the dynamics in the radiallongitudinal plane can be justified as follows. The physics of a non-neutral plasma immersed in a strong external magnetic field is highly anisotropic, and space-charge effects in the direction along the magnetic field couple only weakly with space-charge effects in the plane perpendicular to the magnetic field. In the vertical direction, space-charge defocusing leads to tune depression for the vertical betatron oscillations, which can be compensated by stronger vertical focusing. The dynamics in the radial-longitudinal plane is more complex and the subject of this article. The twodimensional nature of the problem is also confirmed by the results of a large number of PIC simulations, in which it was found that the details of the beam description along the vertical direction, i.e., along the direction of the main applied magnetic field, have a very limited influence on the longitudinal-radial vortex motion [3,9,18].
In the radial-longitudinal plane, i.e., the plane perpendicular to the magnetic field, it is most convenient to study space-charge effects in the frame moving with the coasting beam. The transformation from the laboratory frame to the moving frame is best understood in terms of the polar coordinates ðR; ÈÞ associated with the cyclotron geometry and shown in Fig. 1.
We consider a single-species intense ion beam, and for the simplicity of the notation, the ions are singly charged: q ¼ e. The ion mass is m, so that the cyclotron frequency is ! c ¼ eB=m. The change of reference frame to the frame rotating with the beam is given by the transformation where the primed quantities represent quantities in the moving frame. Using the relations in Eq. (1), we find In the next section, we will use these relations to express the fluid equations in the frame moving with the beam. Before doing so, we briefly discuss the two coordinate systems we will use in the local frame moving with the beam, which are shown in Fig. 2. The first system is a twodimensional Cartesian coordinate system ðx; yÞ with the origin at the beam center, and defined such that ðx; y; zÞ is a right-handed orthogonal coordinate system with the x axis coinciding with the radial direction, i.e., the main radius of the cyclotron, and the y axis coinciding with the longitudinal direction. This is the coordinate system also plotted in red in Fig. 1. In this coordinate system, the ion beam travels in the Ày direction. The second system is a polar coordinate system ðr; Þ centered at the center of the beam, defined such that ðr; ; zÞ is a right-handed orthogonal cylindrical coordinate system and such that x ¼ r cos and y ¼ r sin. It is equipped with the orthonormal basis ðe r ; e Þ, where the unit vectors e r and e are defined in terms of the Cartesian unit vectors e x and e y according to e r ¼ cose x þ sine y and e ¼ À sine x þ cose y . The Cartesian coordinate system will be mostly used for the numerical simulations in Sec. IV, and the polar coordinate system is the most appropriate for the physical interpretation of results. In terms of the polar coordinate system, we will call a ''round beam'' a beam whose density distribution only depends on the coordinate r, and we will call an ''elongated beam'' a beam elongated in the direction of propagation and whose density distribution depends on both r and .

B. Fluid description of the beam
At the beam temperatures and densities relevant to highintensity cyclotrons, the charged particle beam is a nonneutral plasma that can be considered collisionless [3]. Then, by taking the first two moments of the collisionless Vlasov equation, we obtain the following exact fluid equations involving the beam density n, velocity v, and pressure tensor P [19]: In Eq. (4), E is the self-electric field, which is determined by solving Beam center B is the external magnetic field B ¼ Be z since the selfmagnetic field is negligible in the nonrelativistic limit [20]. Using the relations in Eq. (2), we can rewrite Eqs. (3)- (5) in the frame rotating with the beam. We find where as before the primed quantities represent quantities in the moving frame. From now on we will always work in the frame rotating with the beam, so to simplify the notation we will drop the prime symbols in the remainder of the article. Note that in the moving frame, there is a minus in the magnetic field force. It is due to the Coriolis force.
In the frame rotating with the beam, the electrostatic approximation holds [3], so we write E ¼ Àr where is the electrostatic potential. Then, if we call a the characteristic size of the beam (its radius if the beam is circular), N 0 its peak density, and T 0 its peak temperature, renormalizing the fluid quantities according to leads to the following set of nondimensional fluid equations in the frame rotating with the beam: Here, d=dt @=@t þ v Á r and 2 ¼ ! 2 p =! 2 c , where ! p is the plasma frequency at the peak beam density, where D is the Debye length. 2 represents the strength of space-charge effects relative to the confining magnetic force, while 2 represents the strength of temperature effects relative to space-charge effects. The experimental regime of interest for highintensity beams is < 1, so that space-charge effects dominate over temperature effects.
It is clear from Eqs. (9)-(11) that the parameter determining the importance of space-charge effects on the behavior of the beam is 2 . This was observed in a slightly different context in Eq. (12) in [13]. A large majority of cyclotrons satisfy 2 1, and among these most machines even satisfy 2 ( 1. The parameter 2 can then be used to derive an approximate form for the pressure tensor P in the limit of small , as we will do in the next paragraph. It can also be used to perform a multiscale analysis of the resulting beam fluid equations which separates the fast motion associated with betatron oscillations from the slower motion due to space-charge and thermal forces. This is precisely what we do in Sec. III. Equations (9)-(11) cannot be solved as such since we are still facing the problem of closure of the hierarchy of moment equations [19]. However, a key point of our analysis is that with our simplified geometry, a single assumption we make concerning the beam physics is sufficient to address the closure issue. Specifically, we assume that the amplitude of the betatron oscillations in the frame moving with the beam is small compared to the characteristic size of the beam by the ratio , as we will discuss in more detail in Sec. III when we describe the multiscale expansion. In this regime, we prove in the Appendix, starting from the collisionless Vlasov equation, that the pressure tensor P must be gyrotropic to lowest order [21]. That is to say that P can be written as follows: where b ¼ B=B is equal to e z in our geometry, I is the unit 3-by-3 matrix, p ? is the pressure in the radial-longitudinal plane, and p k the pressure along the direction of the magnetic field. The divergence of the gyrotropic pressure tensor is readily calculated: Now, since b Á rb ¼ e z Á re z ¼ 0, in the radiallongitudinal plane the divergence of the pressure tensor takes the form of a pressure gradient to lowest order in : In the multiscale analysis presented in Sec. III, we will need to expand the beam fluid equations to second order in . Since in Eq. (10) the pressure term is multiplied by 2 , we only need an expression for the pressure tensor to zeroth order in . This is precisely what we have in Eq. (14). In the rotating frame, the appropriate fluid equations for the beam thus take the following form: where for the simplicity of the notation we have dropped the ? symbol for p since we only consider the twodimensional radial-longitudinal plane and p k never enters in the analysis. Equations (15)- (17) are the fluid equations we solve in the remainder of this article, using a multiscale analysis. Note that these equations are not fully closed in the sense that an equation for the evolution of the pressure p is missing. As we will show shortly, such an equation is not necessary. Indeed, when the pressure term can be expressed as an exact gradient in the momentum equation, pressure effects do not enter in the final equation for the evolution of the beam density on the space-charge time scale.

A. Multiscale expansion
The time scale associated with the betatron oscillations is 1=! c while, as shown in the previous section, the time scale associated with the expansion of the beam due to space-charge and temperature effects is ð 2 ! c Þ À1 . When 2 ( 1, as it is the case in most isochronous machines, the two time scales are well separated. As an illustration, the difference of time scales in the Small Isochronous Ring at Michigan State University is given by [22] so that the betatron time scale is more than 5 times shorter than the time scale associated with space-charge and temperature effects. We use this scale separation to perform a multiple time scale analysis of Eqs. (15)-(17) [23]. Specifically, each quantity Q is assumed to vary according to the different time scales as follows: Qðr; tÞ ¼ Qðr; t 0 ; t 2 ; t 4 ; . . .Þ ¼ Qðr; t; 2 t; 4 t; . . .Þ: (18) With this formal expansion, we have It is convenient for the rest of the calculation to separate the quantities Q into the sum of a rapidly oscillating partQ due to the betatron oscillations, and a slow monotonic evolution " Q due to space-charge and thermal effects: . . .Þ ¼Qðr; t 0 ; t 2 ; . . .Þ þ " Qðr; t 2 ; . . .Þ: Observe that, by construction, " Q does not depend on t 0 . We will show shortly thatQ is periodic in t 0 .
In the frame moving with the beam, the first nonvanishing contribution to the velocity is due to the betatron oscillations. As mentioned previously, we assume in our analysis that the amplitude of these oscillations is of order a. Since the betatron time scale is 1=! c , this means that the first nonvanishing contribution to the velocity is of order a! c . In other words, the velocity comes in first order in in our ordering and the appropriate expansion for the relevant physical quantities is the following: p ¼ " p 0 þ OðÞ; whereñ 1 andñ 2 represent the first and second order corrections to the density associated with the betatron motion, andṽ 1 andṽ 2 are the first and second order corrections to the velocity associated with the betatron motion.
It is important to see that focusing on beams for which the betatron oscillations are small by compared to the beam size is not inconsistent with the fact that ( 1 in most isochronous cyclotrons, which corresponds to the emittance-dominated regime. When ( 1, betatron oscillations dominate the beam dynamics, but do not alone determine the beam size. Initial conditions and energy spread also play a significant role. For instance, one can inject a beam into an isochronous cyclotron in such a way that particles on the outside have more velocity than particles on the equilibrium orbit, and therefore always remain on the outside of the beam. In the ideal limit, the beam is then in laminar flow. In this article, we consider situations in which the departure from laminar flow is small. We now introduce the expansion (21)- (24) in the set of equations given by Eqs. (15)- (17), and solve order by order in . We start with Poisson's equation, which we only need to zeroth order: Turning to Eq. (15), the mass conservation equation, we see that the first nontrivial equation arises in first order in and is given by This equation describes the evolution of the density on the fast time scale, due to the betatron oscillations. To next order, the mass conservation equation takes the form: Averaging this equation over the fast time scale, we get, because of the periodicity of the oscillating quantities in t 0 , @ " n 0 @t 2 þ r Á ðhñ 1ṽ1 i þ " n 0 " v 2 Þ ¼ 0; where we have introduced the notation Qdt 0 : Equation (28) describes the evolution of the bunch density on the slow time scale. This is where the influence of the electrostatic effects and the beam temperature will enter, and this is therefore all we need from the mass conservation equation. The purpose of the remainder of the article is to solve this equation, first by deriving expressions for hñ 1ṽ1 i and " v 2 in terms of zeroth order quantities, and by then solving the resulting equation numerically. The necessary expressions are obtained from the momentum equation. To first order in , we have and by averaging the second order momentum equation over the fast time scale as was done before, we find At this point, all the relevant equations have been derived. From Eqs. (26), (29), and (30) we can derive expressions forñ 1 ,ṽ 1 , and " v 2 in terms of zeroth order quantities and initial conditions, which we can then insert in Eq. (28) to obtain the desired equation for the evolution of the beam density on the space-charge time scale. This is done in the next subsection.

B. Solving the equations
We start the calculation with the description of the betatron motion to lowest order, i.e.ṽ 1 , and its effect on the bunch densityñ 1 .ṽ 1 is given by Eq. (29), which is just the equation of a particle immersed in a uniform magnetic field. Its well-known solution is v 1 ðr; t 0 ; t 2 ; . . .Þ ¼ u 1 cost 0 þ e z Â u 1 sint 0 ; where u 1 ðr; t 2 ; . . .Þ is the initial betatron velocity, i.e.ṽ 1 at time t 0 ¼ 0; 2; . . . . Now thatṽ 1 is known, Eq. (26) is most conveniently solved by introducing the displacement vector 1 ðr; t 0 ; t 2 ; . . .Þ defined bỹ Substituting forṽ 1 with 1 in Eq. (26), we find @ @t 0 ½ñ 1 þ r Á ð " n 01 Þ ¼ 0 which one can integrate to obtaiñ n 1 ¼Ñ 1 ðrÞ À 1 Á r " n 0 À " n 0 r Á 1 ; whereÑ 1 is the initial beam density. Equation (32) can also be integrated using Eq. (31), and we have 1 ¼ u 1 sint 0 þ u 1 Â e z ðcost 0 À 1Þ ¼ ðṽ 1 À u 1 Þ Âe z : (35) These results to first order in can now be used in Eq. (30) and (28). After a slightly lengthy calculation, we find so that the equation for the evolution of the beam density on the space-charge time scale is given by From the identity r " p 0 Â e z ¼ r Â ð " p 0 e z Þ we see that the pressure term vanishes in Eq. (37) and we finally have Equation (38), combined with Poisson's equation, forms the desired closed set of coupled equations, describing the evolution of the beam density under the effect of space-charge forces. It is interesting to note that temperature effects do not enter in this equation when the pressure tensor is gyrotropic, as it is with the ordering we chose for the amplitude of the betatron oscillations. This is the reason why we do not need to add an equation for the evolution of the pressure to the fluid equations (15)- (17). Equation (38) can be interpreted as the advection of the density profile in the velocity field E Â B=B 2 , the socalled E Â B velocity (there is no minus sign in front of r " 0 because the Lorentz force is in the opposite direction in the moving frame). This is in agreement with Eq. (1) in Ref. [24], and can be explained as follows. In the absence of accelerating gaps, the effect of the betatron oscillations on the density profile averages out to zero on the slow time scale. Thus, the density profile simply follows the slow, averaged motion of the ions. In a homogeneous magnetic field, in the presence of an electric field created by the charges in the bunch themselves, this slow motion is just the E Â B motion.
With this intuitive interpretation, it is easy to predict the early stages of the evolution of the beam for the density profiles discussed in the Introduction. The two situations of interest are illustrated in Figs. 3 and 4, which show the r " 0 Â e z velocity field due to a cylindrically symmetric density distribution and an elongated density distribution, respectively. If we start, as in Fig. 3, with a density profile which is cylindrically symmetric, i.e., such that the beam density is only a function of the distance from the center of the bunch [ " n 0 ¼ " n 0 ðrÞ], then the initial electrostatic potential due to the charges in the beam will also be cylindrically symmetric: " 0 ¼ " 0 ðrÞ. This implies that the initial electric field will be along the e r direction, E ¼ E r e r , and the CERFON et al.
Phys. Rev. ST Accel. Beams 16, 024202 (2013) 024202-6 convective velocity r " 0 Â e z will be along the e direction. We thus expect the density profile to rotate around the center of the bunch. Since this density profile is cylindrically symmetric, a rotation of the profile around the center will leave the profile invariant. In other words, the spacecharge forces in a cylindrically symmetric bunch in a purely drifting region are such that the bunch properties are kept constant. This result agrees with a similar result obtained by Kleeven et al. [25] with a different method and with a very different set of assumptions.
If, however, the initial density distribution is elongated in the direction of propagation of the beam as in Fig. 4, the cylindrical symmetry will be broken, and the E Â B velocity field will distort the beam. Since the electric field is strongest at the extremities of the beam, the distortion will be strongest in these places. This initial situation explains the formation of the spiral arms in the early stages of the beam evolution. In addition, since the electrostatic potential far away from the beam still retains cylindrical symmetry (far away, the beam is seen as a point charge), we expect the E Â B convection to act in such a way as to make the beam conform with this cylindrical symmetry. This is the reason for the formation of the often mentioned galaxy-like shape, and after more turns, the formation of a compact, stable cylindrically symmetric core. In the next section, we solve the system of coupled equations (38) and (39) numerically and confirm these predictions.

IV. NUMERICAL RESULTS
In order to solve Eqs. (38) and (39) numerically, we first rewrite them as follows: In Eq. (40), t represents the numerical time steps. Since they are defined in terms of ! c in the simulations, i.e., in terms of the fast time scale, 2 appears in the advection term. E x is the electric field in the radial direction, satisfying E x ¼ À@ " 0 =@x and E y is the electric field in the longitudinal direction, satisfying E y ¼ À@ " 0 =@y. We assume open boundary conditions for the electric field, so that E is simply due to the charges in the beam.
We solve Eqs. (40) and (41) by discretizing the problem on a regular Cartesian grid. The time integration of Eq. (40) is done with a second order backward difference scheme, except for the first time step which is calculated with a second order trapezoidal scheme. At each time step, we compute the partial derivatives of the density at each grid point with second order centered differences, assuming that the beam density is zero outside the boundary of the computational domain, which is always chosen to be larger than the maximum beam size. We also compute E x and E y at each time step on this grid by numerically solving the two Poisson's equations in Eq. (41). Specifically, we calculate E x and E y at interior points by inverting the standard finite difference second order Laplacian matrix, and use the two-dimensional free space Green's function for Poisson's equation to calculate the electric field at the computational boundary: 0 Â e z (magenta arrows) for a cylindrically symmetric density distribution. To compare with figures in Refs. [3,18], the beam's transport direction is along the negative y direction.
FIG. 4. Velocity field r " 0 Â e z (magenta arrows) for an elongated density distribution. To compare with figures in Refs. [3,18], the beam's transport direction is along the negative y direction. For the simulations we present in this article, we choose the same beam parameters as the ones used in [3] for the study of a 1 mA, 3 MeV coasting beam in PSI injector II. With these parameters, we find 2 $ 0:8. This value is at the upper limit of the validity of the multiscale analysis used in this work. Indeed, one expects the multiple scale analysis derived in Sec. III to be robust as long as < 1, but when is very close to 1, the accuracy of our approximate solution is not as good since it is obtained by neglecting terms of order 3 that are not much smaller than lower order terms. We still choose to apply our theory to the beam dynamics in PSI injector II because it allows us to compare our results with a large body of PIC simulations performed to investigate space-charge effects in this machine [3,8,9], some of which were obtained with very high performance numerical tools. As we will see, the agreement we find between our numerical results and the PIC simulations indicates that the high order terms we neglected in our analysis do not play a significant role in the radiallongitudinal vortex dynamics, even if 2 $ 0:8.
In agreement with [3], we start with an initial charge distribution given by

A. Stability of a round beam
We first demonstrate numerically that a round beam is completely stationary under the effect of space-charge forces, as we discussed in the previous section. We choose 2 x ¼ 2 y ¼ 2:52 mm and run our simulation until the beam has completed 40 turns. The results are shown in Fig. 5. It is clear in this figure that the charge distribution is not affected by the advection in the r 0 Â e z velocity field, and the charge distribution after 40 turns is identical to the initial charge distribution, as expected. In fact, the results shown in Fig. 5 can be seen as a proof of the accuracy of our numerical method. We can now confidently turn to the study of elongated beams and beam spiraling.

B. Elongated beams and beam spiraling
For the study of beam spiraling, we consider an elongated beam with 2 x ¼ 2:52 mm and 2 y ¼ 13:4 mm (15 phase width) as in [3]. The numerical results are shown in Fig. 6 (early times) and Fig. 7 (later times),  [3,18], the beam's transport direction is along the negative y direction. 2 ¼ 0:8 and 2 x ¼ 2 y ¼ 2:52 mm. To compare with figures in Refs. [3,18], the beam's transport direction is along the negative y direction. 2 ¼ 0:8, 2 x ¼ 2:52 mm, and 2 y ¼ 13:4 mm. where, as in Refs. [3,18], we have truncated the beams at 10% of the maximum charge density. By comparing Fig. 6 with the early stages of the beam evolution shown in Ref. [18], we observe that there is a strong similarity between our results and those obtained in simulations with the PICS code. We obtain a slightly stronger deformation than seen with PICS in Ref. [18]. This small discrepancy can be attributed to two factors. First, the magnetic geometry is different in the two simulations: the magnetic field in our simulation is uniform and in the z direction, while Adam used an analytic approximation of the magnetic field in the PSI injector II cyclotron. Second, the shape of the initial charge distribution, which we chose to be identical to that in [3], is slightly different. The good agreement between our results and the PICS simulations is not surprising since PICS is a twodimensional PIC code which relies on the separation between the betatron and space-charge time scales to reduce the computational complexity of the PIC simulations. The agreement also suggests that the details of the magnetic geometry do not play a strong role in the radiallongitudinal vortex motion due to space charges.
Comparing Figs. 6 and 7 with PICN results [18] and with OPAL-CYCL results [3], we conclude that there is qualitative agreement between all descriptions, with the formation of a compact stable core after 40 turns. Furthermore, as in [3], we see that a tenuous low density halo persists, even after 40 turns. However, there remains small differences between [3] and our study. First, the core that is formed after 40 turns is round in our simulations, unlike the slightly elongated core observed in [3]. Second, the beam spiraling effect appears to be stronger in our results. As before, part of the discrepancy can be explained by the difference in the magnetic geometry assumed in the different simulations. It is also a direct consequence of the twodimensional description that we have adopted in our study, which leads to a stronger self-electric field, and thus a larger r 0 Â e z velocity field. This effect has already been observed in the difference between the PICS and PICN codes [9]. Obtaining fully quantitative agreement with OPAL-CYCL on this problem would thus require a three-dimensional fluid treatment, and the inclusion of the details of the PSI injector II magnetic geometry.

V. CONCLUSION
In the vast majority of high-intensity cyclotrons, the betatron time scale remains much shorter than the time scale associated with space-charge and temperature effects. This separation of scales can be used advantageously in analytical and numerical studies of space-charge effects on beam transport. In this article we performed a multiscale analysis relying on this scale separation to study the radial-longitudinal vortex motion of an intense charged particle beam in a regime in which the pressure tensor is gyrotropic. As a result of this analysis, we have been able to show that the evolution of the charge distribution on the space-charge time scale is determined by the nonlinear advection equation of the charge distribution in the E Â B velocity field, where E is the self-electric field due to space charges and B is the externally applied magnetic field. The stability of intense round beams is easily understood in this light. A beam with a charge distribution which is only a function of the beam radius has an associated E Â B velocity field which is always perpendicular to the beam radius in the radial-longitudinal plane. Such a velocity field cannot modify the charge distribution. Beam spiraling in elongated beams is explained in a similarly simple way: in such beams the cylindrical symmetry of the charge distribution is broken, and the E Â B velocity field therefore modifies the charge distribution until the cylindrical symmetry is restored. The spiral shape is a natural consequence of the fact that E is strongest at the extremities of the beam.
Aside from the shape of the beam, the parameter that determines the strength of the spiraling effect is 2 ¼ ! 2 p =! 2 c . In cyclotrons with very intense beams and moderate magnetic fields, this parameter is relatively large ( & 1) and beam spiraling can play a positive role since a dense round stationary core can be formed after a few turns and the spiral halo can be stripped at low energies. In accelerators with less intense beams or higher fields, however, such that ( 1, beam spiraling can have a negative effect since it takes more turns to develop, and high energy beam halos can be formed in this process. The good agreement between our numerical simulations and PIC simulations demonstrates that the potential coupling between the radial-longitudinal motion and the vertical motion is not a key physics element that needs to be included to understand space-charge effects in the radiallongitudinal plane. This suggests that the simplified fluid equations (15)- (17), together with an equation for the evolution of the fluid pressure p, might be very well suited for numerical studies of the novel and interesting regime in which 2 > 1. In this regime the analytic multiscale expansion used in this work breaks down, but the 2D fluid equations presented here could still represent a significant reduction in complexity in numerical simulations as compared to PIC simulations. This topic is currently being investigated, with progress to be reported at a later date.

ACKNOWLEDGMENTS
The authors would like to thank A. Adelmann from the Paul Scherrer Institut for helpful discussions and for providing the beam parameters allowing comparison between this work and previous work on the beam dynamics in PSI injector II.

APPENDIX
In this Appendix we prove, starting from the collisionless Vlasov equation, that the pressure tensor P is gyrotropic to lowest order in when the amplitude of the betatron oscillations is small by compared to the size of the beam.
The collisionless Vlasov equation determining the evolution of the beam distribution function fðr; v; tÞ is Let us compare the relative strength of each term in this equation relative to the magnetic term, in the regime where v $ a! c : We see that all the terms are small by compared to the magnetic term, so that to lowest order in , Equation (A2) implies that to lowest order in the distribution function fðx; y; z; v x ; v y ; v z ; tÞ has the following particular dependence on v x and v y : fðx; y; z; v x ; v y ; v z ; tÞ ¼ f x; y; z; ; v z ; t þ OðÞ: (A3) Thus, to lowest order f is even in v x and v y , which implies that to the same order the off-diagonal components of the pressure tensor vanish: P ¼ p ? I þ ðp k À p ? Þbb þ OðÞ (A4) with p ? ¼ p xx ¼ p yy and p k ¼ p zz .