One-dimensional kinetic description of nonlinear traveling-pulse (soliton) and traveling-wave disturbances in long coasting charged particle beams

This paper makes use of a one-dimensional kinetic model to investigate the nonlinear longitudinal dynamics of a long coasting beam propagating through a perfectly conducting circular pipe with radius $r_{w}$. The average axial electric field is expressed as $\langle E_{z}\rangle=-(\partial/\partial z)\langle\phi\rangle=-e_{b}g_{0}\partial\lambda_{b}/\partial z-e_{b}g_{2}r_{w}^{2}\partial^{3}\lambda_{b}/\partial z^{3}$, where $g_{0}$ and $g_{2}$ are constant geometric factors, $\lambda_{b}(z,t)=\int dp_{z}F_{b}(z,p_{z},t)$ is the line density of beam particles, and $F_{b}(z,p_{z},t)$ satisfies the 1D Vlasov equation. Detailed nonlinear properties of traveling-wave and traveling-pulse (solitons) solutions with time-stationary waveform are examined for a wide range of system parameters extending from moderate-amplitudes to large-amplitude modulations of the beam charge density. Two classes of solutions for the beam distribution function are considered, corresponding to: (a) the nonlinear waterbag distribution, where $F_{b}=const.$ in a bounded region of $p_{z}$-space; and (b) nonlinear Bernstein-Green-Kruskal (BGK)-like solutions, allowing for both trapped and untrapped particle distributions to interact with the self-generated electric field $\langle E_{z}\rangle$. .


I. INTRODUCTION
High-energy accelerators and transport systems [1][2][3][4][5][6] have a wide variety of applications ranging from basic research in high energy and nuclear physics, to applications such as spallation neutron sources, medical physics, and heavy ion fusion. As a consequence, it is increasingly important to develop an improved understanding of collective processes and the nonlinear dynamics of intense charged particle beam systems. While there has been considerable progress in three-dimensional numerical and analytical investigations of the nonlinear Vlasov-Maxwell equations describing intense beam propagation, there is also considerable interest in the development and application of simplified one-dimensional kinetic models to describe the longitudinal nonlinear dynamics of long coasting beams [7][8][9][10][11][12][13][14][15] in linear (linac) or large-major-radius ring geometries. The present paper employs the one-dimensional kinetic formalism recently developed by Davidson and Startsev [14] for a long coasting beam propagating through a perfecting conducting circular pipe with radius r w . In Ref. [14] the average longitudinal electric field is expressed as E z (z, t) = 3 , where e b is the particle charge, g 0 and g z are constant geometric factors that depend on the location of the conducting wall and the shape of the transverse density profile, and λ b (z, t) = dp z F b (z, p z , t) is the line density.
In a previous application of the 1D kinetic formalism developed in Ref. [14], the analyses in Ref. [15] assumed that the longitudinal distribution F b (z, p z , t) corresponded to a so-called waterbag distribution [16][17][18][19], where F b = const. within moving boundaries in the phase space (z, p z ). The weakly nonlinear analysis in Ref. [15] showed that disturbances moving near the sound speed evolve according to the Korteweg-deVries (KdV) equation [20][21][22][23][24]. The classical KdV equation, which arises in several areas of nonlinear physics in which there are cubic dispersive corrections to sound-wave-like signal propagation, also has the appealing feature that it's exactly solvable using inverse scattering techniques.
While the analysis in Ref. [15] reveals many interesting properties of the nonlinear evolution of longitudinal disturbances in intense charged particle beams, it is restricted to the weakly nonlinear regime. In the present analysis, we remove the restriction to the weakly nonlinear regime, and make use of the 1D kinetic model developed in Ref. [14], allowing for moderate to large-amplitude modulation in the charge density of the beam particles. The organization of this paper is the following. In Sec. II, the 1D kinetic model [14] is briefly reviewed (Sec. II A), and exact (local and nonlocal) nonlinear conservation constraints are derived (Sec. II B) for the conservation of particle number, momentum, and energy per unit length of the beam, making use of the nonlinear Vlasov equation for F b (z, p z , t) in Eq. (1), and the expression for E (z, t) in Eq. (2). Removing the assumption of weak nonlinearity made in Ref. [15], Sec. III focuses on use of the fully nonlinear kinetic waterbag model we examine the kinetic model based on Eqs. (9) and (10) [equivalent to Eqs. (1) and (2)] for an even broader class of distribution functions F b (z, p z , t), recognizing that Eqs. (9) and  [25,26], except for the fact that Eq. (10), which connects the effective potential φ (z, t) to the line density λ b (z, t), has a very different structure than the 1D Poisson's equation used in the original BGK analysis.
Depending on the choices of trapped-particle and untrapped-particle distribution functions,

II. THEORETICAL MODEL AND ASSUMPTIONS
This section provides a brief summary of the one-dimensional kinetic g-factor model (Sec. II A) developed by Davidson and Startsev [14] to describe the nonlinear longitudinal dynamics of a long coasting beam propagating in the z-direction through a circular, perfectly conducting pipe with radius r w . The 1D kinetic Vlasov equation for the distribution function F b (z, p z , t) is used (Sec. II B) to derive several important conservation laws (both local and global) corresponding to conservation of particle number, momentum, and energy per unit length of the charge bunch. The results in Secs. II A and II B form the basis for the nonlinear traveling-wave and traveling-pulse solutions studied in Secs. III and IV.

A. Theoretical Model and Assumptions
This paper makes use of a one-dimensional kinetic model [14] that describes the nonlinear dynamics of the longitudinal distribution function F b (z, p z , t), the average self-generated axial electric field E z (z, t), and the line density λ b (z, t) = dp z F b (z, p z , t) , for an intense charged particle beam propagating in the z-direction through a circular, perfectly conducting pipe with radius r w . For simplicity, the analysis is carried out in the beam frame, where the longitudinal particle motion in (z, p z ) phase space is assumed to be nonrelativistic, and the beam intensity is assumed to be sufficiently low that the beam edge radius r b and rms radius R b = r 2 1/2 = x 2 + y 2 1/2 have a negligibly small dependence on line density λ b . Furthermore, properties such as the number density n b (r, z, t) of beam particles are assumed to be azimuthally symmetric about the beam axis (∂/∂θ = 0), where x = r cos θ and y = r sin θ are cylindrical polar coordinates. Finally, the axial spatial variation in the where ∂/∂z k z L −1 z is the inverse length scale of the z-variation. Making use of these assumptions, it can be shown that the one-dimensional kinetic equation describing the nonlinear evolution of the longitudinal distribution function F b (z, p z , t) and average longitudinal electric field E z (z, t) can be expressed in the beam frame correct to order k 2 z r 2 w as [14] ∂ ∂t and Here, e b and m b are the charge and rest mass of a beam particle, and λ b0 = const. is a measure of the characteristic line density of beam particles, e.g. , its average value. Moreover, the constants U 2 b0 and U 2 b2 have dimensions of speed-square, and are defined by where g 0 and g 2 are the geometric factors defined by [14] In obtaining Eqs. (1)-(5), a perfectly conducting cylindrical wall with E z (r = r w , z, t) = 0 has been assumed.
For purposes of illustration, we consider the class of axisymmetric density profiles n b (r, z, t) of the form Here, λ b = dp z F b (z, p z , t) = 2π rw 0 dr rn b (r, z, t) is the line density, r b is the beam edge radius, assumed independent of λ b , and f (r/r b ) is the profile shape function with normaliza- over the interval 0 r < r b , it can be shown that [14] g 0 = ln r 2 From Eqs. (6)-(8), we note that n = 0 corresponds to a step-function density profile; n = 1 corresponds to a parabolic density profile; n 2 corresponds to an even more sharply peaked density profile; and that the precise values of g 0 and g 2 exhibit a sensitive dependence on profile shape [14]. Finally, for the choice of shape function f (r/r b ) = (n + 1) (1 − r 2 /r 2 b ) n , n = 0, 1, 2, · · · , it is readily shown that the mean-square beam radius

B. Conservation Relations
Equations (1) and (2) possess several important conservation laws, both local and global, corresponding to conservation of particle number, momentum, and energy per unit length.
For present purposes we express E z (z, t) = − (∂/∂z) φ (z, t). Equations (1) and (2) then describe the evolution of F b (z, p z , t) and φ (z, t) according to where v z = p z /m b and Here, is a dimensionles measure of the line density λ b (z, t), and λ b0 = const. is the characteristic (e.g. average) value of line density.
It is convenient to introduce the macroscopic moments where is the average axial flow velocity in the beam frame. Note that the effective 'pressure' P b (z, t) and 'heat flow' Q b (z, t) are defined (relative to the average flow velocity V b ) by and where V b (z, t) is the average flow velocity defined in Eq. (12).
We now make use of Eqs. (9) and (10) to derive the local and global conservation laws corresponding to the conservation of particle number, momentum, and energy per unit length of the beam. The subsquent analysis applies to the two classes of beam systems: (a) a very long, finite-length charge bunch (L b r w ) with N b (z → ±∞, t) = 0; and (b) a circulating beam in a large-aspect-ratio (R 0 r w ) ring with periodic boundary condition as the beam circulates around the ring with major raduis R 0 .
(Here, z can be viewed as the arc length around the perimiter of the ring with large radius Number Conservation: From Eqs. (9), (11) and (12), operating on Eq. (9) with λ −1 b0 dp z · · · , and integrating by parts with respect to p z , we obtain where N b (z, t) = λ b (z, t) /λ b0 is the normalized line density, and V b (z, t) is the axial flow velocity [Eqs. (11) and (12)]. Equation (16) is a statement of local number conservation, i.e. , the time rate of change of the local density, ∂N b /∂t, is equal to minus the derivative of the local flux of particles, − (∂/∂z) (N b V b ). If we integrate Eq. (16) over z, applying the boundary conditions described earlier in this section, we obtain which corresponds to the global conservation of the number of beam particles.

Momentum Conservation:
We now operate on Eq. (9) with λ −1 b0 dp z p z · · · , where p z = m b v z , and make use of Eqs. (12) and (13). This gives where − (∂/∂z) ( φ ) is defined in Eq. (10), and we have integrated by parts with respect to p z to obtain Eq. (18) from Eq. (9). Equation (18) can be expressed in an alternate form by making use of Eq. (10) to eliminate e b (∂/∂z) φ and combine Eqs. (12)- (14) to express where V b (z, t) is the average flow velocity, and P b (z, t) is the effective pressure of the beam particles. Substituting Eqs. (10) and (19) into Eq. (18), we obtain Note that Eq. Eq. (20) over z and applying the boundary conditions described earlier in Sec. II gives which corresponds to global momentum conservation.
Energy Conservation: We now operate on Eq. (9) with λ −1 b0 dp z 1 2 m b v 2 z · · · and make use of Eqs. (11)-(13) and p z = m b v z . Integrating by parts with respect to p z , we readily obtain ∂ ∂z From Eqs. (10) and (16), some straightforward algebraic manipulation gives Substituting Eq. (23) into Eq. (22) and rearranging terms, we obtain which corresponds to local conservation of energy. Global energy conservation follows upon integrating Eq. (24) over z, which gives Note that Eq. (25) describes the balance in energy exchange between particle kinetic energy and electrostatic field energy. Moreover, the final two terms in Eq. (25) correspond to electrostatic field energy, and the term proportional to U 2 b0 is positive, whereas the term proportional to U 2 b2 is manifestly negative. Because of the negative sign of the third term in Eq. (25), note that any increase in (∂N b /∂z) 2 averaged over z must compensated by a corresponding increase in the first two terms in Eq. (25).
To summarize, the local conservation laws in Eqs. (16), (18) and (24), and the global conservation laws in Eqs. (17), (21) and (25), provide powerful nonlinear constrains on the evolution of the normalized line density N b , momentum density N b m b V b , and kinetic en- Furthermore, these conservation constraints are exact consequences of the 1D nonlinear Vlasov equation (9) where e b (∂/∂z) φ (z, t) is defined in Eq. (10), and U 2 b0 and U 2 b2 are expressed in terms of the geometric factors g 0 and g 2 in Eq. (3).
Finally, the energy balance equation (22), the momentum balance equation (18), and the continuity equation (16) can be combined to give a dynamical equation for the evolution of the effective pressure P b (z, t) of the beam particles. We make use of where Q b is the effective heat flow defined in Eq. (15). Without presenting algebraic details, some straightforward manipulation of Eq. (22) that make use of Eqs. (16) and (18) then gives To summarize, Eqs. (16), (20) and (27) In the special case where the heat flow contribution (∂/∂z) Q b is negligibly small in Eq. (27), the pressure P b (z, t) evolves approximately according to The continuity equation (16) can be expressed as Combining Eqs. (27) and (28) , we obtain which can be integrated to give the triple-adiabatic pressure relation (P b /N 3 b ) = const. Therefore, for negligibly small heat flow in Eq. (27), the macroscopic fluid model obtained by taking moments of the 1D Vlasov equation (9) closes, and the nonlinear evolution of N b , V b and P b is described by Eqs. (16), (20) and (30).
In Sec. III, we discuss a particular choice of distribution function F b (z, p z , t), corresponding to the so-called waterbag distribution, for which the heat flow Q b (z, t) is exactly zero during the nonlinear evolution of the system. In this case, the closure is exact, and the nonlinear evolution of the system is fully described by Eqs. (16), (20) and (30).

WATERBAG MODEL
The 1D kinetic g-factor model based on Eqs. (1) and (2) can be used to determine the nonlinear evolution of the beam distribution function F b (z, p z , t) for a broad range of system parameters and initial distribution functions. In this section, we examine Eqs. (1) and (2) for the class of exact solutions for F b (z, p z , t) corresponding to the so-called waterbag distribution in which F b (z, p z , t) has uniform density in phase space (Sec. III A). The subclass of coherent nonlinear traveling-wave and traveling-pulse solutions with undistorted waveform are then examined (Sec. III B) for disturbances traveling in the longitudinal direction with constant normalized velocity M = const.

A. Kinetic Waterbag Model
Equations (1) and (2), or equivalently, Eqs. (9) and (10) 1) and (2) support solutions corresponding to sound-wave-like disturbances with signal speed depending on U b0 and the momentum spread of F b , and cubic dispersive modifications depending on U b2 [14].
In this section, we specialize to the class of exact nonlinear solutions for F b (z, p z , t) to Eq. (1) corresponding to the waterbag distribution [15][16][17][18][19] for −∞ < z < ∞ (long coasting beam in linear geometry) or 0 < z < 2πR 0 (large-aspectratio ring with major radius R 0 ). In Eq. (31), the distribution function , assumed single-valued, of course distort nonlinearly as the system evolves according to Eqs. (1) and (2) [or equivalently, Eqs. (9) and (10)]. We integrate across the two boundary curves in Eq. 31 by operating on Eq. (1) with where p z = m b v z . Integrating by parts with respect to p z , and taking the limit → 0 + , we obtain for the nonlinear evolution of the boundary curves where E z is defined in Eq. (2).
For the choice of waterbag distribution in Eq. (31), we calculate several macroscopic fluid quantities [see also Eqs. (11)-(15)] corresponding to line density axial flow velocity beam particle pressure and beam particle heat flow Note that the heat flow is exactly Q b = 0 for the choice of waterbag distribution in Eq. (31).
Making use of the dynamical equations for V − b (z, t) and V + b (z, t) in Eqs. (33) and (34), where E z is defined in Eq. (2), some straightforward algebra shows that and P b (z, t) evolve according to Note from Eqs. (37) and (41) that P b (z, t) can be expressed as where P b0 = const. and λ b0 = const. represent the characteristic (e.g. , average) value of the pressure and line density of the beam particles respectively, and P b0 /λ 3 b0 = 1/12 (m b A) 2 = const. , where A is the constant phase-space density in Eq. (31). By virtue of the fact that the is not surprising that Eqs. (39)-(42) are identical to the macroscopic fluid equations (16), (20) and (30), obtained in Sec. II B, where Eq. (30) has made the assumption of negligible heat For present purposes, we introduce the effective thermal speed U bT associated with the waterbag distribution in Eq. (31) defined by and the nomalized (dimensionless) fluid quantities η (z, t) and U (z, t) defined by In Eq. (44), (U 2 b0 + U 2 bT ) 1/2 is the effective sound speed associated with the geometric factor g 0 and the thermal speed U bT . Furthermore, we introduce the scaled (dimensionless) time variable T and spatial variable Z defined by The fluid description in scaled variables provided by Eqs. (46) and (47) is exactly equivalent to the kinetic description provided by Eqs. (1) and (2) for the choice of waterbag distribution in Eq. (31).

B. Coherent Nonlinear Traveling Ware and Pulse Solutions
Within the context of the present 1D model, Eqs. (46) and (47) can be used to investigate detailed properties of collective excitations over a wide range of system parameters. For example, in the weakly nonlinear regime, for small-amplitude disturbances moving near the sound speed (U 2 b0 + U 2 bT ) 1/2 , Eqs. (46) and (47) can be shown to reduce to the Korteweg-deVries equation [15], which exhibits the generation and interaction of coherent structures (solitons) for a wide range of initial density pertuibations η (Z, T = 0) = 0 [22]. While the analysis in Ref [15] has several interesting features, the results are limited to the weakly nonlinear regime where |η| 1 and |U | 1.
In this paper, we examine Eqs. (46) and (47) ∂ ∂Z Integrating with respect to Z , Eqs. (48) and (49) give Substituting Eq. (52) into Eq. (53), we obtain which is a second-order nonlinear differential equation for the perturbation in line density where V (η) is the effective potential defined by and the dimensionless parameter T , defined by is a measure of the longitudinal thermal speed of the beam particle.
Equations (55) and (56)  In general, the effective potential V (η) in Eq. (57) can be expressed as where In Eqs. (58) and (59), T is restricted to the range 0 < T < 1/3, and M 2 can satisfy M 2 > 1 or M 2 < 1. Examination of Eq. (59) shows that for all allowed values of T and M 2 . Furthermore, it's also clear from Eq. (59) that  In this case, the effective potential has the qualitative shape illustrated in Fig. 6, which has been plotted for the choice of parameters M 2 = 9 and T = 1/30. The physically allowed, localized pulse solutions (soliton solutions) corresponds to the energy level which is the red horizontal line in Fig. 6, and boundary conditions In the special circumstances where M 2 exceeds 1 by a small amount, i.e. , M 2 = 1 + ∆ where 0 < ∆ 1. it is readily shown that Eq. (54) can be approximated for small η by Equation (67) can be solved exactly for η (Z ) = λ b (Z ) /λ b0 − 1 to give Note that the soliton amplitude in Eq. (68) is small for M 2 = 1 + ∆ with ∆ 1. Also, the sech 2 {· · · } pulse shape in Eq. (68) is similar to the soliton pulse shape obtained from the Korleweq-deVries equation in the weakly nonlinear regime [15]. The present analysis of Eqs. (9) and (10) parallels the original Bernstein-Greene-Kruskal (BGK) formulation of BGK solutions to the 1D Vlasov-Poisson equations [25,26], except for the fact that Eq. (10), which connects φ (z, t) to the line density λ b (z, t), has a very different structure than the 1D Poisson equation.
Referring to Eqs. (9) and (10), we introduce the scaled dimensionless variables (Z, P z , T ) defined by where U 2 b0 and U 2 b2 are defined in Eq. (3), and U 2 bT = const. is the longitudinal velocity spread characteristic of the distribution function F b . We futher introduce the dimensionless distribution functionF b (Z, P z , T ) defined bŷ where λ b0 = const. is the characteristic line density of the beam particles, e.g., the average value. From Eqs. (69), (70) and the definition of line density λ b = dp z F b , it follows that the perturbation in line density η = λ b /λ b0 − 1 can be expressed as where the P z integration covers the range −∞ < P z < ∞ in Eq. (71). Transforming variables according to Eqs. (69) and (70), and making use of Eq. (71), it is readily shown that Eqs. (9) and (10) can be expressed in the new variables as and where ψ (Z, T ) is the normalized (dimensionless) potential defined by Equations (72) and Here, P z = V z , andF b (Z , P z , T ) and η (Z , P z , T ) are related by Therefore, the traveling-pulse or traveling-wave solutions that have stationary profile shape in the primed variables (Z , P z , T ) are determined by setting ∂/∂T = 0 in Eqs. (75)-(77).
where ψ (Z ) and η (Z ) solve Eq. (76), and η (Z ) is related toF b (Z , P z ) by Eq. (77). We introduce the energy variable W defined by Then the solution to Eq. (78) forF b (Z , P z ) can be expressed aŝ where Note from Eq. (79) that where + corresponds to V z > 0, and − corresponds to V z < 0. Substituting Eqs. (80) and (81) into Eq. (77) gives which relate the perturbation in beam line density η (Z ) to the potential ψ (Z ) and the distribution functionsF > b (W ) andF < b (W ). Figure 11 shows an illustrative plot of the potential ψ (Z ) as a function of Z . Depending on the values of the energy W and range of Z , there are three classes of particle orbits: (a) particles that are reflected from the potential; (b) particles that are trapped and undergo periodic motion; and (c) passing (untrapped) particles that don't change direction, but pass over the potential maximum, first slowing down and then speeding up during the motion.
For the trapped particles and the reflected particles, it follows thatF > (W ) =F < (W ) so andF On the other hand, for the passing (untrapped) particles,F > U n (W ) andF < U n (W ) can be specified independently, depending on whether the particles have V z > 0 or V z < 0, respectively.
The form of ψ (Z ) shown in Fig. 11 corresponds to a stationary isolated pulse in primed variables, with ψ (Z = ±∞) = 0. By contrast, Fig. 12 shows a plot of ψ (Z ) versus Z for the case where ψ (Z ) has a periodic nonlinear wave structure with ψ (Z + L) = ψ (Z ) .
(86) From Fig. 12, trapped particles with energy W in the range exhibit periodic motion. On the other hand, passing particles with energy W in the range (see Fig. 12) correspond to untrapped particles that pass over the potential ψ (Z ), periodically speeding up and slowing down, but not changing their direction of motion. Furthermore, for the nonlinear periodic waveform for the potential ψ (Z + L) = ψ (Z ) shown in Fig. 12, it follows from Eqs. (76) and (83) that the waveform for the perturbation in line charge also satisfies η (Z + L) = η (Z ). Here, η (Z ) is related to ψ (Z ) and the trapped-particle and untrappedparticle distribution functions by Eq. (83), which gives In Eq. (89), the integration over the trapped-particle distributuionF T r (W ) is over the interval of W corresponding to ψ min < ψ < W < ψ max , and the integration over the untrapped particle distributionF U n (W ) is over the interval of W corresponding to ψ max < ψ < W < ∞.
Equations (76) where W U = const. , A = const. , and W U > ψ max (see Fig. 12). SubstitutingF T r (W ) = 0 and Eq. (90) into Eq. (89) readily gives which can also be expressed as where Note that Eq. (94) has the form of a dynamical equation of motion, with η playing the role of displacement, Z playing the role of time, and V (η) playing the role of an effective potential.
Making use of Eq. (95), it is readily shown that and where the constant of integration in Eq. (97) has heen chosen so that V (η = 0) = 0.
Close examination of Eqs. (93)-(96) shows that Eq. (94) supports oscillatory solutions for When the inequality in Eq. (98) is satisfied, the plot of V (η) versus η has the characteristic shape illustrated in Fig. 13 for the choice of parameter 2W U < 1. Here, V (η) has a minimum at η = 0, and passed through zero at

V. CONCLUSIONS
In this paper, the 1D kinetic model developed in Ref. [14] was used to describe the non-