Approximate periodically focused solutions to the nonlinear Vlasov-Maxwell equations for intense beam propagation through an alternating-gradient field configuration

This paper considers an intense non-neutral ion beam propagating in the z direction through a periodic-focusing quadrupole or solenoidal field with transverse focusing force, 2 kx s xêx 1 ky s yêy , on the beam ions. Here, s bbct is the axial coordinate, gb 2 1 mbc is the directed axial kinetic energy of the beam ions, and the (oscillatory) lattice coefficients satisfy kx s 1 S kx s and ky s 1 S ky s , where S const is the periodicity length of the focusing field. The theoretical model employs the Vlasov-Maxwell equations to describe the nonlinear evolution of the distribution function fb x, y, x0, y0, s and the (normalized) self-field potential c x, y, s in the transverse laboratory-frame phase space x, y, x0, y0 . Here, Ĥ x, y, x0, y0, s 1 2 x02 1 y02 1 1 2 kx s x2 1 ky s y2 1 c x, y, s is the (dimensionless) Hamiltonian for particle motion in the applied field plus self-field configurations, where x, y and x0, y0 are the transverse displacement and velocity components, respectively, and c x, y, s is the self-field potential. The Hamiltonian is formally assumed to be of order e, a small dimensionless parameter proportional to the characteristic strength of the focusing field as measured by the lattice coefficients kx s and ky s . Using a third-order Hamiltonian averaging technique developed by P. J. Channell [Phys. Plasmas 6, 982 (1999)], a canonical transformation is employed that utilizes an expanded generating function that transforms away the rapidly oscillating terms. This leads to a Hamiltonian, H X̃ , Ỹ , X̃ 0, Ỹ 0, s 1 2 X̃ 02 1 Ỹ 02 1 1 2 kf X̃2 1 Ỹ 2 1 c X̃ , Ỹ , s , correct to order e3 in the “slow” transformed variables X̃, Ỹ , X̃ 0, Ỹ 0 . Here, the transverse focusing coefficient in the transformed variables satisfies kf const, and the asymptotic expansion procedure is expected to be valid for a sufficiently small phase advance (s , p 3 60±, say). Properties of axisymmetric beam equilibrium distribution functions, F b H 0 , with ≠ ≠s 0 ≠ ≠Q, are calculated in the transformed variables, and the results are transformed back to the laboratory frame. Corresponding properties of the periodically focused distribution function fb x, y, x0, y0, s are calculated correct to order e3 in the laboratory frame, including statistical averages such as the mean-square beam dimensions, x2 s and y2 s , the unnormalized transverse beam emittances, ex s and ey s , the self-field potential, c x, y, s , the number density of beam particles, nb x, y, s , and the transverse flow velocity, Vb x, y, s . As expected, the beam cross section in the laboratory frame is a pulsating ellipse for the case of a periodic-focusing quadrupole field or a pulsating circular cross section for the case of a periodic-focusing solenoidal field.


I. INTRODUCTION
Periodic focusing accelerators and transport systems [1][2][3][4][5] have a wide range of applications ranging from basic scientific research to applications such as heavy ion fusion, tritium production, spallation neutron sources, and nuclear waste treatment, to mention a few examples [6][7][8][9].Of particular interest, at the high beam currents and charge densities of practical interest, are the combined effects of the applied focusing field and the intense selffields produced by the beam space charge and current on determining detailed equilibrium, stability, and transport properties [1].Through analytical studies based on the nonlinear Vlasov-Maxwell equations, and numerical simulations using particle-in-cell models and nonlinear perturbative simulation techniques, considerable progress has been made in developing an improved understanding of the collective processes and nonlinear beam dynamics characteristic of high-intensity beam propagation in periodicfocusing and uniform-focusing transport systems .However, despite the extensive literature on intense beam equilibrium and stability properties, until the present paper, the Kapchinskij-Vladimirskij (KV) beam equilibrium [10,11], including its recent generalization to a rotating beam in a periodic-focusing solenoidal field [21][22][23], has been the only known periodically focused equilibrium solution to the nonlinear Vlasov-Maxwell equations for an intense beam propagating through an alternating-gradient quadrupole or solenoidal field configuration.While allowing for high space-charge intensity, the KV distribution is nonetheless of very limited practical interest, particularly because the (monoenergetic) distribution function has a highly-inverted (and unphysical) distribution in phase space and the corresponding density profile is exactly uniform in the beam interior.
It is, therefore, very important to develop a framework based on the nonlinear Vlasov-Maxwell equations [12,21] that is able to investigate the equilibrium and stability properties of a far more general class of periodically focused beam distribution functions.In a recent calculation [34], Channell has developed a third-order Hamiltonian averaging technique for investigating solutions to the nonlinear Vlasov-Maxwell equations for systems subject to a periodic external force.Following the Von Zeipel procedure, the formalism [34] uses a canonical transformation given by an expanded generating function to transform away the rapidly oscillating terms [35][36][37][38] and end up with a Hamiltonian H that depends only on "slow" variables.The purpose of the present analysis is to apply this averaging technique to intense beam propagation through a periodicfocusing lattice.The asymptotic expansion procedure is expected to be valid [34] for sufficiently small phase advance (s & 60 ± , say).
To briefly summarize, the present analysis considers a high-intensity non-neutral beam of positive ions (with charge 1Z b e and rest mass m b ) propagating in the z direction with characteristic average axial momentum g b m b b b c and directed kinetic energy ͑g b 2 1͒m b c 2 .The beam propagates through an applied field that produces a transverse focusing force, 2͓k x ͑s͒x êx 1 k y ͑s͒y êy ͔, on the beam particles.Here, V b b b c const is the average axial velocity, g b ͑1 2 b 2 b ͒ 21͞2 is the relative mass factor, c is the speed of light in vacuo, s b b ct is the axial coordinate, the ion motion in the beam frame is assumed to be nonrelativistic, and the lattice functions, k x ͑s͒ and k y ͑s͒, have axial periodicity length S const.Both the cases of a periodic-focusing quadrupole field [Eq.(9)] and a periodic-focusing solenoidal field [Eq.(11)] are considered in the present analysis.Furthermore, the analysis assumes negligible axial momentum spread, and the starting point is the nonlinear Vlasov-Maxwell equations (3) and (4) for the distribution function f b ͑x, y, x 0 , y 0 , s͒ and (normalized) self-field potential c͑x, y, s͒ in the transverse phase space ͑x, y, x 0 , y 0 ͒ in the laboratory frame [12,21].Here, the Hamiltonian for single-particle motion in the laboratory frame is given in dimensionless units by [Eq.( 6)] Ĥ͑x, y, x 0 , y 0 , s͒ where k x ͑s 1 S͒ k x ͑s͒ and k y ͑s 1 S͒ k y ͑s͒ are the (oscillating) lattice functions.The Hamiltonian Ĥ is formally assumed to be of order e, a small dimensionless parameter proportional to the characteristic strength of the focusing field as measured by the lattice coefficients k x ͑s͒ and k y ͑s͒.The organization of this paper is the following: The assumptions and theoretical model are summarized in Sec.II, including the nonlinear Vlasov-Maxwell equations for the distribution function f b ͑x, y, x 0 , y 0 , s͒ and self-field potential c͑x, y, s͒ in the laboratory frame.In Sec.III, we make use of Channell's third-order Hamiltonian averaging technique [34] to transform from laboratory-frame variables ͑x, y, x 0 , y 0 ͒ to a new Hamiltonian H ͑ X, Ỹ, X0 , Ỹ 0 , s͒ in the slow variables ͑ X, Ỹ, X0 , Ỹ0 ͒ correct to order e 3 .The formalism employs a canonical transformation given by an expanded generating function to transform away the rapidly oscillating terms [35][36][37][38].This leads to a Hamiltonian in the transformed variables of the form [Eq. ( 79)] H ͑ X, Ỹ , X0 , Ỹ 0 , s͒ where k f const.Of course, an important by-product of the generating function analysis is the determination of the coordinate transformation that relates the laboratory-frame variables ͑x, y, x 0 , y 0 ͒ to the slow variables ͑ X, Ỹ, X0 , Ỹ0 ͒ [Eqs.(87)-( 90)].The major simplification associated with transforming to the slow variables ͑ X, Ỹ, X0 , Ỹ 0 ͒ is immediately evident from the expression for H ͑ X, Ỹ, X0 , Ỹ0 , s͒.In particular, the focusing coefficient k f is both constant (independent of s) and isotropic in the transverse plane.This should be contrasted with the expression for the Hamiltonian Ĥ͑x, y, x 0 , y 0 , s͒ in the laboratory frame, where the focusing coefficients k x ͑s͒ and k y ͑s͒ are rapidly oscillating functions of s.In Sec.IV, following a discussion of the nonlinear Vlasov-Maxwell equations for F b ͑ X, Ỹ, X0 , Ỹ 0 , s͒ and c͑ X, Ỹ, s͒ in the transformed variables, we present several examples of axisymmetric equilibrium solutions, i.e., distribution functions F 0 b ͑H 0 ͒ with ≠͞≠s 0 and ≠͞≠Q 0, corresponding to beam equilibria with circular cross sections in the transformed variables [12,21].Of particular note is the class of distribution functions that satisfy ≠F 0 b ͑H 0 ͒͞≠H 0 # 0, which can be shown to be stable [25,26].Finally, in Sec.V, we exploit the inverse coordinate transformation, X͑x, y, x 0 , y 0 , s͒, Ỹ ͑x, y, x 0 , y 0 , s͒, etc., to determine properties of the periodically focused distribution function f b ͑x, y, s 0 , y 0 , s͒ in the laboratory frame correct to order e 3 , consistent with the class of constant-radius circular cross-section beam equilibria F 0 b ͑H 0 ͒ in the transformed variables.A wide range of important physical quantities are determined, including the distribution function f b ͑x, y, x 0 , y 0 , s͒; statistical averages such as the transverse mean-square beam dimensions, ͗x 2 ͘ ͑s͒ and ͗ y 2 ͘ ͑s͒, and the unnormalized transverse emittances, e x ͑s͒ and e y ͑s͒; and macroscopic properties such as the number density of 074401-2 074401-2 beam particles, n b ͑x, y, s͒ R dx 0 dy 0 f b ͑x, y, x 0 , y 0 , s͒, and the self-field potential, c͑x, y, s͒, etc.
Finally, in the third-order averaging technique developed in Sec.III, it should be emphasized that the Hamiltonian is formally assumed to be of order e, a small dimensionless parameter proportional to the characteristic strength of the focusing field [see Eqs.(6) and (15)] as measured by the lattice coefficients k x ͑s͒ and k y ͑s͒.To assure transverse confinement of the beam particles, the space-charge potential c͑x, y, s͒ in Eq. ( 6) is, of course, smaller than or comparable in size to the applied focusing potential, ͑1͞2͒ ͓k x ͑s͒x 2 1 k y ͑s͒y 2 ͔, and the kinetic energy contribution, ͑1͞2͒ ͑x 02 1 y 02 ͒, is allowed to be comparable in size to the applied focusing potential in the sense of a maximal ordering analysis.In this regard, treating the single-particle Hamiltonian to be of order e ø 1, where e is proportional to the focusing-field strength, is similar to the assumption made in standard analyses of the particle dynamics in intense charged particle beams at moderate values of phase advance [1][2][3][4][5].For completeness, in Sec.V D we provide a semiquantitative estimate of the range of validity of the asymptotic analysis in Secs.III and IV by relating the small parameter e to the focusing-field strength and the phase advance for the case of a sinusoidal quadrupole focusing lattice, k q ͑s͒ kq sin͑2ps͞S͒.

II. VLASOV-MAXWELL DESCRIPTION AND BASIC ASSUMPTIONS
In the present analysis, we consider a thin, intense non-neutral ion beam with characteristic radius r b and average axial momentum g b m b b b c propagating in the z direction through a periodic focusing field with axial periodicity length S. Here, r b ø S is assumed, ͑g b 2 1͒m b c 2 is the directed axial kinetic energy of the beam ions, g b ͑1 2 b 2 b ͒ 21͞2 is the relativistic mass factor, V b b b c is the average axial velocity, 1Z b e and m b are the ion charge and rest mass, respectively, and c is the speed of light in vacuo.The axial momentum spread of the beam ions is assumed to be negligibly small, and the ion motion in the beam frame is assumed to be nonrelativistic.We introduce the scaled time variable s b b ct and the (dimensionless) transverse velocities x 0 dx͞ds and y 0 dy͞ds.Then, within the context of the assumptions summarized above, the nonlinear beam dynamics in the transverse laboratory-frame phase space ͑x, y, x 0 , y 0 ͒ is described self-consistently by the nonlinear Vlasov-Maxwell equations for the distribution function f b ͑x, y, x 0 , y 0 , s͒ and the normalized self-field potential c͑x, y, s͒ Z b ef͑x, y, s͒͞g 3  b m b b 2 b c 2 , where f͑x, y, s͒ is the electrostatic potential.For a thin beam ͑r b ø S͒, we take the applied transverse focusing force on a beam particle to be of the form where ͑x, y͒ is the transverse displacement from the beam axis and the s-dependent focusing coefficients satisfy where S const is the axial periodicity length.The Vlasov-Maxwell equations for f b ͑x, y, x 0 , y 0 , s͒ and c͑x, y, s͒ can then be expressed as [12,21] and Here, n b ͑x, y, s͒ R dx 0 dy 0 f b ͑x, y, x 0 , y 0 , s͒ is the number density of the beam ions, and the constants, K b and N b , are the self-field perveance and the number of beam ions per unit axial length, respectively, defined by The nonlinear Vlasov-Maxwell equations ( 3) and ( 4) can be used to investigate detailed beam propagation and stability properties [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] over a wide range of system parameters and choices of periodic lattice functions, k x ͑s͒ and k y ͑s͒.As a general remark, it is important to note that the characteristics of the Vlasov equation (3) correspond to the single-particle equations of motion in the applied field plus self-generated fields.For example, the coefficient of ≠͞≠x is dx͞ds x 0 , the coefficient of ≠͞≠x 0 is dx 0 ͞ds 2k x ͑s͒x 2 ≠c͞≠x, etc.Moreover, the laboratory-frame Hamiltonian Ĥ for transverse singleparticle motion consistent with Eqs. ( 3) and ( 4) is given (in dimensionless units) by Ĥ͑x, y, x 0 , y 0 , s͒ For Ĥ specified by Eq. ( 6), Hamilton's equations, dx Ќ ͞ds ≠ Ĥ͞≠x 0 Ќ and dx 0 Ќ ͞ds 2≠ Ĥ͞≠x Ќ , then give the equations of motion for the transverse displacement, x Ќ ͑s͒ x͑s͒ êx 1 y͑s͒ê y , of an individual beam ion in the laboratory frame.In Sec.III, we will make use of Channell's third-order Hamiltonian averaging technique [34] to transform away the rapidly oscillating terms [35][36][37][38] in Eq. ( 6) and end up with a Hamiltonian H that depends only on slow variables ͑X, Y , X 0 , Y 0 ͒.
In subsequent sections, we will consider two classes of periodic-focusing lattices.The first corresponds to an applied alternating-gradient quadrupole magnetic field [12], with coupling coefficient defined by where B 0 q ͑s͒ ϵ ͑≠B q x ͞≠y͒ ͑0,0͒ ͑≠B q y ͞≠x͒ ͑0,0͒ .The second corresponds to a periodic-focusing solenoidal magnetic field [12,21], with coupling coefficient defined by where B 0 z ͑s͒ ϵ ͑≠B z ͞≠s͒ ͑0,0͒ .An important distinction between the two cases is evident.For a periodic quadrupole lattice, the average of k q ͑s͒ over one lattice period S is zero, R S 0 ds k q ͑s͒ 0, and the periodic solutions to Eqs. ( 3) and (4) typically correspond to elliptical cross-section beams with oscillating (as a function of s) major and minor transverse dimensions, a͑s͒ and b͑s͒.On the other hand, for a periodic-focusing solenoidal field, the average of k s ͑s͒ over one lattice period S is nonzero, R S 0 ds k s ͑s͒ S ks fi 0, and periodic solutions to Eqs. ( 3) and (4) typically correspond to circular cross-section beams with oscillating root-mean-square (rms) beam radius, r b ͑s͒.Furthermore, for the case of a solenoidal focusing field [Eq.(11)], the nonlinear Vlasov-Maxwell equations ( 3) and ( 4) are valid in a frame of reference rotating about the beam axis at the Larmor frequency [21], Following the third-order canonical transformation in Sec.III to the new Hamiltonian H in the slow variables ͑X, Y , X 0 , Y 0 ͒, in Secs.IV and V we examine the nonlinear Vlasov-Maxwell equations in the transformed variables and utilize the back-transformation to laboratory-frame variables ͑x, y, x 0 , y 0 ͒.In this regard, for specific choices of equilibrium distribution function F 0 b ͑H 0 ͒ with ≠͞≠s 0 in the transformed variables, it is important to determine key physical properties of the (periodically focused) ion beam distribution function f b ͑x, y, x 0 , y 0 , s͒ in the laboratory frame.For future reference, in laboratory-frame variables, we denote the statistical average of a phase function x͑x, y, x 0 , y 0 , s͒ by where N b R dx dy dx 0 dy 0 f b const is the number of beam ions per unit axial length.One key property of the beam distribution function f b ͑x, y, x 0 , y 0 , s͒ is the density profile defined by n b ͑x, y, s͒ Z dx 0 dy 0 f b ͑x, y, x 0 , y 0 , s͒ .
Other important properties include the rms beam radius, r b ͑s͒, the rms x and y dimensions of the beam, a͑s͒ and b͑s͒, the unnormalized total transverse beam emittance, e͑s͒, and the unnormalized x-and y-transverse beam emittances, e x ͑s͒ and e y ͑s͒.These quantities are defined by , where the statistical averages, ͗x͘, are defined according to Eq. (12).

III. CANONICAL TRANSFORMATION OF HAMILTONIAN AND PARTICLE COORDINATES TO SLOW VARIABLES
In this section, we make use of Channell's third-order Hamiltonian averaging technique [34] to transform from laboratory-frame variables ͑x, y, x 0 , y 0 ͒ to the slow variables ͑X, Y , X 0 , Y 0 ͒, with a new Hamiltonian H ͑X, Y , X 0 , Y 0 , s͒.The formalism employs a canonical transformation given by an expanded generating function [34] to transform away the rapidly oscillating terms [35][36][37][38].We formally express the laboratory-frame Hamiltonian H͑x, y, x 0 , y 0 , s͒ as H͑x, y, x 0 , y 0 , s͒ e Ĥ͑x, y, x 0 , y 0 , s͒ where Ĥ is defined in Eq. ( 6), and e is a small dimensionless parameter.In Eq. ( 15), the applied focusing potential V ͑x, y, s͒ is expressed as where U͑x, y͒ is the steady (s-independent) contribution, and Ṽ ͑x, y, s͒ is the rapidly oscillating part.For future reference, from Eqs. ( 6), (9), and (11) we express where the oscillating focusing coefficients are defined by kx ͑s͒ k x ͑s͒ 2 kx and ky ͑s͒ k y ͑s͒ 2 ky , and the 074401-4 074401-4 average focusing coefficients are defined by kx S 21 R S 0 ds k x ͑s͒ and ky S 21 R S 0 ds k y ͑s͒.For an alternating-gradient quadrupole field with R S 0 ds k q ͑s͒ 0 [Eq.( 9)], it follows that kx 2 ky 0 , kx ͑s͒ 2 ky ͑s͒ ϵ k q ͑s͒ , ( and therefore U q ͑x, y͒ 0. On the other hand, for a periodic-focusing solenoidal field [Eq.(11)] with ks S 21 R S 0 ds k s ͑s͒ fi 0, it follows that kx ky ϵ ks fi 0 , kx ͑s͒ ky ͑s͒ ϵ ks ͑s͒ k s ͑s͒ 2 ks , and therefore U sol ͑x, y͒ is generally nonzero.
We now collect together the results in Eqs. ( 29)- (33) and substitute them into the expression for the transformed Hamiltonian H ͑X, Y , X 0 , Y 0 , s͒ P e n H n ͑X, Y , X 0 , Y 0 , s͒ in Eq. (22).This readily gives where the oscillatory potentials Ṽ0 , Ṽ1 , Ṽ2 , c1 , and c2 are defined in Eqs.(30) and (31).The main objective of the present analysis is to transform to new coordinates ͑X, Y , X 0 , Y 0 ͒ such that the transformed Hamiltonian H ͑X, Y, X 0 , Y 0 , s͒ P ǹ1 e n H n ͑X, Y , X 0 , Y 0 , s͒ is slowly varying.Thus far, the generating function S͑x, y, X 0 , Y 0 , s͒ P ǹ1 e n S n ͑x, y, X 0 , Y 0 , s͒ has been arbitrary and unspecified.We now make use of this freedom to choose ͕S n ͖ in such a way that the transformed Hamiltonian H ͑X, Y, X 0 , Y 0 , s͒ is slowly varying correct to third order in the expansion parameter e.The analysis will involve s integrations over the periodic lattice functions kx ͑s 1 S͒ kx ͑s͒ and ky ͑s 1 S͒ ky ͑s͒.For future reference, it is convenient to introduce the definitions of several key quantities that occur in the subsequent analysis.We now solve Eq. ( 34) order by order, beginning with order e.

Canonical transformation to order e
Setting the coefficient of the terms of order e equal to zero in Eq. (34) gives for the first-order transformed Hamiltonian H 1 ͑X, Y , X 0 , Y 0 , s͒ where U͑X, Y͒ is the steady confining potential defined in Eq. ( 30), c͑X, Y , s͒ is the slowly varying self-field potential, and the lowest-order oscillatory confining potential Ṽ0 ͑X, Y , s͒ is defined in Eq. (30).To assure that H 1 is slowly varying, we choose the first-order generating function S 1 in Eq. ( 36) so that the final two terms on the righthand side of Eq. (36) exactly cancel.Integrating from s 0, this readily gives where a x ͑s͒ and a y ͑s͒ are defined in Eq. (35).Because ≠S 1 ͞≠s 2 Ṽ0 ͑X, Y, s͒, the expression for H 1 ͑X, Y , X 0 , Y 0 , s͒ in Eq. ( 36) reduces to which is slowly varying because of the choice of S 1 in Eq. (37).From Eqs. ( 26), (28), and (37), it also follows that the first-order transverse displacement coordinates ͑x 1 , y 1 ͒ and velocity coordinates ͑x 0 1 , y 0 1 ͒ are given by and Eqs. ( 39) and (40) lead to several simplifications in the subsequent analysis.In particular, from Eqs. ( 30), ( 32), (39), and (40), it follows that the first-and second-order contributions to the oscillatory focusing-field potential are given by and the first-and second-order contributions to the oscillatory self-field potential are given by c1 ͑X, Y , X 0 , Y 0 , s͒ 0 , c2 ͑X, Y , X 0 , Y 0 , s͒ µ In Eqs.(41) and (42), the second-order perturbed orbits x 2 ͑X, Y, X 0 , Y 0 , s͒ and y 2 ͑X, Y, X 0 , Y 0 , s͒ are yet to be determined.

Canonical transformation to order e 2
We now make use of Ṽ1 0 c1 and x 1 0 y 1 and set the coefficient of e 2 equal to zero in Eq. (34).This gives for the second-order transformed Hamiltonian where ͑x 0 1 , y 0 1 ͒ is defined in Eq. (40), and the coefficients a x ͑s͒ and a y ͑s͒ are defined in Eq. (35).We rewrite Eq. (43) in the equivalent form 074401-7 074401-7 where ͗a x ͘ ϵ S 21 R S 0 ds a x ͑s͒ and ͗a y ͘ ϵ S 21 R S 0 ds a y ͑s͒.To assure that H 2 ͑X, Y , X 0 , Y 0 , s͒ is slowly varying, we choose the second-order generating function S 2 in Eq. (44) such that where the oscillatory coefficients b x ͑s͒ and b y ͑s͒ are defined in Eq. (35).Substituting Eq. (45) into Eq.(44) then gives for the second-order transformed Hamiltonian H 2 ͑X, Y , X 0 , Y 0 , s͒ where the coefficients ͗a x ͘ and ͗a y ͘ are constants (independent of s).Furthermore, from Eqs. ( 24), ( 25), (27), and (45), the second-order transverse displacement coordinates ͑x 2 , y 2 ͒ and velocity coordinates ͑x 0 2 , y 0 2 ͒ are given by and As a general remark, from the definitions of the oscillatory coefficients, b x ͑s͒ and b y ͑s͒, in Eq. ( 35), we note that and that b x ͑s 0͒ 0 b x ͑s S͒ and b y ͑s 0͒ 0 b y ͑s S͒.

Canonical transformation to order e 3
Returning to the expression for the transformed Hamiltonian H in Eq. ( 34), we set the coefficient of e 3 equal to zero and make use of the definitions of Ṽ2 ͑X, Y, X 0 , Y 0 , s͒ and c2 ͑X, Y, X 0 , Y 0 , s͒ in Eqs. ( 41) and (42).This gives, for the third-order Hamiltonian Making use of the expressions for ͑x 0 1 , y 0 1 ͒, ͑x 2 , y 2 ͒, and ͑x 0 2 , y 0 2 ͒ in Eqs. ( 40), (47), and (48), it is straightforward to show that Eq. ( 50) can be expressed as where the s-dependent factors b x ͑s͒, b y ͑s͒, a x ͑s͒, and a y ͑s͒ are defined in Eq. ( 35) in terms of the periodic lattice functions kx ͑s 1 S͒ kx ͑s͒ and ky ͑s 1 S͒ ky ͑s͒.For the applications of interest here, not only are the averages R S 0 ds kx ͑s͒ 0 and R S 0 ds ky ͑s͒ 0, but the lattice functions kx ͑s͒ and ky ͑s͒ are assumed to have odd half-period symmetry (see examples in Fig. 1) with kx ͑s 2 S͞2͒ 2 kx ͓2͑s 2 S͞2͔͒ , ky ͑s 2 S͞2͒ 2 ky ͓2͑s 2 S͞2͔͒ .
(52) Some straightforward integration by parts that makes use of the definitions in Eq. (35) shows that the averages ͗b x ͘ and ͗b y ͘ can be expressed as and the averages ͗a x ͘ and ͗a y ͘ can be expressed as 074401-8 074401-8 It therefore follows from Eqs. ( 52) -(54) that whereas ͗a x ͘ and ͗a y ͘ are generally nonzero.We now return to the third-order Hamiltonian in Eq. ( 51) and choose the generating function S 3 to exactly cancel all rapidly oscillating terms on the right-hand side of Eq. (51).Because ͗b x ͘ 0 ͗b y ͘, we pick Here, d x ͑s͒ ϵ a 2 x ͑s͒ 2 2 kx ͑s͒b x ͑s͒ and d y ͑s͒ ϵ a 2 y ͑s͒ 2 2 ky ͑s͒b y ͑s͒, and ͗d x ͘ ϵ S 21 R S 0 ds d x ͑s͒ and ͗d y ͘ ϵ S 21  R S 0 ds d y ͑s͒.Substituting Eq. ( 56) into Eq.( 51) then gives for the slowly varying third-order Hamiltonian For future reference, we further simplify the expressions for ͗d x ͘ and ͗d y ͘.Making use of Eq. ( 35), kx ͑s͒ da x ͞ds and db x ͞ds a x ͑s͒ 2 ͗a x ͘, gives Making use of the fact that b x ͑s͒ and b y ͑s͒ vanish at s 0 and s S, Eq. ( 58) readily gives the compact representations Finally, making use of Eqs. ( 24), ( 25), (27), and (56), the third-order transverse displacement coordinates ͑x 3 , y 3 ͒ and velocity coordinates ͑x 0 3 , y 0 3 ͒ are given by and Here, use has been made of the expressions for S 1 ͑X, Y , s͒ and ͑x 2 , y 2 ͒ in Eqs. ( 37) and ( 47), and the O͑e 3 ͒ contributions to x 0 3 and y 0 3 from eS 1 ͑X, Y, s͒ eS 1 ͑x 2 e 2 x 2 , y 2 e 2 y 2 , s͒, have been included by Taylor expansion in Eq. ( 24).

B. Third-order transformed Hamiltonian and coordinate transformation
The averaging approach developed in Sec.III A represents a powerful formalism for determining the thirdorder slowly varying Hamiltonian H ͑X, Y , X 0 , Y 0 , s͒ 38), (46), and (57), we obtain, correct to third order in e, where ͗a x ͘, ͗a y ͘, ͗d x ͘, and ͗d y ͘ are defined in Eqs. ( 35) and ( 59), and we have set the expansion parameter e 1.
In Eq. ( 62), U͑X, Y ͒ is the steady focusing potential defined by U͑X, Y͒ ͑1͞2͒ ͑ kx X 2 1 ky Y 2 ͒, and c͑X, Y , s͒ is the slowly varying self-field potential in the transformed variables.It is useful to introduce the average focusing coefficients k fx and k fy defined by Rearranging terms in Eq. ( 62), and making use of Eqs. ( 59) and (63), it follows that Eq. ( 62) can be expressed in the equivalent form For completeness, it should be pointed out that if we introduce the additional canonical transformation (known as a fiber transformation) to variables ͑ X, Ỹ, X0 , Ỹ 0 ͒ defined by then the transformed Hamiltonian in Eq. ( 64) can also be expressed as For future reference in Sec.IV, we now simplify the expression for the transformed Hamiltonian H defined in Eq. ( 64) for the two cases corresponding to: (a) the alternating-gradient quadrupole focusing field in Eqs. ( 9) and ( 18) and ( b) the periodic-focusing solenoidal field in Eqs.(11) and (19).

Transformed Hamiltonian for an alternating-gradient quadrupole field
In this case, kx ͑s͒ 2 ky ͑s͒ k q ͑s͒ and kx 2 ky S 21 R S 0 ds k q ͑s͒ 0, and it follows that a x ͑s͒ 2a y ͑s͒ a q ͑s͒ ϵ Z s 64) and (67), for a periodic quadrupole lattice with k q ͑s 1 S͒ k q ͑s͒ and S 21  R S 0 ds k q ͑s͒ 0, it follows that the slowly varying Hamiltonian H q ͑X, Y, X 0 , Y 0 , s͒ is given correct to third order in e by the expression (68) Here, k fq ϵ ͑3͞S͒ R S 0 ds͓a 2 q ͑s͒ 2 ͗a q ͘ 2 ͔ is the average quadrupole focusing coefficient, and use has been For purposes of illustration, listed in Table I are the values of the lattice functions defined in Eq. ( 67) for the choice of a periodic-focusing quadrupole lattice with k q ͑s͒ kq sin͑2ps͞S͒, where kq const.Here, we have introduced the dimensionless amplitude defined by l q ϵ kq S 2 ͞2p and the lattice wave number defined by 074401-10 074401-10 k s S͞2p.An identical set of values is obtained from Eq. ( 69) for a periodic-focusing solenoidal lattice with ks ͑s͒ ks sin͑2ps͞S͒.

Transformed Hamiltonian for a periodic-focusing solenoidal field
In this case, from Eqs. ( 11) and ( 19), ks ͑s͒ k s ͑s͒ 2 ks , where ks S 21  R S 0 ds k s ͑s͒ fi 0 and S 21 R S 0 ds 3 ks ͑s͒ 0, and there is a high degree of symmetry about the beam axis because kx ͑s͒ ky ͑s͒ ϵ ks ͑s͒.In particular, from Eqs. ( 35), (59), and (63), we find From Eqs. ( 64) and ( 69), for a periodic-focusing solenoidal field with k s ͑s 1 S͒ k s ͑s͒ and S 21  R S 0 ds k s ͑s͒ ks fi 0, it follows that the slowly varying Hamiltonian H s ͑X, Y , X 0 , Y 0 , s͒ is given correct to third order in e by Here, use has been made of U s ͑X, Y ͒ ͑1͞2͒ ks ͑X 2 1 Y 2 ͒, and k fs ϵ ͑3͞S͒ R S 0 ds͓a 2 s ͑s͒ 2 ͗a s ͘ 2 ͔ is the average focusing coefficient associated with the oscillating lattice coefficient ks ͑s 1 S͒ ks ͑s͒.

Function
Value Function Value

Coordinate transformation for a periodic-focusing solenoidal field
Finally, we make use of Eqs. ( 39), (40), (47), (48), (60), and (61) to evaluate x͑X, Y , X 0 , Y 0 , s͒ X 1 ex 1 1 e 2 x 2 1 e 3 x 3 1 • • •, etc., for a periodic-focusing solenoidal field with k s ͑s 1 S͒ k s ͑s͒ and S 21  R S 0 ds k s ͑s͒ ks fi 0. Correct to order e 3 , making use of the symmetries in Eq. ( 69) and setting e 1, we readily obtain and Here, the coefficients in Eqs. ( 73) and (74) are defined in Eq. ( 69).In addition, the terms in Eqs. ( 73) and (74) proportional to a s ͑s͒, b s ͑s͒, R s 0 ds b s ͑s͒, R s 0 ds͓d s ͑s͒ 2 ͗d s ͔͘, and a s ͑s͒b s ͑s͒ are of order e, e 2 , e 3 , e 3 , and e 3 , respectively.Finally, as expected, it is evident from Eqs. ( 73) and (74) that there is a high degree of symmetry in the x and y motions for the case of a periodic-focusing solenoidal field.
In concluding this section, for the case of a periodicfocusing solenoidal field, we demonstrate an important check on Eqs. ( 73) and (74) related to the conservation 074401-12 074401-12 of canonical angular momentum [21].In laboratory-frame variables (actually Larmor-frame variables) the normalized canonical angular momentum is defined by , and making use of x 1 0 y 1 , the canonical angular momentum can be expressed as Some straightforward algebra that makes use of Eqs. ( 73)-(75) gives where XY 0 YX 0 P Q is the canonical angular momentum in the slow variables.An important conclusion is immediately evident from Eq. (76).We denote X R cosQ and Y R sinQ, where R Equation (76) then reduces to correct to order e 3 .Therefore, as expected, when c͑R, Q, s͒ is axisymmetric in the transformed variables with ≠c͞≠Q 0, it follows that P u P Q const ͑independent of s͒ , (78) corresponding to conservation of canonical angular momentum [12,21].

IV. NONLINEAR VLASOV-MAXWELL EQUATIONS IN THE SLOW VARIABLES
In this section, we examine properties of the nonlinear Vlasov-Maxwell equations for F b ͑X, Y , X 0 , Y 0 , s͒ and c͑X, Y, s͒ in the slow phase space variables (Sec.IV A) and present several examples of equilibrium solutions F 0 b ͑H 0 ͒ with ≠͞≠s 0 (Sec.IV B).The coordinate transformations in Eqs. ( 71) and (72) (periodic quadrupole field) and in Eqs. ( 73) and (74) (periodic solenoidal field) are then used in Sec.V to examine statistical averages and key properties of the periodically focused ion beam distribution f b ͑x, y, x 0 , y 0 , s͒ in laboratory-frame variables.For completeness, the linearized Vlasov-Maxwell equations in the transformed variables are presented in Sec.IV C.

A. Transformed Hamiltonian and nonlinear Vlasov-Maxwell equations in the slow variables
For present purposes, it is convenient to work with the slow variables ͑ X, Ỹ, X0 , Ỹ0 ͒ which are related to ͑X, Y , X 0 , Y 0 ͒ by the fiber transformation [38] in Eq. (65).
On the other hand, for the case of a periodic-focusing solenoidal field, ͑ X, Ỹ, X0 , Ỹ 0 ͒ and k f are defined by where use has been made of Eqs.(65), (69), and (70).
The major simplification associated with transforming to the slow variables ͑ X, Ỹ, X0 , Ỹ 0 ͒ is immediately evident from Eq. (79).In particular, the focusing coefficient k f occurring in Eq. ( 79) is both constant (independent of s) and isotropic in the transverse plane.This should be contrasted with the expression for the Hamiltonian Ĥ͑x, y, x 0 , y 0 , s͒ in the laboratory frame defined in Eq. ( 6), where the focusing coefficients k x ͑s͒ and k y ͑s͒ are rapidly oscillating functions of s.
For the Hamiltonian defined in Eq. ( 79), the singleparticle equations of motion are given by and the nonlinear Vlasov equation for F b ͑ X, Ỹ, X0 , Ỹ 0 , s͒ can be expressed as 074401-13 074401-13 where k f const is defined in Eqs. ( 80) and (81).Comparing Eqs. ( 82) and (83), note that the characteristics of the Vlasov equation ( 83) correspond to the single-particle equations of motion in the transformed variables.For example, the coefficient of ≠͞≠ X is d X͞ds X0 , the coefficient of ≠͞≠ X0 is d X0 ͞ds 2k f X 2 ≠c͞≠ X, etc.The slowly varying self-field potential c͑ X, Ỹ, s͒ occurring in Eq. ( 83) is determined self-consistently in terms of the distribution function which should be compared with Eq. ( 4).
The nonlinear Vlasov-Maxwell equations ( 83) and (84) can be used to investigate detailed equilibrium and stability properties in the slow variables ͑ X, Ỹ, X0 , Ỹ 0 ͒ over a wide range of system parameters [12,21], including beam intensity ͑K b ͒, focusing-field strength ͑k f ͒, and choices of equilibrium distribution function F 0 b ͑H 0 ͒, consistent with the assumption that the phase advance is sufficiently small [34] to assure good convergence of the averaging technique leading to Eq. ( 79).Of course, to determine properties of the (periodically focused) beam in the laboratory frame, use will be made of the back-transformation to the laboratory-frame coordinates ͑x, y, x 0 , y 0 ͒ defined in Eqs.(71) and (72) (periodic-focusing quadrupole field) or in Eqs. ( 73) and (74) (periodic-focusing solenoidal field).
For simplicity, in the subsequent analysis of Eqs.(83) and (84), we employ free boundary conditions in which the conducting wall is assumed to be infinitely far removed from the ion beam in the transverse plane.Two points are especially noteworthy in this regard.First, while the variables ͑ X, Ỹ͒ are spacelike and the variables ͑ X0 , Ỹ 0 ͒ are velocitylike in a formal analysis of Eqs. ( 83) and (84), it is clear that the backtransformation to the laboratory-frame coordinates defined in Eqs. ( 71) and (72), or in Eqs. ( 73) and (74), inexorably mixes the dependence of ͑x, y, x 0 , y 0 ͒ on the variables ͑ X, Ỹ, X0 , Ỹ 0 ͒.Second, a rigid conducting boundary in the laboratory frame will typically have a pulsating (s-dependent) shape in the transformed variables.
For example, consider the coordinate transformation for a periodic-focusing quadrupole field given in Eqs. ( 71) and (72).Correct to order e 2 , Eq. (71) gives x ͓1 2 b q ͑s͔͒ X and y ͓1 1 b q ͑s͔͒ Ỹ .Therefore, a circular cross-section conducting wall with constant radius ͑x 2 1 y 2 ͒ 1͞2 r w const in the laboratory frame corresponds to a pulsating conducting wall with elliptical cross section, X2 ͞a 2 w ͑s͒ 1 Ỹ2 ͞b 2 w ͑s͒ 1, in the transformed variables, where a 2 w ͑s͒ r 2 w ͓͞1 2 b q ͑s͔͒ 2 and b 2 w ͑s͒ r 2 w ͓͞1 1 b q ͑s͔͒ 2 .As noted earlier, the subsequent analysis in Secs.IV and V effectively assumes that the conducting wall is infinitely far removed from the beam ͑r w !`͒.
The nonlinear Vlasov-Maxwell equations ( 83) and (84) support a broad class of equilibrium solutions ͑≠͞≠s 0͒ in which the equilibrium distribution function F 0 b depends on the variables ͑ X, Ỹ, X0 , Ỹ͒ exclusively through the Hamiltonian H 0 ; i.e., Here, H 0 ͑ X, Ỹ, X0 , Ỹ 0 ͒ is defined in Eq. ( 92).Substituting Eq. (93) into Eq.( 83), it is readily shown that is an exact consequence of Eq. ( 92), where use can be made of the chain rule for differentiation to express Here, we have expressed ͑≠͞≠ X͒c 0 ͑ R͒ ͑ X͞ R͒ ͑≠͞≠ R͒c 0 ͑ R͒, etc.Because ≠c 0 ͞≠Q 0, the canonical angular momentum P Q X Ỹ0 2 Ỹ X0 is also an exact single-particle constant of the motion ͑dP Q ͞ds 0͒ in the transformed variables.Therefore, more generally speaking, the equilibrium distribution function F 0 b ͑H 0 , P Q ͒ could also depend explicitly on P Q as well as H 0 [12,21].Such beam equilibria are typically rotating and will not be considered in the present analysis.
There is clearly enormous latitude in specifying the functional form of the equilibrium distribution function F 0 b ͑H 0 ͒ in the transformed variables [12].Once the form of F 0 b ͑H 0 ͒ is specified, however, the corresponding equilibrium self-field potential c 0 ͑R͒ is to be calculated self-consistently from Eq. (84).For ≠͞≠Q 0 and ≠͞≠s 0, Eq. (84 where H 0 is defined in Eq. ( 92), and is the radial density profile in the transformed variables.
Because H 0 depends explicitly on c 0 ͑ R͒ [Eq.( 92)], the Maxwell equation ( 95) is generally a nonlinear differential equation for the self-field potential c 0 ͑ R͒.Expressing Ũ ͑1͞2͒ ͑ X02 1 Ỹ 02 ͒ and H 0 Ũ 1 k f R2 ͞2 1 c 0 ͑ R͒, and converting the velocity integration range in Eqs.(95) and (96) according to R •, the equilibrium density profile n 0 b ͑ R͒ in the transformed variables can be expressed in the equivalent form Other equilibrium properties are also readily calculated in terms of F 0 b ͑H 0 ͒.For example, in the transformed variables, because H 0 is an even function of X0 and Ỹ 0 , the average local flow velocity in the transverse plane is equal to zero; i.e., V 0 Moreover, the effective perpendicular temperature T 0 Ќb ͑ R͒ in the transformed variables is defined (in energy units) by where The general class of equilibrium distribution functions described by Eqs. ( 92) and (93) corresponds to an intense charged particle beam with circular cross section confined in the transverse plane by a uniform focusing force (k f const).This class of distribution functions F 0 b ͑H 0 ͒ has been extensively analyzed in the literature [12,16,17,21].For present purposes, we summarize here several key properties of the equilibrium and give specific examples 074401-16 074401-16 of beam equilibria F 0 b ͑H 0 ͒ in the transformed variables.These results will be very useful in Sec.V when we transform back to the laboratory frame where the beam properties are periodically focused as a function of s.

Statistical averages
The statistical average of a phase function x͑ X, Ỹ, X0 , Ỹ 0 , s͒ over the equilibrium distribution function F 0 b ͑H 0 ͒ in the transformed variables is defined in the usual manner by where H 0 is defined in Eq. ( 92), and R͒ is the number of beam particles per unit axial length.Because H 0 ͑ X, Ỹ, X0 , Ỹ0 ͒ is an even function of X, Ỹ , X0 , and Ỹ 0 , it follows that the statistical average of any odd power of X, Ỹ, X0 , or Ỹ0 , or products thereof, is equal to zero.For example, it follows that ͗ X Ỹ0 ͘ 0 0 ͗ Ỹ X0 ͘ 0 , etc. Similarly, the rms beam radius R b0 and unnormalized beam emittance e 0 in the transformed variables are defined by where R b0 and e 0 are constants (independent of s) because ≠F 0 b ͑H 0 ͒͞≠s 0. Because of the high degree of symmetry of H 0 , it also follows that Finally, some straightforward algebra that makes use of Eqs. ( 98) and (101) shows that That is, e 2 0 is directly proportional to the perpendicular pressure P 0 Ќb ͑ R͒ n 0 b ͑ R͒T 0 Ќb ͑ R͒, averaged over the radial cross section of the beam.

Radial force balance and envelope equation for the rms beam radius R b0
The formal expression for the perpendicular pressure 98) can be used to derive the equation for equilibrium radial force balance on a beam fluid element in the transformed variables.Taking the derivative of Eq. ( 98) with respect to R, it is readily shown that where we have integrated by parts with respect to Ũ and assumed ) which will be recognized as the equation for local radial force balance on a beam fluid element in the transformed variables.Solving Eq. (95) for ≠c 0 ͑ R͒͞≠ R in terms of the radial density profile n 0 b ͑ R͒, Eq. (105) can also be expressed as [21] ≠ The local force balance Eq. ( 106) can be used to derive a global radial force balance equation that relates the emittance e 0 , the focusing coefficient k f , and the rms beam radius R b0 .To briefly summarize, we operate on Eq. ( 106) with 2p R 0 d R R2 • • • and integrate by parts with respect to R, assuming P 0 Ќb ͑ R ! `͒ 0 n 0 b ͑ R ! `͒.Some straightforward algebraic manipulation that makes use of Eq. ( 103) and N b 2p R 0 d R Rn 0 b ͑ R͒ readily gives the global force balance condition [21,24] √ where 2 is the self-field perveance.Equation (107), valid for general choice of F 0 b ͑H 0 ͒, plays the role of an envelope equation for the rms beam radius R b0 and represents a powerful constraint condition on beam equilibrium properties [21].As 074401-17 074401-17 expected, if we make the identification R b0 r b0 ͞ p 2, Eq. ( 107) is similar in form to the familiar envelope equation for the outer radius r b0 of a uniform-density KV beam equilibrium [10,11] in the smooth-beam approximation ͑dr b0 ͞ds 0͒.For specified values of k f , K b , and e 2 0 , note that Eq. ( 107) can be solved for the mean-square beam radius to give As expected, we find from Eq. ( 108) that R 2 b0 increases with increasing beam intensity ͑K b ͒, increasing beam emittance ͑e 0 ͒, and decreasing focusing-field strength ͑k f ͒.

Phase advance s 0
It is convenient to introduce the effective phase advance s 0 over one lattice period S defined by s 0 e 0 R S 0 ds͞2R 2 b0 e 0 S͞2R 2 b0 , where R 2 b0 const is the mean-square beam radius defined in Eq. ( 108).This gives where s 0y ͓s 0 ͔ K b !0 ϵ p k f S is the vacuum phase advance defined in the limit of negligible beam intensity, As noted in Sec.I, the averaging technique developed in Sec.III is expected to provide good convergence properties [34] provided the phase advance s 0 is sufficiently small (s 0 , 60 ± p͞3, say).
It is important to note from Eq. (109) that s 0 ͞s 0y decreases monotonically from unity as the normalized beam intensity K b ͞2 p k f e 0 is increased.That is, selffield effects (as measured by K b ) depress the phase advance s 0 from its vacuum value s 0y .

Density inversion theorem and condition for transverse confinement
As noted earlier for the specified equilibrium distribution function F 0 b ͑H 0 ͒, when the expression for the density profile n 0 b ͑ R͒ in Eq. ( 97) is substituted into Eq.( 95), the resulting equation for the self-field potential c 0 ͑ R͒ is generally nonlinear.Without loss of generality, we take the zero of potential to be c 0 ͑ R 0͒ 0 and denote the on-axis value of beam density in the transformed variables by n 0 b ͑ R 0͒ ϵ nb .Integration of Eq. (95) from R 0 Careful examination of Eqs. ( 95) and ( 97) then shows that a necessary condition for a radially confined beam equilibrium with n 0 b ͑ R ! `͒ 0 is given by where v2 pb 4p nb Z 2 b e 2 ͞g b m b is the on-axis plasma frequency squared [12], and use has been made of the definition The inequality in Eq. ( 110) is simply a statement that the focusing force (proportional to k f b 2 b c 2 ) must exceed the repulsive spacecharge force (proportional to v2 pb ͞2) for there to be transverse confinement of the beam particles.
A further important result is evident from the expression for n 0 b ͑ R͒ in Eq. ( 97).We introduce the effective total potential V ͑ R͒ defined by V ͑ R͒ ͑1͞2͒k f R2 1 c 0 ͑ R͒.Then, taking the derivative of n 0 b ͑V ͒ with respect to V in Eq. ( 97) gives Assuming F 0 b ͓ Ũ 1 V ͑ R͔͒ Ũ!` 0 and integrating by parts with respect to Ũ in Eq. ( 111) gives for the distribution function F 0 b ͑H 0 ͒.Equation ( 112) is known as the density inversion theorem [1,21].In particular, for specified density profile n 0 b ͑ R͒, we make use of Eq. (95) to determine the self-field potential c 0 ͑ R͒ and evaluate the effective potential V ͑ R͒ ͑1͞2͒k f R2 1 c 0 ͑ R͒.Solving then for R͑V ͒, assumed to be monotonic, we evaluate ≠n 0 b ͞≠V ͑≠n 0 b ͞≠ R͒ ͑≠ R͞≠V ͒ in Eq. ( 112), which determines the equilibrium distribution function F 0 b ͑H 0 ͒ in the transformed variables.

Kinetic stability theorem
An important kinetic stability theorem [25,26] can be demonstrated from the nonlinear Vlasov-Maxwell equations ( 83) and (84) for the distribution function F b ͑ X, Ỹ, X0 , Ỹ 0 , s͒ and self-field potential c͑ X, Ỹ, s͒.In particular, we express F b F 0 b ͑H 0 ͒ 1 dF b ͑ X, Ỹ, X0 , Ỹ 0 , s͒ and c c 0 ͑ R͒ 1 dc͑ X, Ỹ, s͒ and make use of the global conservation constraints for total energy U͑s͒ and generalized entropy S͑s͒ satisfied exactly by Eqs. ( 83) and (84); i.e., 074401-18 074401-18 Here, dU͞ds 0 dS͞ds, and G͑F b ͒ is a smooth, differentiable function of F b with G͑F b !0͒ 0. Without presenting algebraic details, it can be shown [25,26] that a sufficient condition for stability is that the equilibrium distribution function F 0 b ͑H 0 ͒ be a monotonically decreasing function of energy H 0 ; i.e., That is, whenever Eq. ( 114) is satisfied, the system is stable, and the perturbations dc and dF b do not amplify.The stability theorem in Eq. ( 114) is a very powerful result and is valid nonlinearly (finite-amplitude perturbations) as well as for small-amplitude perturbations.For example, Eq. ( 114) implies that a beam with thermal equilibrium [1,18,21] distribution F 0 b ͑H 0 ͒ [Eq.( 115)] is stable and can propagate quiescently over large distances.On the other hand, a Kapchinskij-Vladimirskij beam equilibrium [Eq.( 116)] has an inverted population in H 0 , and there is (in principle) free energy available to cause the perturbations dc and dF b to amplify [10][11][12][13][14].

Examples of self-consistent beam equilibria
For future reference, we briefly consider several examples of beam equilibria, F 0 b ͑H 0 ͒, in the transformed variables ͑ X, Ỹ, X0 , Ỹ 0 ͒.Specifically, we consider the following choices of F 0 b ͑H 0 ͒: Thermal Equilibrium: [1,18,21] Here, nb and TЌb are positive constants with dimensions of density and temperature (energy units), respectively, Hamiltonian defined in Eq. ( 92), and U͑x͒ is the unit step function defined by U͑x͒ 1 for 0 # x , 1 and U͑x͒ 0 for x .1.We take the onaxis self-field potential to be c 0 ͑ R 0͒ 0 and identify nb n 0 b ͑ R 0͒ with the on-axis beam density.For each choice of F 0 b ͑H 0 ͒ in Eqs. ( 115)-( 117), the self-field potential c 0 ͑ R͒ is determined self-consistently in terms of the beam density 95).Finally, for the general class of beam equilibria, F 0 b ͑H 0 ͒, the transverse temperature profile T 0 Ќb ͑ R͒ is defined by Eq. (98).
A detailed evaluation of beam equilibrium properties for the choices of distribution functions in Eqs. ( 115)-( 117) is presented elsewhere [21], and essential results are summarized in Table II.For each example, the inequality k f b 2 b c 2 .v2 pb ͞2g 2 b is required to assure radial confinement of the beam particles, where v2 pb 4pZ 2 b e 2 nb ͞g b m b is the on-axis plasma frequency squared.Moreover, in each case, the rms beam radius R b0 and unnormalized beam emittance e 0 defined in Eq. ( 101) are related by the global force balance constraint in Eq. ( 107), and the unnormalized emittance e 0 can be expressed as the average over perpendicular pressure given in Eq. ( 103).
It is evident from Table II that the equilibrium profiles for the density n 0 b ͑ R͒ and perpendicular temperature T 0 Ќb ͑ R͒ differ significantly for the three choices of equilibrium distribution functions in Eqs. ( 115)-(117).First, for the choice of thermal equilibrium distribution in Eq. ( 115), we note from Table II that the equilibrium density profile exhibits a highly nonlinear dependence on the self-field potential c 0 ͑ R͒, which must generally be determined by numerical integration of Eq. ( 95).The corresponding density profile, 2c 0 ͑ R͔͖͒, is generally bell shaped and radially diffuse, assuming a maximum value ͑ nb ͒ at R 0, and decreasing monotonically to zero   II and Eq. ( 95) is found to be radially very broad [21] in units of the thermal Debye length; i.e., R b0 ¿ l D ϵ ͑g 2 b TЌb ͞4p nb Z 2 b e 2 ͒ 1͞2 .In this case, the density is approximately constant in the beam interior with n 0 b ͑ R͒ Ӎ nb const, and n 0 b ͑ R͒ drops rapidly to exponentially small values over a few Debye lengths at the beam surface.For the choice of equilibrium distribution function in Eq. ( 115), it also follows that the transverse temperature profile is uniform over the beam cross section, with T Ќb ͑ R͒ TЌb const.Moreover, because ≠F 0 b ͑H 0 ͒͞≠H 0 # 0 for the choice of distribution function in Eq. ( 115), it follows from Eq. ( 114) that the equilibrium is stable [25,26].
By contrast, at any beam intensity, the choice of (monoenergetic) equilibrium distribution function in Eq. ( 116) gives a step-function density profile, with n 0 b ͑ R͒ nb const for 0 # R , r b0 p 2 R b0 , and n 0 b ͑ R͒ 0 for R .r b0 .In this case, the beam has a "sharp" outer boundary at radius r b0 determined self-consistently from Moreover, from Table II, unlike the thermal equilibrium case, the transverse temperature profile is parabolic, with T 0 Ќb ͑ R͒ TЌb ͑1 2 R2 ͞r 2 b0 ͒ in the beam interior ͑0 # R , r b0 ͒.Finally, a most important feature of Eq. ( 116) is that F 0 b ͑H 0 ͒ has a highly inverted population in energy, which is (singularly) peaked at H 0 TЌb ͞g b m b b 2 b c 2 .Therefore, as expected, there is free energy available to drive collective instabilities [10][11][12][13][14] for the choice of equilibrium distribution function in Eq. ( 116), at least at sufficiently high beam intensity.Finally, from Table II, the choice of waterbag equilibrium distribution [16,17,21] in Eq. ( 117) also gives a density profile n 0 b ͑ R͒ with sharp outer boundary at radius r b0 .In this case, r b0 is determined self-consistently from where I 0 ͑x͒ is the modified Bessel function of the first kind of order zero, and l D ͑g 2 b TЌb ͞4p nb Z 2 b e 2 ͒ 1͞2 is the thermal Debye length.Unlike the KV beam equilibrium, however, we note from Table II that the density profile n 0 b ͑ R͒ decreases monotonically from the value nb at R 0 to zero at R r b0 .Moreover, the transverse temperature profile has exactly the same radial shape as the density profile, with T 0 Ќb ͑ R͒ ͑ TЌb ͞ nb ͒n 0 b ͑ R͒.Similar to the case of a thermal equilibrium beam, the distribution function in Eq. ( 117) satisfies ≠F 0 b ͑H 0 ͒͞≠H 0 # 0, and the waterbag equilibrium is expected to be stable [25,26] by virtue of Eq. (114).

V. PERIODICALLY FOCUSED BEAM PROPERTIES IN THE LABORATORY FRAME
As discussed in Sec.IV, in the slow variables ͑ X, Ỹ, X0 , Ỹ 0 ͒, a wide variety of equilibrium and stability properties can be calculated in a straightforward manner because the focusing force is constant (k f const) in the transformed variables, and the simplest class of equilibrium distribution functions, F 0 b ͑H 0 ͒, correspond to beams with circular cross section.When transformed back to the laboratory frame, however, the beam is periodically focused, and its properties are generally s dependent.In this section, we carry out the back-transformation to the laboratory frame and calculate several properties of the beam, such as (a) the distribution function f b ͑x, y, x 0 , y 0 , s͒ (Sec.V A), (b) statistical averages such as the mean-square transverse beam dimensions, ͗x 2 ͘ ͑s͒ and ͗ y 2 ͘ ͑s͒, and the unnormalized emittances, e x ͑s͒ and e y ͑s͒ (Sec.V B), and (c) macroscopic properties of the beam in the laboratory frame, such as the density profile n b ͑x, y, s͒ (Sec.V C).Throughout Sec.V, extensive use will be made of the inverse coordinate transformation, X͑x, y, x 0 , y 0 , s͒, Ỹ͑x, y, x 0 , y 0 , s͒, etc., defined correct to order e 3 in the Appendix [Eqs.(A1) and (A2) for a periodic-focusing quadrupole field and Eqs.(A3) and (A4) for a periodic-focusing solenoidal field].In calculating the statistical averages in Sec.V B, we will also make use of the forward coordinate transformation, x͑ X, Ỹ, X0 , Ỹ0 , s͒, y͑ X, Ỹ, X0 , Ỹ 0 , s͒, etc., defined correct to order e 3 in Eqs.(87)-(90).In this regard, it is important to keep in mind the relative ordering of the various terms in Eqs. ( 87) -(90) and Eqs.(A1)-(A4).For example, a q ͑s͒ and ͗a q ͘ are of order e, b q ͑s͒ is of order e 2 , etc. [see Eq. ( 91)].
Therefore, Eqs. ( 121) and (122) together with the coordinate transformations in Eqs.(A1)-(A4), map the equilibrium distribution function F 0 b ͑H 0 ͒, which is uniformly focused and has circular cross section in the transformed variables, into a pulsating, periodically focused distribution function in the laboratory frame.
The result in Eq. ( 121), together with the associated definitions in Eq. ( 122) and the Appendix, make accessible for the first time a broad class of high-intensity, periodically focused distribution functions that are analytically tractable, in addition to the familiar Kapchinskij-Vladimirskij equilibrium.Therefore, it is anticipated that Eq. ( 121) together with the results in Sec.IV B will be very useful in providing input data for numerical simulation studies based on the nonlinear Vlasov-Maxwell equations, as well as experimental studies of beam matching into periodic-focusing channels.

Periodic focusing quadrupole field
To illustrate the application of Eq. ( 125) to a periodicfocusing quadrupole field, we make use of Eqs. ( 87) and (125) to evaluate ͗x 2 ͘ ͑s͒.This readily gives 100)] and because ͓ R s 0 ds b q ͑s͔͒ 2 is of order e 6 [Eq.( 91)], we obtain correct to order e 3 .Here, We conclude from Eqs. ( 126) and (127) that the beam cross section in the laboratory frame corresponds to a pulsating ellipse (in an rms sense) with minor axes, a͑s͒ and b͑s͒.It should also be kept in mind that Eqs. ( 126) and (127) apply to the entire class of equilibrium distributions, F 0 b ͑H 0 ͒, in the transformed variables, and that self-field effects are allowed to be arbitrarily intense, consistent with radial confinement of the beam particles by the focusing field [Eq.( 110)].
It is important to recognize the implications and limitations of Eq. (129); that is, the transverse emittances, e x ͑s͒ and e y ͑s͒, are conserved quantities when the back-transformation to the laboratory frame is carried out.First, and very important, Eq. ( 129 121)-( 123)], at least to order e 3 .That is, Eq. ( 129) pertains to the transverse emittances associated with periodically focused beam equilibria in the laboratory frame.Of course, if the system is perturbed about equilibrium, i.e., F b ͑ X, Ỹ, X0 , Ỹ 0 , s͒ F 0 b ͑H 0 ͒ 1 dF b ͑ X, Ỹ , X0 , Ỹ 0 , s͒, then the perturbations evolve self-consistently according to Eqs. ( 83) and (84), and there will be a corresponding change in the total laboratory-frame transverse emittances e x ͑s͒ and e y ͑s͒ [see Eqs. ( 14) and (120)] associated with the changes in F b and c.It is well known that such variations in the laboratory-frame emittances about equilibrium values can be sizable [2,5], particularly if the equilibrium distribution F 0 b ͑H 0 ͒ is unstable, and there is a significant redistribution of particles in phase space.Second, and also important, in the equilibrium case it is important to keep in mind that Eq. ( 129) is an approximate result obtained in the context of the asymptotic analysis in Secs.III and IV (correct to order e 3 ), and there are undoubtedly corrections to Eq. ( 129) of order e 4 or smaller.Finally, referring ahead to Sec.V C, there is an important comparison to be made with Sacherer's classic analysis [39] of periodically focused intense beam propagation.For the case of constant emittances, e x and e y , analysis of Sacherer's rms envelope equations [39] shows that beam density profiles with an oscillatory elliptical cross section constitute self-consistent periodically focused solutions for intense beam propagation through a periodic quadrupole lattice.In Sec.V C, for a periodic quadrupole lattice, we find that the laboratory-frame density profile n b ͑x, y, s͒ corresponding to the beam equilibrium F 0 b ͑H 0 ͒ indeed has an oscillatory elliptical cross section [Eq.( 136)].In addition, however, the asymptotic analysis presented here demonstrates that the transverse emittances are constant, at least to order e 3 .

Periodic focusing solenoidal field
We now summarize several key results for statistical averages in the laboratory frame for the case of a periodicfocusing solenoidal field.In this case, we make use of Eq. ( 125), together with the coordinate transformations for x͑ X, Ỹ, X0 , Ỹ 0 , s͒, y͑ X, Ỹ, X0 , Ỹ0 , s͒, etc., defined in Eqs.(89) and (90).The high degree of symmetry in the x and y motions in the solenoidal focusing field is, of course, reflected in the statistical averages.Paralleling the analysis for the quadrupole focusing case, and making use of Eqs.(89), (125), and ͗ X X0 ͘ 0 0 ͗ Ỹ Ỹ0 ͘ 0 , we readily obtain correct to order e 3 .Here, ͗ X2 ͘ 0 ͗ Ỹ2 ͘ 0 R 2 b0 ͞2 follows from Eq. (102).From Eq. ( 130), the mean-square radius, r 2 b ͑s͒ ͗x 2 1 y 2 ͘ ͑s͒, in the laboratory frame is given by Because b s ͑s͒ is of order e 2 [Eq.( 91)], we can also express ͓1 2 b s ͑s͔͒ 2 Ӎ 1 2 2b s ͑s͒ correct to order e 3 in Eq. ( 131).From Eqs. ( 130) and ( 131), for a periodicfocusing solenoidal field, we conclude that the beam cross section is circular and that the rms beam radius r b ͑s͒ oscillates with periodicity length S.
Other statistical averages of practical interest are also readily calculated from Eqs. ( 89), (90), and (125).Without presenting details, some straightforward algebraic manipulation that makes use of Eq. ( 102) and ͗xx 0 1 yy 0 ͘ 2 ͑s͒ ͓a s ͑s͒ 2 ͗a s ͔͘ 2 ͗ X2 1 Ỹ2 ͘ 2 0 , correct to order e 3 .Equations ( 131) and ( 132) can be used to calculate the total transverse emittance in the laboratory frame, defined by e 2 ͑s͒ 4͓͗x 2 1 y 2 ͘ ͗x 02 1 y 02 ͘ 2 ͗xx 0 1 yy 0 ͘ 2 ͔.Because b s ͑s͒ is of order e 2 , we note that 1 2 b 2 s ͑s͒ Ӎ 1 correct to order e 3 .Therefore, from Eqs. ( 131) and (132) we obtain correct to order e 3 .Therefore, from Eq. ( 133), the transverse emittance e͑s͒ is a conserved quantity (independent of s) correct to order e 3 .Similar to the quadrupole focusing case, the results summarized in Eqs.(130)-( 133) have a wide range of applicability.In particular, the results apply to the entire class of equilibrium distribution functions F 0 b ͑H 0 ͒, and the self fields are allowed to have arbitrary intensity.

C. Macroscopic profiles in the laboratory frame
Equation (121), supplemented by the definition of H 0 in Eq. ( 122), provides a closed expression for the laboratory-frame distribution function f b ͑x, y, x 0 , y 0 , s͒ in terms of the transformed distribution function F b ͑ X, Ỹ, X0 , Ỹ0 , s͒ for the case where F b F 0 b ͑H 0 ͒ corresponds to an equilibrium distribution in the transformed variables.As such, Eqs. ( 121) and (122), supplemented by the orbit transformations for X͑x, y, x 0 , y 0 , s͒, Ỹ͑x, y, x 0 , y 0 , s͒, etc., defined in Eqs.(A1)-(A4) can be used for a direct evaluation of a wide variety of macroscopic profiles in the 074401-23 074401-23 laboratory frame, such as the beam density profile, n b ͑x, y, s͒ R dx 0 dy 0 f b ͑x, y, x 0 , y 0 , s͒, the transverse flow velocity, V b ͑x, y, s͒ n 21  b b b c R dx 0 dy 0 ͑x 0 êx 1 y 0 êy ͒f b ͑x, y, x 0 , y 0 , s͒, etc.
As expected, the density contours in the laboratory frame have an elliptical cross section for a periodic-focusing quadrupole field [Eq.( 136)] and a circular cross section for a periodic-focusing solenoidal field [Eq.( 137)].
Moreover, in both cases, it can be shown from Eqs. ( 136) and (137) that the number of particles per unit axial length is conserved 136) and (137) are particularly attractive representations of the density profile n b ͑x, y, s͒ because the s-dependent distortion of the profiles appears explicitly in the definitions of R͑x, y, s͒.It should be pointed out, however, for continuously varying profiles n 0 b ͑ R͒, that Eqs. ( 136) and (137) can also be Taylor expanded locally about a radius r ͑x 2 1 y 2 ͒ 1͞2 in the laboratory frame, treating b q ͑s͒ and b s ͑s͒ as small parameters (of order e 2 ).For example, in the periodic quadrupole case, R͑x, y, s͒ ͓x 2 ͑͞1 2 b q ͒ 2 1 y 2 ͑͞1 1 b q ͒ 2 ͔ 1͞2 Ӎ r 1 b q ͑s͒ ͑x 2 2 y 2 ͒͞r for small b q ͑s͒.The expression for n b ͑x, y, s͒ in Eq. (136) can then be approximated by The self-field potential c͑x, y, s͒ in the laboratory frame is determined self-consistently in terms of the density profile n b ͑x, y, s͒ by integrating the Maxwell equation (4).For the case of a periodic-focusing quadrupole field, the density profile has the form given in Eq. (136).
We introduce the minor dimensions defined by a 2 ͑s͒ ͓1 2 b q ͑s͔͒ 2 R 2 b0 and b the mean-square radius associated with the equilibrium distribution F 0 b ͑H 0 ͒ in the transformed variables defined in Eq. (101).Substituting Eq. (136) into Eq.( 4) then gives where R͑x, y, s͒͞R b0 is defined by In Eq. ( 140), we have introduced the scaled radial variable R͑x, y, s͒͞R b0 in the argument of n 0 b ͑ R͞R b0 ͒ without loss of generality.Taking c, ≠c͞≠x, and ≠c͞≠y to be equal to zero at ͑x, y͒ ͑0, 0͒, the exact solutions to Eq. ( 140) for the self-field force components in the laboratory frame, F s x 2≠c͞≠x and F s y 2≠c͞≠y, are given by [39] where T ͑x, y, s, j͒ is defined by For a specified functional form of n 0 b ͑ R͞R b0 ͒, Eq. (142) can be used to calculate the detailed dependence of 2≠c͞≠x and 2≠c͞≠y on ͑x, y, s͒.
The transverse self-field force is even simpler to determine for the case of a periodic-focusing solenoidal field because of the azimuthal symmetry of c͑r, s͒ in the laboratory frame.In this case we introduce the mean-square radius r 2 b ͑s͒ ͓1 2 b s ͑s͔͒ 2 R 2 b0 defined in Eq. ( 131) and express R͑x, y, s͒͞R b0 ͓͑x  145) can be used to determine the corresponding self-field force in the laboratory frame which is a periodic function of s.

Transverse flow velocity V b ͑ ͑ ͑x, y, s͒ ͒ ͒
For a specified equilibrium distribution function, F 0 b ͑H 0 ͒, in the transformed variables, the corresponding laboratory-frame distribution function, f b ͑x, y, x 0 , y 0 , s͒, defined in Eq. ( 121) can be used to calculate other macroscopic properties, such as the transverse temperature profiles, flow velocity components, etc.We illustrate this with a direct calculation of the transverse flow velocity defined (in dimensional units) by For present purposes, we consider the case of a periodicfocusing solenoidal field, where the inverse coordinate transformation, X͑x, y, x 0 , y 0 , s͒, Ỹ ͑x, y, x 0 , y 0 , s͒, etc., occurring in the definition of H 0 in Eq. ( 122) is defined in Eqs.(A3) and (A4).Of particular importance when calculating the velocity moments in Eq. ( 146) are the symmetries associated with the kinetic energy term, ͑1͞2͒ ͓ X02 ͑x, y, x 0 , y 0 , s͒ 1 Ỹ 02 ͑x, y, x 0 , y 0 , s͔͒, in Eq. (122 ͓ X02 ͑x, y, x 0 , y 0 , s͒ 1 Ỹ 02 ͑x, y, x 0 , y 0 , s͔͒ To calculate the transverse flow velocity defined in Eq. ( 146), we make use of Eqs. ( 121), (122), and (147), keeping in mind the relative size of the various terms in Eq. ( 147) [see Eq. ( 91)] and approximating (for example) shows that the transverse flow velocity is purely radial with where V rb ͑r, s͒ is defined by correct to order e 3 .Here, x êx 1 y êy r êr , where êr cosu êx 1 sinu êy is a unit vector in the radial direction, and we have neglected terms proportional to b s ͑s͒ ͓ R s 0 ds b s ͑s͔͒, which are of order e 5 .The transverse flow velocity can be calculated in a similar manner for a periodic-focusing quadrupole field, although the flow pattern is more complicated than in Eq. (148) because of the elliptical cross section of the beam.

D. Range of validity of asymptotic expansion procedure
To conclude Sec.V, we summarize the illustrative conditions required for validity of the present asymptotic expansion procedure for the case of a periodic-focusing quadrupole lattice with the sinusoidal waveform considered in Table I.Here, the strength of the focusing field is measured by the dimensionless parameter l q ϵ kq S 2 ͞2p, and the characteristic vacuum phase advance is defined by s oy p k fq S, where k fq ϵ ͑3͞2͒l 2 q ͞S 2 [see Table I and Eq. ( 109)].For this choice of focusing lattice, it follows that l q and s oy are related by Therefore, l q , ͑2͞3͒ 1͞2 0.82 corresponds to s oy , 1 (vacuum phase advance less than 60 ± ).We now consider the relative size of the various terms in the coordinate transformation relating ͑x, y, x 0 , y 0 ͒ to ͑ X, Ỹ, X0 , Ỹ0 ͒ in Eqs.(87) and (88).Here, we estimate the characteristic (maximum) values by jxj ϳ j Xj ϳ jyj ϳ j Ỹj ϳ r b , jx 0 j ϳ j X0 j ϳ jy 0 j ϳ j Ỹ 0 j ϳ r b ͞S , (151) jcj ϳ K b , where r b is the characteristic (rms) beam radius, K b is the dimensionless self-field perveance defined in Eq. ( 5), and use has been made of Maxwell's equation (84) and N b ϳ pr 2 b n b to estimate jcj ϳ K b .The correction terms to x X and y Ỹ in Eq. (87) then stand in the ratio jb q j: j R s 0 ds b q j S .
Similarly, the correction terms to x 0 X0 and y 0 Ỹ 0 in Eq. (88) stand in the ratio jb q j:Sja q j:Sja q b q j:S É Z s 0 ds͓d q 2 ͗d q ͔͘ We now make use of the entries in Table I and the estimate jcj ϳ K b in Eq. ( 151) to compare the size of the various terms in Eqs. ( 152) and (153).We readily find that the correction terms in Eq. ( 152) stand in the ratio l q 2p : 2l q ͑2p͒ 2 . (154) Similarly, the correction terms in Eq. ( 153) stand in the ratio l q 2p :l q : l 2 q 2p : l 2 q 2p : l q ͑2p͒ 2 ?
From Eq. ( 107), the largest value of the self-field perveance K b allowed in the limit of negligibly small transverse emittance is K b ϳ k fq r 2 b , which gives K b S 2 ͞r 2 b ϳ k fq S 2 ͑3͞2͒l 2 q .Therefore, in the limit of intense spacecharge field, the ratio of terms in Eq. (155) reduces to l q 2p :l q : l 2 q 2p : l 2 q 2p : 3 2 Therefore, from Eqs. ( 154) and (156), we conclude that the key small parameter required for validity of the present asymptotic analysis is e l q ͞2p , 1, at least for the case of a sinusoidal quadrupole focusing field considered in Table I.From Eq. (150), this corresponds to s oy ͞2p , ͑3͞2͒ 1͞2 , which leads to the conjecture (Sec.I) that the phase advance s 0 should be smaller than 60 ± ( p͞3).The important practical test of the range of validity awaits detailed comparison with experiment and numerical simulations.

VI. CONCLUSIONS
In this paper, we have developed and applied a thirdorder Hamiltonian averaging technique for investigating solutions to the nonlinear Vlasov-Maxwell equations for the case of an intense ion beam propagating through a periodic-focusing quadrupole field or a periodic-focusing solenoidal field.The formalism used a canonical transformation given by an expanded generating function to transform away the rapidly oscillating terms and end up with a Hamiltonian H that depends only on slow variables.The assumptions and theoretical model were summarized in Sec.II, including the nonlinear Vlasov-Maxwell equations for the distribution function f b ͑x, y, x 0 , y 0 , s͒ and self-field potential c͑x, y, s͒ in the laboratory frame.In Sec.III, we made use of Channell's third-order Hamiltonian averaging technique [34] to transform from laboratory-frame variables ͑x, y, x 0 , y 0 ͒ to a new Hamiltonian H ͑ X, Ỹ, X0 , Ỹ 0 , s͒ in the slow variables ͑ X, Ỹ, X0 , Ỹ0 ͒ correct to order e 3 .The formalism employed a canonical transformation given by an expanded generating function to transform away the rapidly oscillating terms.This led to a Hamiltonian H ͑ X, Ỹ, X0 , Ỹ 0 , s͒ in the transformed variables of the form given in Eq. ( 79), where k f const.An im-portant by-product of the generating function analysis was the determination of the coordinate transformation that relates the laboratory-frame variables ͑x, y, x 0 , y 0 ͒ to the slow variables ͑ X, Ỹ, X0 , Ỹ 0 ͒ [Eqs.(87)-(90)].The major simplification associated with transforming to the slow variables ͑ X, Ỹ, X0 , Ỹ 0 ͒ is immediately evident from the expression for H ͑ X, Ỹ, X0 , Ỹ0 , s͒ in Eq. (79).In particular, the focusing coefficient k f is both constant (independent of s) and isotropic in the transverse plane.This should be contrasted with the expression in Eq. ( 6) for the Hamiltonian Ĥ͑x, y, x 0 , y 0 , s͒ in the laboratory frame, where the focusing coefficients k x ͑s͒ and k y ͑s͒ are rapidly oscillating functions of s.In Sec.IV, following a discussion of the nonlinear Vlasov-Maxwell equations for F b ͑ X, Ỹ, X0 , Ỹ 0 , s͒ and c͑ X, Ỹ, s͒ in the transformed variables, we presented several examples of axisymmetric equilibrium solutions, i.e., distribution functions F 0 b ͑H 0 ͒ with ≠͞≠s 0 and ≠͞≠Q 0 corresponding to constant-radius beam equilibria with a circular cross section in the transformed variables [12,21].Of particular note is the class of distribution functions that satisfy ≠F 0 b ͑H 0 ͒͞≠H 0 # 0, which can be shown to be stable [25,26].Finally, in Sec.V, we exploited the inverse coordinate transformation, X͑x, y, x 0 , y 0 , s͒, Ỹ ͑x, y, x 0 , y 0 , s͒, etc., to determine properties of the periodically focused distribution function f b ͑x, y, x 0 , y 0 , s͒ in the laboratory frame, correct to order e 3 , consistent with the class of constant-radius circular cross-section beam equilibria F 0 b ͑H 0 ͒ in the transformed variables.A wide range of important physical quantities were determined, including the distribution function f b ͑x, y, x 0 , y 0 , s͒; statistical averages such as the transverse mean-square beam dimensions, ͗x 2 ͘ ͑s͒ and ͗ y 2 ͘ ͑s͒, and the unnormalized emittances, e x ͑s͒ and e y ͑s͒; and macroscopic properties such as the number density of beam particles, n b ͑x, y, s͒ R dx 0 dy 0 f b ͑x, y, x 0 , y 0 , s͒, the self-field potential, c͑x, y, s͒, etc.Finally, in Sec.V D, we summarized the illustrative conditions required for validity of the present asymptotic expansion procedure for the case of a periodic-focusing quadrupole lattice with sinusoidal waveform (Table I).The important practical test awaits detailed comparison with experiment and numerical simulations.

ACKNOWLEDGMENTS
This research was supported by the Department of Energy and the Short Pulse Spallation Source Project and LANSCE Division of Los Alamos National Laboratory, and in part by the Office of Naval Research.

FIG. 1 .
FIG. 1. Examples of periodic-focusing lattice functions with odd half-period symmetry [Eq.(52)] corresponding to (a) a sinusoidal lattice function, kx ͑s͒ k sin͑2ps͞S͒, with k const, and ( b) a periodic step-function lattice with amplitude k const and filling factor h.

b c 2 2 v2 pb ͞2g 2 b ͒͞ k f b 2 b c 2 d ø 1 ,
b0 in Eq. (116) for 0 # R , r b0 ϵ p 2R b0 ; for 0 # R , r b0 ϵ p 2R b0 ; (zero, otherwise) (zero, otherwise) 3. Waterbag distribution I0͑rb0͞lD ͒2I0͑ R͞lD ͒ I0͑rb0͞lD ͒21 TЌb n 0 b ͑ R͒ nb Determined in Eq. (117) for 0 # R , r b0 ; from Eq. (103) (zero, otherwise) DAVIDSON, HONG QIN, AND PAUL J. CHANNELL 074401 (1999) with n 0 b ͑ R ! `͒ 0. At sufficiently low beam intensity with v2 pb ͞2g 2 b ø k f b 2 b c 2 , it is found that the density profile is approximately Gaussian, with n 0 b ͑ R͒ Ӎ nb exp͕2g b m b b 2 b c 2 k f R2 ͞2 TЌb ͖.On the other hand, at very high beam intensity with ͑k f b 2 the density profile evaluated numerically from Table ) pertains to the class of equilibrium distributions F 0 b ͑H 0 ͒, , constant-radius beam equilibria in the transformed variables (Sec.IV B) and transform back to matched, periodically focused solutions with period S in the laboratory frame [Eqs.(

correct to order e 3 .͓r 2 2 .
Equation (138)  shows quite naturally the quadrupole distortion of the density profile in the laboratory frame by the periodic-focusing quadrupole field.Similarly, in the periodic solenoidal case, R͑x, y, s͓͒͑x 2 1 y 2 ͒͑͞1 2 b s ͒ 2 ͔ 1͞2 Ӎr 1 b s ͑s͒r, and the expression for n b ͑x, y, s͒ in Eq. (137) can be approximated by n b ͑x, y, s͒ n 0 b ͑r͒ 1 b s ͑s͒ n 0 b ͑r͔͒ ,(139)correct to order e 3 .The (pulsating) profile in Eq. (139) of course remains axisymmetric in the laboratory frame for the case of a periodic-focusing solenoidal field.Self-field potential c͑ ͑ ͑x, y, s͒ ͒ ͒

Z
r͞r b ͑s͒.Substituting Eq. (137) into the Maxwell equation (4) for c͑r, s͒ then gives r͞r b ͑s͒ 0 dT Tn 0 b ͑T ͒ .(145) Similar to the periodic quadrupole case, once the functional form of n 0 b ͑ R͒ is determined self-consistently for a specified equilibrium distribution function, F 0 b ͑H 0 ͒, in the transformed variables, Eq. ( where we have estimated j≠c͞≠ Xj ϳ jcj͞r b , etc.

TABLE II .
Equilibrium properties for various choices of F 0 b ͑H 0 ͒.