Frequency map analysis of a three-dimensional particle in the core model of a high intensity linac

We consider the dynamical properties of a particle-core model for a uniformly ﬁlled triaxial ellipsoid in a periodic lattice of a high intensity linac. The mismatched oscillation modes are analytically computed in the smooth approximation and are compared with the numerical results of a tracking program. The study of the phase space in the mismatched case is performed by the frequency map analysis. In particular, we can analyze the effect of the nonlinear resonances between the envelope modes and the single particle sincrobetatron frequencies. A chaoticity criterion based on the frequency map analysis allows one to compute the stability region around the beam core. An estimate of the transport and its enhancement due to mismatch is provided by tracking orbits at the border of the stability region.


I. INTRODUCTION
One of the big challenges of the new generation of high power proton linacs is the control of beam losses down to a very low percentage. These losses are associated with the presence of a beam halo, populated by very few particles at large distance from the average beam dimensions. It is generally recognized that one of the main mechanisms generating beam halo is given by the nonlinear single particle resonances driven by the space charge, with an enhancement related to mismatched beam [1].
The above-mentioned effects can be better understood considering the frequency map analysis (FMA) [2] used for a test particle, and the beam core is assumed to follow known dynamics (particle-core model). In a previous article [3] we simulated the 2D dynamics of a mismatched beam propagating in a FODO channel. This numerical experiment pointed out the effect of the strong resonance n 1 n 2 on the particle transport in the phase space and the advantage of using different transverse tunes. In a linac the particle dynamics is intrinsically 3D, being the three tunes comparable. We assume an ellipsoidal uniform charge distribution. Neglecting this particular case does not correspond to a self-consistent solution of the Poisson-Vlasov problem. In a previous article we applied the FMA to a cylindrical symmetric ellipsoidal beam bunch located in a solenoid focusing channel and with rf focusing effects but without acceleration [4]. In such a case the space charge term for a bunch with cylindrical symmetry is calculated defining a form factor [5].
In this article we extend the FMA to a generic triaxial ellipsoid: the introduction of a generalized form factor allows one to compute the envelope modes in the nonsymmetric case. In this way by using the FMA a parametric analysis of a magnetic lattice can be explored with the advantages of using three different tunes [6]. The dy-namics of a 3D particle-core model was also studied in Ref. [7] using a different approach. We work out the details of the analysis for a magnetic line which is formed by two quadrupoles and two rf gaps according to the parameters given in Table I: In the sequel we refer to this structure as the FOGODOGO cell. This cell corresponds to the period of the ISCL superconducting linac investigated in the Italian research program TRASCO [8], which studies the feasibility of an accelerator driven system. In the present approach we neglect the effects of the acceleration and consider a steady state beam at the average energy. In Fig. 1 the longitudinal and transverse beam envelopes are plotted along the period. For this linac the number of particles per bunch is 5 3 10 8 . The sensitivity needed in the simulations with these parameters is related to the acceptable losses of 1 W͞m, that corresponds to relative losses of 0.7 3 10 26 per meter.
The particle distribution out of the radio frequency quadrupole (i.e., at the linac injection) is supposed to be almost homogeneous, except for the queue mainly due to the Debye screen that extends approximately to 1.1 beam envelope in the three dimensions; the population in the queue is Ӎ 3 3 10 25 of the total particles in the bunch. In order to fulfill our requirement of 1 W͞m, less than 3% of these particles, subject to nonlinear space charge forces, should reach large amplitudes and impinge on the beam pipe. To understand the dynamical transport at large amplitudes due to nonlinear resonances, we focus the FMA on the study of dynamics of particles in the queue region.
The FMA is performed on different plane sections of the phase space for the 3D particle-core model and gives a global description of the dynamics where resonant and chaotic regions are clearly detected. In this way it is possible to see the effects of each resonance on the halo formation as well as estimate a "dynamic aperture" for the short-term stability of the test particle. Moreover, taking advantage of the particle-core model, we can concentrate on the dynamics of the particles in the unstable regions with significant statistics and low computational cost. This paper is organized as follows: In Sec. II we discuss the computation of the electric field for a general uniform ellipsoidal distribution introducing the generalized form factor. In Sec. III we define the envelope and single particle equations of the particle-core model, and we explicitly compute the envelope modes in smooth approximation. In Sec. IV we perform the FMA for the TRASCO ISCL lattice and we compare the results of our approach with tracking simulations of the particle-core model. In particular, we study the possible contribution to the halo formation of the particles in the region within 1.1 beam envelopes.

II. ELECTRIC FIELD OF AN ELLIPSOIDAL CHARGE DISTRIBUTION
We consider a charge distribution uniformly distributed into a 3-axial ellipsoid of the following equation: According to Newton's potential theory, the electric potential of the distribution reads where Q is the total charge and e 0 is the vacuum dielectric constant. If x is an internal point of the ellipsoid, the lower integration limit x is set x 0; otherwise, for an external point x x͑ x͒ where x͑ x͒ is the positive root of the equation Equation (3) defines the axis of a confocal ellipsoid passing through the point x. In the sequel we introduce the notation a i q a 2 i 1 x for the axis of the confocal ellipsoid. The electric field E is given by where ͑i, j, k͒ defines any permutation of the indices ͑1, 2, 3͒. In order to give a simpler form to the integral in Eq. (4) we change the variable y â i ͞ q a 2 i 1 l, and after some algebraic passages using the equalityâ 2 j 2â 2 i a 2 j 2 a 2 i Eq. (4) reads where we define the following form factor: As is well known, the electric field is linear inside the ellipsoid (1). In the case of a cylindrical symmetric ellipsoid, the form factor (6) reduces to the form factor reported in the literature f͑p͒ F͑p, p͒ [5]. Equation (5) can be interpreted in the following way: At any external point x, the electric field generated by an ellipsoidal uniform charge distribution is equivalent to the electric field generated by a confocal uniformly charged ellipsoid passing through the point x. The form factor (6) is a symmetric function with respect to the argument and satisfies the equality which is a consequence of the Gaussian theorem for the electric field.

III. 3D PARTICLE-CORE MODEL
We consider a test particle of charge e and mass m moving in the FOGODOGO cell. The particle is under the influence of the space charge force of a 3D uniformly charged ellipsoidal bunch. The beam is bunched by rf cavities with longitudinal electric field E 0 and wave number k 2p͞bl. In the reference frame of the laboratory the equations of motion read where 9 denotes the derivative with respect to the arc length of the reference orbit, r q x 2 1 1 x 2 2 is the distance from the reference orbit, b and g are the relativistic factors, and I 0 ͑z͒ and I 1 ͑z͒ are the modified Bessel functions of order 0 and 1; K͑s͒ and K 3 ͑s͒ ekE 0 ͑s͒͑͞mg 3 b 3 c 2 ͒ are the quadrupole and the longitudinal focusing strength, respectively, which are periodic functions of period equal to the cell length L. The parameter with I c pe 0 mc 3 ͞e 7.8MA, is the space charge parameter for a beam current I.
Equations (8) have to be coupled with the envelope equations to determine the functionsâ i ͑s͒. We define the edge emittances e i 5 q where ͗ ͘ is the average over to the charge distribution. By a direct calculation we have ͗x 2 i ͘ a 2 i ͞5, and the envelope equations read The systems (10) and (8) define the 3D particle-core model. The uniform ellipsoidal distribution is not a solution of the Poisson-Vlasov system because the corresponding stationary distribution in the phase space is singular; nevertheless, it allows one to keep the Hamiltonian character of the 3D particle-core model similar to the dynamics of the full particle system (Liouville problem). The envelope equations have a periodic solution defining the matched beam. Introducing the Poincaré map of system (10), the matched solution corresponds to an elliptic fixed point. An analytical approach to the study of an equilibrium solution of the system (10) is possible in the smooth approximation; i.e., one substitutes the focusing strengths with their average values K i n 2 0i ͞L 2 , where n 0i ͑i 1, 2, 3͒ are the zero current tunes. In this case the periodic orbit reduces to a fixed point a i of the system and it is possible to compute analytically both the envelope's frequencies and the corresponding eigenvectors. As a consequence, the dependence from the relevant parameters of the model can also be analytically studied. To compute the envelope frequencies we need to linearize Eqs. (10) and work out the derivative of the form factor (6): an explicit calculation gives where we introduce the function The function G͑p, q͒ has some properties resulting from Eq. (12) and the Hamiltonian character of the envelope equations ͑1 2 p 2 ͒G͑p, q͒ 1 ͑1 2 q 2 ͒G͑q, p͒ 3F͑p, q͒ 2 1 , Considering Da i a i 2 a i , which stands for the displacement with respect to matched solutions and Eq. (10), the linearized envelope equations in the smooth approximation read where n i are the depressed tunes. The symmetric matrix H ij can be explicitly computed according to where the diagonal elements are The right-hand side of Eq. (14) is a negative defined matrix that can be diagonalized to compute the eigenvectors and the eigenvalues through Cardano's formula. In the case of an axially symmetric ellipsoid ͑ā 1 ā 2 ͒ the matrix H simplifies and it can be written as a function of the form factor f͑p͒ F͑p, p͒ [6] where p ā 3 ͞ā 1 , and we defined g͑p͒ G͑p, p͒ The calculation of the eigenvalues and the eigenvectors of Eq. (14) is explicitly performed by using the quantities The eigenvalues read and the corresponding eigenvectors are where the mode mixing angle f is defined The eigenvalues and the eigenvectors computed in the smooth approximation are useful for the analysis of the Fourier spectrum of the orbits in the FOGODOGO cell.

IV. FREQUENCY MAP ANALYSIS OF MODULATED MAPS
For a given perturbed d degrees of freedom Hamiltonian system, the Kol'mogorov-Arnol'd-Moser theory [9] proves the existence of invariant tori in the phase space which are characterized by the frequencies n i ͑i 1, . . . , d͒ associated with the quasiperiodic orbits lying on them. Therefore a map exists from the set of invariant tori to the frequency space, the frequency map (FM). For an almost integrable system the FM is smooth. The theoretical results can be extended to perturbed symplectic maps even if the dynamics is externally modulated by a quasiperiodic signal. More precisely, we consider a symplectic map z n11 M ͑ z n ; n a͒ that depends periodically on the modulation frequencies a ͑a 1 , . . . , a s ͒ where z ͑x 1 , x 0 1 , . . . , x d , x 0 d ͒ are the canonical variables. It can be proven that if the map M is almost integrable there exist invariant tori in the extended phase space R 2d1s associated with the frequencies ͑ n, a͒ ͑ n [ R d ͒ that satisfy a Diophantine condition [10]. Since the modulation frequencies a are fixed, each invariant torus is characterized by the n frequencies and there still exists a FM from the set of invariant tori into a frequency space R d . We note that the dynamics of the modulated map (22) feels the effects of the nonlinear resonances between the internal frequencies n and the modulation frequencies a so that the density of resonances is increased in the phase space and the chaoticity may appear due to the phenomenon of resonances overlapping.
The basic idea of the FMA is the numerical computation of the FM in order to study its smoothness. If one considers a uniform grid of points in a phase space section and computes the frequencies n associated with each orbit passing through a grid point, one gets a picture of the FM of the system. The smoothness property of the FM implies that the grid should appear smoothly deformed in the frequency space if the orbits belong to invariant tori. But when the FM is extended to the resonant or chaotic regions of the phase space the points no longer define a grid in the frequency space. More precisely, in the case of a modulated symplectic map (22) the initial points belonging to a resonant region will be mapped on a resonant plane so that a resonant region in the frequency space is defined by a local maximum of point density centered at the resonant plane (23), whose amplitude is proportional to the width of the resonance. The points corresponding to a chaotic region will appear as a fuzzy cloud due to the high sensitivity of the dynamics to the initial condition and to the finite precision of the numerical computation. As a consequence, the FMA provides information on the location of regular, resonant, and chaotic regions in the phase space and it can be used to study the effect of reso-nance overlapping on the stability of particle dynamics. The Poincaré map of the FOGODOGO cell in the presence of a mismatched beam is a modulated symplectic map with a equal to the envelope frequencies. In Fig. 2 we show the histogram of a number of points in the frequency space for a FM computed by using a grid of 5 3 10 3 points along the x 1 axis for the FOGODOGO cell in the presence of a matched (left-hand picture) and a 10% mismatched (right-hand picture) beam. The local peaks of the density correspond to the main resonances, whereas sudden local changes in the density are features of the chaotic regions. We observe that the appearance of resonances between the particle frequency n 1 and the envelope frequencies a enlarges the chaotic region. The spatial step of the initial grid defines the scale of the minimal structure of the phase space that can be detected by the FMA.
The problem of numerically calculating the FM for a modulated symplectic map is solved if one can identify the particle frequencies n in the Fourier spectrum of an orbit. This can be performed by projecting the orbit in the coordinate plane ͑x j , x 0 j ͒ ͑j 1, . . . , d͒ and looking for the frequency corresponding to the main Fourier component; this procedure is correct if both the nonlinear effects and the amplitude of the external modulation do not perturb the linear motion too much. This is the case of the Poincaré map of the considered FOGODOGO cell in the presence of a mismatched beam: in Fig. 3 we plot the fast-Fourier transform (FFT) of the projection of an orbit on the ͑x 1 , x 0 1 ͒ coordinate plane in the case of a matched (left-hand side) and 10% mismatched beam for a test particle near the border of the beam core. Even if in the mismatched case the Fourier spectrum is much richer due to the presence of envelope frequencies a, the particle frequency n 1 is clearly identified by the maximal amplitude in the spectrum.
The numerical computation of the FM needs an accurate evaluation of the frequencies for a great number of points. This task can be performed if one develops a very efficient algorithm to compute the frequencies of a regular orbit with a relatively small number of iterations. We used a Hanning filter together with an interpolation algorithm on the FFT [11]. For regular orbits (quasiperiodic orbits) this method allows a precision proportional to 1͞N 4 , where N is the iteration number and has the same computational load of the FFT [11]. The iteration number can be related to the weakness of chaos that is detected in the phase space by FMA.

V. NUMERICAL METHODS AND RESULTS
We integrate the particle-core equations (8) and (10) by using a leap-frog symplectic scheme which alternates linear transformations, kick maps for space charge, and nonlinear rf forces. 100 kicks per period were used. The form factor F͑p, q͒ is computed by a linear interpolation on a grid of points for the value of p, q. This method allows one to increase the velocity with a relative precision equivalent to 10 25 . The matched solution for the envelope equations is calculated by a bisection method.
We checked the tracking program by comparing the results for the envelope frequencies and the particle tunes analytically computed in the smooth approximation with the numerical values of the actual FOGODOGO cell. In Fig. 4 we plot the envelope frequencies a i as a function of the space charge parameter m. We observe that the difference between the smooth approximation and the actual lattice increases as we increase the beam intensity; this is due to the contribution of the beam oscillations in the presence of the space charge force. Table II reports the particle frequencies and the envelope frequencies corresponding to the nominal case of Table I. As is generally understood, the halo formation for a mismatched beam is mainly driven by the resonances 2n i a j ͑i, j 1, 2, 3͒ [1]. In Fig. 5 we plot the ratios n 1 ͞a j and n 01 ͞a j (left-hand picture) and the ratios n 3 ͞a j and n 03 ͞a j (right-hand picture) for each envelope frequency a j ͑j 1, 2, 3͒. If the value 1͞2 is in the area between the two curves, we expect the resonance to be excited in the phase space. We observe that for high current values, the 1͞2 resonance is crossed by each particle frequency, but we also need information about the position and the strength of the resonances in the phase space in order to prove the presence of large chaotic regions resulting from resonances overlapping in the mismatched cases [7]. Moreover, the above-mentioned curves do not take into account the nonlinear effects given by rf fields. Figure 6 shows the nonlinear tune shifts along the transverse (left) and the longitudinal (right) axis in presence or absence of space charge effect as a function of the distance from the beam core. We can see that only at large amplitude (5 times the beam envelopes), where the contribution of  the space charge force is negligible, the rf nonlinearity introduces a small negative tune shift. As a consequence we can disregard the effect of the rf nonlinear field on the beam core. When one considers a mismatched beam, a direct visualization of the phase space of the Poincaré map of the FOGODOGO cell is not useful since the excitation of the envelope modes continuously modulates the phase space structure. In such cases a description of the phase space is achieved by using the FMA. To perform the FMA we compute the Poincaré map of the lattice at the middle of the focusing quadrupole ͑ x 0 0͒; then we consider a family of orbits whose initial conditions are distributed on a uniform grid on a plane section of the bunch. We use the ͑x 1 , x 2 ͒ and the ͑x 1 , x 3 ͒ sections and 10 000 points distributed on concentric ellipses around the beam core. We analyze a region of the phase space covering up to 3 times the beam envelopes. Each orbit has been iterated for 2048 periods divided into two sets of 1024 points. We sepa-rately compute the particle tunes for each set and evaluate the norm of the difference D n of the tunes. This difference is used to apply a chaoticity criterion [12], to distinguish between regular and chaotic orbits. Indeed in the case of regular orbits which have a stationary Fourier spectrum, the accuracy of our algorithm would imply NDn ø 1. On the contrary, for the chaotic orbits the Fourier spectrum is not stationary and the criterion (24) allows one to detect the fluctuations between the frequencies of the first and the second set of iterations. The choice (24) for the threshold in the chaoticity criterion is related to the limited number of periods ͑Ӎ10 2 ͒ that define the linac lattice. A weak chaotic orbit could satisfy the criterion (24), but the diffusion time would be much longer than the number of periods in the linac. Moreover, the request of very low tolerance in the losses (#10 26 the total number of particles in a bunch) can be controlled by the density of grid points in the phase space sections. A Monte Carlo tracking with 10 5 particles has not detected any orbit that performs a fast diffusion in a regular region according to the criterion (24) for the chosen density [12]. Therefore we can define the border for the stable regions with a relatively small number of points ͑Ӎ5000͒ in the initial grid. The diffusion in the chaotic regions can be measured by a tracking that takes advantage of the information of the FMA. The frequencies are also used to perform the FMA. In Fig. 7   low order resonances in the frequency space whose number is increased due to the presence of envelope frequencies. Because of the resonance overlapping phenomenon the chaotic region at the border of the beam core is clearly enlarged. The plots of the initial conditions (bottom parts), corresponding to resonance regions, give a picture of the chaotic area around the beam core that allows the orbits to diffuse up to large amplitudes.
The tracking results confirm that the resonance overlapping is one of the principal causes of the reduction of the stability region around the beam core. We compute the dynamic aperture by evaluating the number of regular orbits in the initial grid of points that satisfy the criterion (24). We consider a sequence of transverse elliptical shells around the beam core with 100 points on each shell. In Fig. 9 we plot the number of regular points for each shell as a function of the initial radius R 0 normalized to the beam envelopes, The sudden drop of the number of regular orbits defines the dynamic aperture, which is larger in the matched case than in the mismatched case. The FMA together with the application of the stability criterion gives detailed information  on the effects of the nonlinear resonances and the extension of the chaotic regions. We make a first check of the FMA by tracking 10 5 particles randomly distributed in a transverse annulus of thickness 10% of the beam envelopes for the mismatched case. We have defined the transverse maximum radius as a function of the iteration number n, where the "max" is taken of the point distribution. The results are reported in Fig. 10: in the left-hand part we plot the maximum transverse amplitude after 10 3 iterations; it shows a halo formation up to an amplitude of 2.5 beam envelopes. In the right-hand part we plot the initial conditions of the orbits which can contribute to halo formation ͑R . 1.4͒ together with the initial conditions satisfying the stability criterion (24) (black dots). The two sets of points  9. (Color) Number of regular orbits on transverse shells around the beam core as a function of initial radius normalized to the transverse beam envelope size: each shell contains 100 initial conditions. The left-hand picture refers to the matched case, whereas the right-hand picture refers to the 10% mismatched case.  1.4) to detect the chaotic particles that can contribute to the halo formation. In the right-hand graph, we plot the initial conditions corresponding to the orbits in the chaotic region (small red dots) together with the profile of the beam (blue line) and the threshold ellipse R . 1.4. We plot with different markers the orbits that satisfy R . 2 (yellow dots) to show that they are spread out in the unstable region due to the chaotic character of the dynamics. In the right-hand graph, we also plot the points satisfying the stability criterion (24) (black dots) which do not overlap with the unstable region.
definitely do not overlap: this fact proves that the FMA is able to detect the stable regular regions. We also compute R͑n͒ on a population of 5 3 10 4 particles randomly distributed on five ellipsoidal shells of thickness 0.02 times the beam envelope in each direction; we consider 500 iterations for each orbit. In Fig. 11 we plot R͑n͒ for each shell of initial conditions. Two different effects can be observed: in the case of regular orbits a small increase of the R value is given by the smear of the orbits; in the case of chaotic orbits, R becomes a stepwise function due to the presence of partial topological barriers in the phase space, at the borders of fast diffusion regions. In the matched case we see the effect of chaoticity in the last two external shells. This shows the existence of a transverse regular region of order 1.08a 01 around the beam core.
In the mismatched case (10% of the quadrupole mode), R reaches large values (Ӎ2 beam envelopes) after a few hundred iterations in the second shell (between 1.02 and 1.04 beam envelopes). These results show that the mismatch of the beam drastically reduces the stable region around the beam core and enhances the halo formation phenomenon. This effect can be explained by the resonance overlapping in the phase space detected by the FMA. In Fig. 12, in order to identify the most unstable regions around the beam core, we plot the initial conditions corresponding to the maximum value of R͑n͒ for n # 500.  The isolated points correspond to regular orbits, whereas the cluster of points indicates the presence of chaotic regions where R has a sensitive dependence on the initial conditions. We note that the location of the most unstable regions can be far from the sections ͑x 1 , x 2 ͒ and ͑x 1 , x 3 ͒ used in the FMA. This gives the warning that two plane sections could not be enough for an exhaustive representation of the phase space in a 3D particle in core model. We compare the dynamic aperture of the stability criterion in the transverse section (see Fig. 9) with the dynamic aperture of the tracking results shown in Fig. 11. The FMA estimate is slightly optimistic in the matched case since we consider the ͑x 1 , x 2 ͒ plane section only, and we neglect the longitudinal-transverse coupling. Conversely in the mismatched case, the FMA estimate of the dynamic aperture is in accordance with the 3D tracking results since the main coupling through the envelopes frequencies is included.

VI. CONCLUSION
The FMA turns out to be a powerful tool to study the phase space of the 3D particle-core model in the mismatched case. The resonances excited in phase space are clearly detected and their extension can be measured. Even if we limit the analysis to the plane sections of the phase space, the FMA can be used to get information on the global behavior of the 3D dynamics. A single chaoticity criterion allows one to distinguish between chaotic and regular regions; the diffusion in phase space can be accurately measured by using a tracking code for the orbits in the chaotic region.