High intensity ion beams in rf undulator linac

The possibility of using a radio frequency undulator field to accelerate a high intensity ion beam in a linac is discussed. Such an accelerator can be realized using the periodical interdigital H-type resonator structure. The accelerating force is produced by an electric field which is a combination of two or more spatial harmonics, none of them being synchronous with the ion beam. The value of this force is proportional to the squared charge. The equations of motion in Hamiltonian form are derived by means of smooth approximation. The analysis of the 3D effective potential function allows finding the conditions of the beam focusing and acceleration. Two ways to increase ion beam intensity are considered: (i) to enlarge beam cross section; (ii) to neutralize the beam space charge by accelerating ions with opposite charge signs within the same bunch. The basic results are confirmed by a numerical simulation.


I. INTRODUCTION
Production of high intensity ion beams in a linear accelerator (linac) is a challenging task of contemporary nuclear physics.Such accelerators can be used as neutron generators and can also be employed in nuclear energetic, thermonuclear synthesis as well as in other applications.In a conventional radio frequency (rf) linac the beam is accelerated by a synchronous wave of the rf field.The radio frequency quadrupole (RFQ) structures are usually used in the buncher of the linac.The current transmission coefficient in the RFQ can be limited by the losses due to small channel aperture and powerful influence of the space charge fields.Therefore, the maximum proton beam current in the RFQ is 120 -150 mA [1].Another limitation of the RFQ structure is the low rate of the energy gain (usually not greater than 300-400 keV=m).
An alternative method of ion acceleration in electromagnetic fields without a synchronous wave was presented in [2].Some analytical studies have already been published in [3,4].The acceleration mechanism is similar to the acceleration mechanism in an inverse free electron laser (IFEL), where the electron beam is accelerated by a ponderomotive force [5,6].In IFEL the accelerated gradient equals the product of undulator field amplitude (B or E) and electromagnetic wave amplitude (E v ).In our case, the accelerating force is driven by a combination of two nonsynchronous waves which are supplied by two undulators.
Here we define an undulator as any structure producing a periodic-distribution electromagnetic field for which the phase velocity differs significantly from that of the beam.This innovative type of linac was called an undulator linear accelerator (UNDULAC).
There are three different types of undulators that can be used to design the required configuration of accelerating fields-magnetic (UNDULAC-M), electrostatic (UNDULAC-E), and rf undulator (UNDULAC-RF).As it has been shown in [7], one of the undulators must be of the rf type, the second one being, optionally, of magnetic, electrostatic or radio frequency types.The accelerating structure of UNDULAC-rf can be an interdigital H-type (IH) periodic resonator with drift tubes.It is simpler than RFQ and extends the limit of the beam current and the rate of energy gain; also it increases the transmission coefficient [8].
As the main factor limiting beam intensity in ion accelerators is the space charge force, there exist two ways to increase ion beam intensity: (i) to enlarge beam cross section and (ii) to use space charge neutralization.We will study these methods for UNDULAC-RF.

II. EQUATION OF MOTION IN SMOOTH APPROXIMATION
The first method of increasing output beam intensity involves the enlargement of the beam cross section.Ribbon beams can be used for this purpose.Acceleration of ribbon beams is possible in UNDULAC-E [2] and UNDULAC-RF [3,4].A relatively simple IH periodic resonator with two rows of electrodes connected to the vanes can be employed in UNDULAC-RF (Fig. 1).The rf field in a periodic structure can be represented as an expansion of spatial harmonics under the assumption that the period of the structure is a slowly varying function of the longitudinal coordinate.The analysis of the ribbon beam dynamics in a polyharmonic field is an intricate problem.Because of the fast oscillations of particles and the nonlinear dependence of the rf field components on the transverse coordinates in ribbon beam accelerators, it is not possible to apply the standard linear approximation of the field near the axis.In the quasistatic approximation, the potential of the rf field in the slot channel can be expressed as where h n 2n=D is the longitudinal wave number of the nth harmonic and is the phase advance per period of structure D. Function U n x; y satisfies the equation ( ? is the transverse component of the Laplacian) and determines the dependence of the potential versus transverse coordinates.We assume that the beam velocity differs greatly from the phase velocity for all rf field harmonics, n !=ch n , n 0; 1; . . ., but is close to the phase velocity of the combined wave which is a superposition of an nth and pth harmonic, v 2!=ch p h n .In that case, an analytical solution for the equation of motion is not possible using the traditional methods but can be provided by the means of the averaging method.The momentum and the coordinates of the particle can be represented as a sum of slowly varying and rapidly oscillating components: p p p, r r r.Supposing that the amplitude p of the fast momentum oscillations is much lower than the amplitude p of the slowly varying components, we can write the equation of the motion in a smooth approximation as it was described in [2 -4]: ( where U eff has the meaning of the effective potential function of the combined wave field in the resonator.This function describes the interaction of particles with the polyharmonic field and depends only on the slowly varying transverse coordinates R ? ; (where 2x= and 2y=) and the slowly varying phase of the particle in the field of the combined wave Here !t is dimensionless time, 2z= is dimensionless longitudinal coordinate, and ! are the wavelength and frequency of the rf field, h v h p h n =2 is the longitudinal wave number of the combined wave field.The combined wave field can accelerate ion beam as a ponderomotive force accelerates electrons in the IFEL.As it was shown in [2], the effective potential function U eff is proportional to h Ẽ Ẽi.In UNDULAC-RF U eff can be represented as where the first term, defines the longitudinal acceleration and transverse defocusing.Here e i q 2W 0 E i is the dimensionless amplitude of rf field harmonics, W o mc 2 is the rest energy, q is the charge of the particle, and n;v h n h v h v .From the equation of motion, one can find the Hamiltonian of the system in the smooth approximation: The equation of motion written in form ( 3) and ( 4) allows a simple analysis of both longitudinal and transverse motion of particles.Specifically, using this approach, one can find the transverse focusing conditions, the frequencies of the transverse and phase oscillations of particles, and the connection between the longitudinal and transverse motion.

III. ION BEAM DYNAMICS IN UNDULATOR LINAC
Let us consider the simplest case when only two spatial rf field harmonics are present in the system (n 0 and p 1) and the electric field E has only the longitudinal component at the axis ( 0).It is reasonable to introduce a quasisynchronous (reference) particle with the average velocity s which is equal to v .If the phase advance of the rf field per period of the structure is (-mode), then the longitudinal wave numbers are h 0 =D, h 1 3=D, h s 2=D, and the phase velocity of the combined wave is v D=.The expression for the effective potential function U eff can be written as 2e 0z e 1z ÿ e 0? e 1? sin2'g: Combining expressions ( 3) and ( 5), we get an equation for the particle velocity at the axis: where e eff e 0 e 1 = s is the effective amplitude of the combined wave.
To provide efficient bunching and acceleration of the beam, the undulator linac must consist of two subsections: a gentle buncher subsection (L b ) and an acceleration subsection (L a ).We consider the dynamics in an undulator linac with the assumption that the amplitudes e 0 and e 1 of the rf field harmonics, as well as phase ' s of the reference particle, are functions of the longitudinal coordinate in the gentle buncher subsection and are constant in the accelerating subsection.In UNDULAC-RF the modulation period of the beam is half of the rf field period.For the combined wave, two phase ranges should be considered: ÿ=2 < ' < =2 and =2 < ' < 3=2.In a gentle buncher, the reference particle phase, ' s , increases linearly from ÿ=2 to ÿ3=8 (in the first range), and from =2 to 5=8 (in the second range) and the rf field harmonic amplitudes continuously increase from the initial value of E 0;1;in E 0;1 z 0. In the acceleration subsection these parameters are constant values.Such separation allows significant improvement of the phase capture.The velocity of the reference particle varies along the accelerator axis according to Solving this equation, one can easily find the laws of variation of the period of the rf structure D s , as well as of the phase velocities of the zero rf field harmonic, 0;s 2 s , and the first rf field harmonic, 1;s 2 s =3.From Eq. ( 4), one can find a separatrix for the particles at the axis.
The vertical size () and the phase acceptance (') of the separatrix vary due to the fact that the amplitude of the combined wave and phase ' s increases.Figure 2 shows the variation of the vertical size of the separatrix of the combined wave as the function of the longitudinal coordinate z (curve 1) for E 1 =E 0 0:3 (a) and 0.6 (b).In this figure we also show the same curves obtained for the case when the beam interacts with the zero harmonic alone, i.e., under the assumption that velocity is close to 0;s (curve 2) and when the beam interacts with the first harmonic only, i.e., under the assumption that velocity is close to 1;s (curve 3).Curve 4 in Fig. 2 shows the change of the longitudinal velocity of the particle in the smooth approximation.The calculation was carried out for a beam of deuterium ions for E o 150 kV=cm, initial ion energy of W in 150 keV ( s;in 0:013), buncher subsection length L b 1 m, and acceleration subsection length L 2 m.The separatrices for the combined wave (curve 1), for the zero and the first rf field harmonics (curves 2 and 3, respectively) for 0:3 and 0.6 in the middle of the buncher subsection are represented in Fig. 3.
In the smooth approximation the phase trajectories of all accelerated particles are always within the separatrix of the combined wave.The particles can be lost (escape from the acceleration process) only if the adiabatic conditions are disturbed because of a rapid change in the phase velocity and amplitude of the combined wave during bunching.However, there are other reasons for losses.Figure 2 shows that, when the first harmonic amplitude is large (if E 0 150 kV=cm, must exceed 0.3), the separatrices of the combined wave and the first harmonic may partially or completely overlap in the middle of the buncher subsection.For a small ( < 0:3), fast oscillations of the particle velocity may by far exceed the vertical size of the separatrix (Fig. 4, curves 1-3), Figure 4 also shows the longitudinal velocity of particles with rapid oscillating components (curve 4).Even if the separatrices do not overlap, the phase trajectory of the particle with a large amplitude of the fast oscillation may be captured into the separatrix of the first harmonic.In these two cases described above, the velocity of particle may approach the phase velocity of the first harmonic of rf field.When this happens, the particle can be captured and accelerated by the first harmonic or even escape the acceleration process.
The averaging method fails in this case.Since these effects are observed for some coefficient , one may expect that an optimal harmonic amplitude ratio exists when the losses of the accelerated ions are minimal.The above expression for the 3D effective potential function (4) allows us to find (in smooth approximation) the condition for transverse stability of particles in the undulator linac.Indeed, in the absolute minimum of U eff , the conditions for longitudinal and transverse stability of the beam are established simultaneously.In the simple linear approximation (h x x 1 and h y y 1), the transverse focusing condition can be derived analytically.If the phase advance of the rf field per period of the structure is , then it is easy to show that focusing exists at any harmonic amplitude ratio .
Figure 5 shows the sections of effective potential function U eff 0; 0; for x 0 and y 0 (curve 1), U eff 0; y; 0 for x 0 and ' ÿ ' s 0 (curve 2), and U eff x; 0; 0 for y 0 and 0 (curve 3).The curves are produced for E o 150 kV=cm; 0:6; h x =h y 1=25; and ' s ÿ=2, s 0:013 [Fig.5(a), input of the buncher] and ' s ÿ3=8, s 0:018 [Fig.5(b), output of the buncher].In the buncher input, U eff has two minima and one maximum along the x axis.As the particle velocity grows, the local maximum smoothens out and U eff obtains one absolute minimum.The presence of the intermediate maximum does not lead to a considerable redistribution of the density of particles in the beam, since the frequency of transverse oscillations (across the width of the ribbon beam) is very small.
In this section we studied the single particle approach and the influence of the space charge field was not taken into account.But the space charge field decreases the depth of the well U eff and the value of the current transmission coefficient.The influence of the Coulomb interaction of ions will be studied entirely by means of a computer simulation.

IV. TWO ION BEAMS DYNAMICS IN UNDULATOR LINAC AND SPACE CHARGE NEUTRALIZATION
The other possible method of increasing beam intensity is the space charge neutralization in undulator linac.In a conventional rf linac (RFQ and drift tube linac), the intensity of the ion beam can be made twice as high by simultaneous acceleration of ions with opposite charge signs (H , H ÿ or D , D ÿ ) [9].The accelerating force is proportional to the charge of the ion and oppositely charged ions are bunched in the different phases of the accelerating wave.Two bunches (one with a positive and another one with a negative charge) become separated and weakly interact with each other after the initial part of the buncher.In this case the phase separation of the bunch is large and the space charge neutralization cannot be achieved.
In contrast to conventional linacs, the space charge neutralization can be realized in undulator linacs [10].This conclusion follows from (3) and ( 6).The rapidly oscillating components of the momentum and the coordinate of particle, p and r, have opposite signs and the effective potential function, U eff , depends on the square of particle charge, i.e., averaged motion of positively and negatively charged ions is identical.In UNDULAC the ions with opposite signs are accelerated within the same separatrix.Therefore, although the bunch centers oscillate in antiphase, special separation of the positive and negative ions is small.

V. BEAMDULAC CODE
Analytical methods cannot be used to completely investigate beam dynamics.Numerical simulation is necessary to verify the results of any analytical study and space charge field can only be accurately treated by a numerical simulation.We have been developing the BEAMDULAC code for these simulations since 1999 [8].The BEAMDULAC computes the beam dynamics in linear accelerators and transport lines.2D and 3D versions were developed for axisymmetric structures and for ribbon beams, respectively.
The well-known cloud-in-cell (CIC) method is utilized for an accurate treatment of the space charge effects in the BEAMDULAC code.That is especially important in the case of high-intense charged particle beam.The equation of motion for each particle is solved in the external and the interparticle Coulomb fields.The charge is deposited on the grid points using the CIC technique.To determine the potential of the Coulomb field, the Poisson equation is solved on the grid with the periodic boundary conditions at both ends of the domain in the longitudinal direction.The aperture of the channel is represented as an ideal conducting surface with a rectangular or a circular cross section.This allows consideration of the shielding effect, which is sufficiently important for transverse focusing of ribbon beams.The fast Fourier transform (FFT) algorithm is used to solve the Poisson equation on a 3D grid.The obtained Fourier series for the space charge potential can be analytically differentiated, and thus each component of the Coulomb electrical field can be found as a series with known coefficients.The Coulomb repulsion force is the main factor limiting the beam current in linac.In our code, the space charge field can be calculated with the same precision as the Coulomb potential without numerical differentiation.
The external fields in BEAMDULAC code can be represented by means of three different methods.The field can be represented analytically or as a series of spatial harmonics.The field amplitude is represented as a polynomial coefficient series in the second case.The code can also be used for beam dynamics study in ''real field'' which is defined on a 2D or 3D grid by electrodynamics simulation codes or experimental results.Time is used as an independent variable and the standard fourth-order Runge-Kutta method is applied for integration of the equation of motion.
The last modification of the code includes multi-ion beam dynamics.This allows computation of the Coulomb field for the dual beam and investigation of simultaneous motion of positive and negative ions.The Poisson equation is solved using the conventional FFT algorithm.The distribution of particles on a 3D grid is calculated first.Then the Fourier series coefficients for charge are defined and the algebraic equations connecting Fourier coefficients for the charge and the potential are solved on the grid.Finally, the Fourier synthesis of the Coulomb potential and the space charge field components was done.The modification of space charge distribution calculation and the algebraic equation were provided for simulation of self-consistent dynamic.The new BEAMDULAC-2B version can be used for simulation of beam dynamics in the case when the beam contains the ions with a different charge-to-mass ratio.

VI. NUMERICAL SIMULATION OF SINGLE D ÿ ION BEAM DYNAMICS IN UNDULAC-RF
The numerical simulation was done for a ribbon beam of deuterium D ÿ ions.The results of such a simulation are detailed in [8].Let us briefly consider the results of the dynamics simulation in UNDULAC-RF where longitudinal rf field for the mode is used.The simulation was done assuming the following parameters: initial energy of deuterium ions W in 100 keV ( in 0:01), length of accelerator channel is 2 m, the cross section of the accelerator channel is 2a 2b 0:8 20 cm 2 , the wave length is 2 m.The effective amplitude of combined wave was chosen to be equal to E 2 0 0:3230 kV=cm 2 .In this case the output beam energy is equal to 1 MeV.
The current transmission coefficient is equal to 95% for paraxial injected zero-current beam (2l 2t 1 0:04 cm 2 ).The optimal value of is equal to 0.25-0.4,which coincides with the value found analytically.This ratio can be easily realized.We have obtained that K t mightily decreased if the beam size is larger than the critical value 2l 2t 5 0:3 cm 2 ).The particle losses are caused by fast oscillations of particle phases and longitudinal velocities.Current transmission coefficient K t is equal to 85%-90% for beam size of 2l 2t 5 0:3 cm 2 .As it was found the optimal length of the bunching section is equal to half of the accelerator length [8].The rate of energy gain is equal to 700-800 keV=m in the accelerating subsection of UNDULAC.This value is 2 times as much as the rate of energy gain in RFQ.
Numerical simulation of single beam dynamics shows that the current limit in UNDULAC-RF is equal to I max 0:3 A and the beam cross section is equal to 2l 2t 5 0:3 cm 2 (see Fig. 6).The output normalized transverse emittance is 2 times as much as the input emittance.For example, E y z 0 1 mm mrad and E y z 2 m 1:8 mm mrad when the space charge field is taken into account.

VII. SIMULATION OF DUAL BEAM DYNAMICS
The study of dynamics for dual deuterium D and D ÿ beam was done using the new code BEAMDULAC-2B.The code modification was discussed above (see Sec.

FIG. 2 .
FIG. 2. (Color) The variation of the vertical size of the separatrix of the combined wave as a function of the longitudinal coordinate z (curve 1), zero (curve 2), and first (curve 3) rf field harmonics and of the longitudinal particle velocity in the smooth approximation (curve 4) for 0:3 (a) and 0.6 (b).

FIG. 4 .
FIG. 4. (Color) The effect of the first rf field harmonic to the beam dynamics.The variation of the vertical size of the separatrix of the combined wave (curve 1), zero (curve 2), and first (curve 3) rf field harmonics and the longitudinal velocity of the particles with a rapid oscillating component (curve 4) are shown.
V). Five thousand macroparticles and grid range 16 16 16 for Coulomb field calculation was used for numerical simulation.About 4000 time steps are necessary for simulation of beam dynamics in the UNDULAC-RF section.Let us represent some results of the simulation of two beam dynamics.The input and output dual beam parameters are shown in Fig. 7.The input and output normalized transverse emittances in x; x and y; y planes are shown in Figs.7(a) and 7(b).The oscillations of phases for the mass center are plotted in Fig. 7(c) for both particle types.The output parameters are shown by the blue points and solid black lines for D ions and by cyan points and dotted black lines for D ÿ .The initial beam transverse emittance is shown in Figs.7(a) and 7(b) by red points and a solid brown line.The output phase spectra for D and D ÿ ions are shown in Fig. 7(d).

Figure 8
shows the longitudinal phase spaces for different z coordinates and illustrates the beam bunching.It is clear from the figures that D and D ÿ ions are accelerating within the same bunch in UNDULAC as was proposed.In the phase space the trajectories for positive and negative ions are oscillating in the opposite directions.

FIG. 9 .
FIG. 9. (Color) The current transmission coefficient (a) and the transverse beam emittance (b) versus the total initial beam flux.