A procedure to obtain a near equilibrium phase space distribution function has been derived for beams with space charge effects in a generalized periodic focusing transport channel. The method utilizes the Lie transform perturbation theory to canonically transform to slowly oscillating phase space coordinates. The procedure results in transforming the periodic focusing system to a constant focusing one, where equilibrium distributions can be found. Transforming back to the original phase space coordinates yields an equilibrium distribution function corresponding to a constant focusing system along with perturbations resulting from the periodicity in the focusing. Examples used here include linear and nonlinear alternating gradient focusing systems. It is shown that the nonlinear focusing components can be chosen such that the system is close to integrability. The equilibrium distribution functions are numerically calculated, and their properties associated with the corresponding focusing system are discussed.

## I. INTRODUCTION

There is wide interest in the study of periodically focused beams with space charge intensities high enough that self-field effects become important.^{1} Interest has also arisen in the study and applications of nonlinear focusing, some of which involve intense beams.^{2–11} Equilibrium and quasi equilibrium distributions are useful in understanding the transport of intense beams with respect to stability analysis, in preparing beams that are well matched or to perform numerical simulations of beams that are perturbed from this matched condition.^{12} Various equilibrium properties of beams have already been studied^{13–23} for linear focusing systems. The only known exact beam equilibrium for periodic focusing systems is the Kapchinskij-Vladmirskij (KV) distribution,^{24} which was recently generalized for coupled lattices.^{25} This distribution function is unphysical with the density having a discontinuity at the edge and being exactly uniform in the interior of the beam. Moreover, the KV distribution involves expressing the distribution as a *delta* function of the Courant-Snyder (CS)^{26} invariants or the CS invariants generalized to coupled lattices,^{27} both of which are valid for only linear focusing systems.

Perturbation methods have been employed to perform canonical transformations in finding an averaged Hamiltonian.^{23,28,29} The averaged Hamiltonian allows one to look for more realistic distribution functions that are close but not exactly at equilibrium. These methods were developed primarily for application to linear focusing systems. In this paper, we utilize the Lie transform method to perform such a perturbation analysis. The Lie transform approach exploits the invariance of Poisson brackets under canonical transformations. This becomes useful when dealing with Hamiltonians more complex than linear focusing systems, such as ones used in this paper. In the Lie method, transforming between an averaged coordinate frame and one with higher order oscillating components may be easily accomplished with the help of the so called “push forward” and “pull back” operators. In the procedure described in this paper, we first proceed to obtain a solution in the averaged coordinate system yielding a function independent of “s,” the axial distance or the time like variable. Following this, the transformation operators are used to recover the oscillation of the beam along “*s*,” which arises due to the periodicity of the focusing system. The distribution functions obtained are then used to compute other associated quantities, namely, the number density and flow velocity functions.

Hamiltonian perturbation theory has been used to determine a configuration for nonlinear lattices^{30} with improved integrability. It was shown there that a condition exists, wherein this nonlinear lattice can be designed in such a way that the averaged Hamiltonian is azimuthally symmetric in the absence of space charge forces. In this paper, we extend this analysis to include space charge forces and in the process show that the same conditions of optimum integrability derived in Ref. 30 are retained in the presence of a quasi equilibrium distribution. We restrict ourselves to computing quasi equilibrium distributions in systems that have such an averaged, azimuthally symmetric Hamiltonian. It has been previously shown^{30–32} that maximum integrability is desirable because it leads to improved confinement, or what is known as increased dynamic aperture. This concept formed the basis of the proposal of the “Integrable Optics Test Accelerator” (IOTA).^{33} The combination of optimum integrability and a distribution close to equilibrium is favorable for greater acceptance of the transport system as well as a stable beam with respect to mismatches, which would otherwise cause initial emittance growth.

This paper is organized as follows. The underlying concept and a brief outline of the analysis is discussed in Sec. II. In Sec. III, this analysis is applied to a specific class of periodic focusing systems, and in the process, explicit expressions for the distribution function, number density, and flow velocity distribution are derived, and finally the results are applied to alternate gradient focusing systems consisting of quadrupoles and/or sextupoles. In Sec. IV, a set of dimensionless parameters are described that fully describe the system scaled with respect to the length of the lattice period. Additionally, a set of conditions are defined that bring about an equivalence between a linear and a nonlinear focusing system. The number density and flow velocity distributions are numerically computed for a thermal equilibrium case, and their properties are discussed for a linear focusing case in Sec. V and an equivalent nonlinear focusing system in Sec. VI. All the results of this paper are summarized in the final section.

## II. DESCRIPTION OF THE CANONICAL TRANSFORMATION PROCEDURE

This section gives a brief outline of the canonical transformation procedure to be used in the paper. We will be dealing with a system that represents a thin long beam with a constant axial velocity and no longitudinal oscillations. Under such conditions, the axial distance may be treated as the independent variable. The Lie transformation is then defined with respect to a phase space function *w* such that it satisfies the following Poisson bracket relationship:

where $Z=(P,Q)$ is a phase space vector representing the generalized positions and momenta of the system, *w* is the Lie generating function, *s* is the independent variable, in this case the axial distance, and *ϵ* is a continuously varying parameter such that $Z(\u03f5=0)=z$, where $z=(p,q)$, the original phase space vector. The above relationship resembles Hamilton's equation with respect to a “Hamiltonian” *w*, and “time” *ϵ*. This guarantees that the transformation is canonical for all values of *ϵ*.

The Lie operator *L* is defined such that it performs a Poisson bracket operation with respect to *w*. Symbolically, this is represented as

A transformation operator *T* is defined such that its role is to replace the variables of a function by the new canonical variables. For example, for a function $g(z)$,

For the identity function, this is simply

*T* and $T\u22121$ are sometimes referred to as the “push forward” and “pull back” operators, respectively.

The method used in this paper is based on previous work on Lie transformation theory for performing Hamiltonian perturbation analysis.^{30,34–39} It involves using *ϵ* as an order parameter, where the original Hamiltonian is first expressed as

This Hamiltonian may then be transformed to

where the terms *K _{n}* are determined up to a desired order. The parameter

*ϵ*is used only for the purpose of keeping track of the terms and the respective order they belong to. Thus,

*ϵ*can be set to unity in the final expressions.

The analysis begins with an expansion^{36} of the operators *T*, and *L*, and the function *w* such that

where

The expressions for the perturbation expansion terms in orders of *ϵ* for *T* have been derived in Ref. 39 and are given by

The formulation involved in obtaining the transformed Hamiltonian, which is performed order by order involves determining terms of the transformed Hamiltonian *K* and the Lie generating function *w*. The relationships between *K* and *w*, derived in Ref. 39, are given here up to third order, which are

It may be noted that in these equations, *K* is initially expressed as a function of ($q,p$), the original variables. However, based on the canonical invariance of the whole formulation, ($q,p$) may be simply replaced by ($Q,P$), the transformed variables in the final expression for *K*.

In the study of self consistent systems, such as beams with space charge effects, it is often required to solve for the phase space density $f(q,p)$, which satisfies the Vlasov equation

In many cases, such as systems with electrostatic space charge forces, it is possible to decompose *H* as

where *H ^{e}* is the term arising from external forces and

*H*arises from space charge forces and is itself a function of

^{sc}*f*.

In Vlasov-Poisson systems, *H ^{sc}* is proportional to the electrostatic potential $\varphi $ which satisfies Poisson's equation

where *q* is the charge of the particle and *ϵ*_{0} is the free-space permittivity. The ordering parameter, *ϵ*^{m}, has been inserted to allow for what is known as maximal ordering. Since we are solving for beam equilibria, this factor will be selected so that the lowest-order space charge potential *H ^{sc}* is of the same order as the lowest-order averaged confining forces. This is necessary for an equilibrium condition as otherwise the electrostatic forces would cause beam divergence.

If we define a function

then it may be easily verified from Eq. (3) that

From the property of canonical invariance of the Vlasov equation (see, for example, Ref. 40), we then have

Using the perturbation expansion terms of the transformation operator *T* given by Eq. (7), we can in principle expand *f* in terms of *F* using Eq. (17) in order to obtain the phase space distribution in the original coordinate system. This gives

where

From the analysis in Sec. III, it will be clear that the different order terms represent oscillations at different spatial frequency scales in “*s*,” with the lowest order term being independent of “*s*.” The ordering convention is chosen such that the zeroth term in *f* corresponds to the lowest order nonzero term in the expansion. We solve for *F* to be independent of “*s*,” corresponding to an equilibrium. Thus, it is clear that *F* will not have terms oscillating in “*s”* at multiple scales. As a result, we can allow *F* to belong to a single order. Other systems that produce multiple terms in *F* have been analyzed in the past, for example, the self consistent, quasi-static response to ponderomotive effects.^{41}

Combining Eqs. (19) and (15), $\varphi $ may be expressed as a decomposition of its different order terms such that

where

The different terms in the electrostatic potentials represent contributions arising from the different orders in the perturbation of the charge density. Thus, the lowest nonzero term in $\varphi $ corresponds to the lowest nonzero term in *f*.

The motivation to carry out the procedure outlined in this section is based on the presumption that it is possible to perform a canonical transformation from $(q,p,s)$ to $(Q,P,s)$ such that the resulting phase space function, $F(Q,P,s)$ is easier to solve for than $f(q,p,s)$. It is then more convenient to first solve for *F* and then transform back to *f* rather than solving directly for *f*. In this study, we solve for an equilibrium, where *F* is independent of *s*. This is achieved through a transformation that depends on *s*. If the transformation is perturbative, which is true in this context, *f* may be obtained only up to a desired order. It has been shown in Appendix A that the terms higher than the lowest nonzero term in *f* do not change the total number of particles in the distribution. Thus, the total charge is conserved even upon truncation of the series expansion in *f*.

## III. APPLICATION TO PERIODIC FOCUSING SYSTEMS

The analysis described in Sec. II is now applied to the problem of finding a solution for a specific class of systems. We do this in two steps. In the first step, we obtain the averaged Hamiltonian $K(Q,P)$ which is independent of “*s*”. Thus, any distribution function that is a function of *K* is an equilibrium in the $(Q,P)$ space. In the second step, we show how to recover the “*s*” dependence of the beam by transforming the distribution function $F(Q,P)$ to $f(q,p,s)$. Using the expressions for *f*, we then derive expressions for the number densities and the velocity vectors of the flow in position space.

### A. Obtaining the averaged Hamiltonian

For the sake of brevity and generality, we shall express the Hamiltonian in this section as

We perform the analysis on an explicit form for the Hamiltonian of such a system expressed in cylindrical coordinates later in this section. We thus restrict our analysis to systems whose Hamiltonians can be decomposed into a “Kinetic energy” term $A(q,p)$, an “external potential” term $U(q,s)$, and a “space charge term” $\varphi (q,s)$. In Cartesian coordinates, $A$ is simply $1/2(px2+py2\cdots )$, while in other curvilinear coordinate systems an equivalent expression for $A$ typically depends on $q$. We make further assumptions on *A* and *U* at a later stage. The term *ξ* is a constant factor that multiplies the space charge potential and will be expressed explicitly in Sec. III C.

In the analysis to follow, we assume that the lowest nonzero term of the Hamiltonian is proportional to *ϵ*. This ordering enables one to perform a canonical transformation to slowly varying coordinates and has been used before.^{30,40,42} Based on this ordering scheme we have $H0=0$ and

and further, $H2=\xi \varphi 2,\u2009H3=\xi \varphi 3$, etc. We shall assume that *U* is periodic over *s* such that

The value of *m* appearing in Eq. (15) is yet to be determined. Besides, it may be noted that $\varphi n=0$ for all *n* < *m*. Thus, we already know that that *m* > 0 since $H0=0$, otherwise there would be a zeroth order contribution by the electrostatic potential to the Hamiltonian. With this ordering scheme, the perturbation expansion results in a powers series of the ratio between the fast and slow oscillation time periods. The “fast” time period in this case being *S* and the “slow” one will be a periodicity corresponding to the as yet undermined averaged Hamiltonian. Applying the zeroth order transformation given by Eq. (12a) and using $H0=0$ simply gives

For convenience, we decompose $U(q,s)$ into $\u3008U\u3009$ a term averaged over the period *S* and a purely oscillatory term $U\u0303$, such that

where $\u3008\cdots \u3009=1/S\u222bss+S(\cdots )ds$ denotes an average over the period *S*. For the purpose of this paper, we impose another condition on *U*, which is,

This condition is not unusual for a periodic focusing system in a linac with equal focusing and defocusing components such as an alternating gradient system consisting of quadrupole magnets. A storage ring would have an asymmetry in the average external focusing depending upon the configuration of the dipoles. In such cases, one may need to retain a nonzero $\u3008U\u3009$.

The expression for the first order term of the Hamiltonian given by Eq. (24) will have to be used in Eq. (12b) to obtain the terms *K*_{1} and *w*_{1}. This leads to

It may be noted that this equation has two unknown quantities, *K*_{1} and *w*_{1}. Choosing one will determine the other accordingly. We intend to eliminate fast oscillating terms in the transformed Hamiltonian and the term *K*_{1} is chosen accordingly. This rule in determining *K* will be followed for all higher orders too. Additionally, it is necessary that the expression for *w* be purely oscillatory in order to prevent the perturbation scheme from becoming secular (unbounded) in *s.*^{39} Thus, we must require that *K* should cancel the terms with a nonzero average over *S* in the expressions of *w* for every order. We know that by definition $\u3008U\u0303\u3009=0$ and so it cannot be a part of *K*_{1}. Thus, its easy to see that there is no external focusing term in the transformed Hamiltonian for this order. Note that this would not be true if Eq. (28) was not satisfied. In this problem, we are seeking an equilibrium solution and thereby will require that the lowest order space charge term be of the same order as the lowest order external potential appearing in the transformed Hamiltonian. As a result, we need to set *m* > 1 and so $\varphi 1=0$. By applying all the above conditions, we then obtain

Inserting this into Eq. (29) and integrating with respect to *s* yields

where the Roman numerical superscript represents an integration over *s* with an additive constant chosen such that

This is equivalent to an indefinite integral of a Fourier expansion of the purely oscillatory function $U\u0303$. This condition on *w*_{1} will be extended to all higher orders of *w*. Similarly, a superscript “II” represents a double integral with the same condition and so on. This convention will be used throughout this paper. In addition, a subscript of $(q,p)$ hereafter represents a gradient with respect to the corresponding phase space variable. For example, the phase space vector, $Ap=\u2202A/\u2202p$. Further, $Apq$ corresponds to the second order tensor in phase space, $\u22022A/\u2202p\u2202q$ and so on.

Applying the first order results to the second order equation, Eq. (12c) gives

The quantity $U\u0303qI\xb7Ap$ has a zero average over the period *S* because it is linear with respect to $U\u0303$ while the term *A*, by definition has no oscillating components. Imposing the restriction that *K*_{2} does not include oscillating terms, and *w*_{2} must be purely oscillatory, we see that *K*_{2} does not contain an external force term. Based on the fact that an equilibrium distribution is sought, we thus need to have $\varphi 2=0$, and consequently *m* > 2. Thus, we obtain

Integrating Eq. (33) to obtain *w*_{2} gives

where, as mentioned earlier, the superscript $II$ represents a double integral with an additive constant such that $\u3008U\u0303II\u3009=0$, leading to the necessary condition $\u3008w2\u3009=0$.

We now proceed to evaluate terms of the next order in *ϵ*. Using the expressions already derived for the lower order terms and applying them to Eq. (12d) gives

As usual, *K*_{3} is chosen so that it cancels all components which remain nonzero when averaged over the period *S* in Eq. (36). Imposing this condition, and recognizing the terms that do average to zero over the period *S*, we get

where we used $\u3008U\u0303qII\xb7App\xb7U\u0303q\u3009=\u2212\u3008U\u0303qI\xb7App\xb7U\u0303qI\u3009$ from integration by parts. Here, we require that the electrostatic potential term be a part of *K*_{3} because it contains a non zero external focusing component. This leads us to *m* = 3. It remains to be shown that $\varphi 3$ does not contain fast oscillating components. Applying Eq. (20) and recognizing the fact that *T*_{0} is the identity operator, we have

where we seek *F* to be independent of *s*. Consequently, from Eq. (22) with *m* = 3, we get

It is clear from Eq. (39) that $\varphi 3$ does not average to zero over the period *S*, since *F* in this case is independent of *s*. Thus, we need to choose $\varphi 3$ to be a part of *K*_{3}.

Evaluating *w*_{3} explicitly will involve integrating all terms of Eq. (36) in the manner similar to previous cases, that is, imposing the condition that the integral has a zero average over fast oscillations. Let us first define a term $G\u0303$ such that

It is easy to see from this expression that $G\u0303$ is periodic over *S* with a zero average, or $\u3008G\u0303\u3009=0$. The function *w*_{3} after integrating Eq. (36) in the usual way, can now be expressed more compactly as

It is straightforward to verify that the required condition, $\u3008w3\u3009=0$ is satisfied.

Thus far, we have performed the task of deriving the generating function *w* and the transformed Hamiltonian *K* up to third order. Performing the analysis up to third order was necessary in order to get the fundamental effects of the actual system. We state that this is sufficient for the specific application of the paper. Obtaining higher order terms would involve carrying out the same procedure, while the expressions become longer.

The transformed Hamiltonian may be expressed up to third order by

where, as mentioned earlier, the parameter *ϵ* is now set to one. Combining the terms given by Eqs. (30), (34), and (37), respectively, we get

### B. Deriving the distribution function and some other associated quantities

In this subsection, we derive expressions for the distribution function *f*. Using this, we obtain expressions for the number density and the flow velocity vector both as functions of *q* and *s*. It is easier to relate the latter two quantities with the physical nature of the external potential than the phase space distribution function.

The averaged Hamiltonian given by Eq. (43) has an external potential independent of *s* allowing one to compute an equilibrium distribution, which is a steady state solution of Eq. (18) given by

Transforming *F* to *f* up to the desired order based on the procedure described in the Sec. III A would then yield a near equilibrium distribution function in the original coordinate system. The zeroth order term is given by Eq. (38) where it is now established that *F* is independent of *s*. We may determine *f*_{1} by using Eq. (20), along with Eqs. (31) and (11b). We get

We can continue this to one more higher order by using the transformation given by Eq. (11d) with the generating functions, given by Eqs. (31), (35), and (41). This gives

Although expressions for all the required terms in the above equation are derived in this paper, the expression for *f*_{3} shall not be expanded here explicitly for the sake of brevity and since we do not use it.

The number density *n*, or the number of particles per unit volume in position space is the integral of the distribution function over momentum space. The zeroth order term, after referring to Eq. (38), is simply

Using Eq. (45), it is clear that the first order term vanishes as shown below

where we have performed an integration by parts and it is assumed that $F\u21920$ as $p\u2192\u221e$. We again perform an integration by parts to determine *n*_{2} along with the following assumptions, (1) $F\u21920$ as $p\u2192\u221e$ rapidly enough that any product between *F* and a polynomial of $p$ would vanish as $p\u2192\u221e$. This would be true for any physical distribution. The second assumption would be, (2) the “Kinetic energy” term *A* is a second order function in momentum, and so any second derivative of *A* with respect to a momentum component is independent of the momentum variable. Using Eq. (46), we then obtain

Using the relationship $f0(p,q)=F(p,q)$ from Eq. (38) in Eq. (22) and integrating the expression in momentum space, we see that the lowest order nonzero space charge potential term satisfies

We now proceed to calculate the flow velocity vectors or equivalently the particle flux as a function of position. This gives an idea of the motion of the particles affected by the external potential. This motion is separate from the thermal motion of the individual particles, which averages to zero in the case of a uniform and isotropic temperature distribution. We know from Hamilton's equations, and the form given by Eq. (23), that is, $dq/ds=q\u0307=Ap$. The velocity vector $v$ is then given by

While introducing the normalization factor *n*_{0}, we have made use of the fact that $n1,n2,n3\u2026$ do not contribute to total mass because $f1,f2,f3\u2026$ integrate to zero as shown in Appendix A.

From Eq. (30), it can be seen that the quantity *A* is first order in *ϵ* and so the zeroth order term of $v,\u2009v0$ vanishes. This is consistent with the fact that *n*_{0} does not evolve with respect to *s* and so the net particle flux at all values of $q$ must vanish for this order. The same should be true for $v1$ because $n1=0$. This can be verified from the following argument. By matching the orders in *ϵ*, we obtain

where we know that $f0=F$. Since the external potential in this system is purely a function of position, the equilibrium distribution *F* is symmetric about $p$. On the other hand, *A* is second order in $p$ and so $Ap$ is asymmetric about $p$. Thus, the product of these two quantities is also asymmetric about $p$ and must vanish when integrated over $p$. We should expect that $v2$ is nonzero, which may be determined as follows:

where Eq. (45) has been made use of and the usual conditions applied to *F* as $p\u2192\u221e$ while integrating by parts. It is easy to verify that the continuity equation for this order

is satisfied by substituting the above expression for $v2$ and that for *n*_{2} given by Eq. (50), along with the fact that $v0=0$.

### C. Application to alternating gradient focusing systems

We work with a system consisting of a long coasting beam subjected to an alternate gradient focusing system, composed of quarupoles or a combination of quadrupoles and sextupoles. The Hamiltonian for a particle in such a beam can be represented in cylindrical coordinates as shown in Ref. 30 as

where the position on a plane transverse to the direction of propagation is described by the coordinates $(r,pr)$, the radial position and momentum and $(\theta ,l)$, the azimuthal angle and angular momentum. Also, *m*_{0} is the rest mass, *q* is the charge and *v _{s}* is the axial velocity of the particle, and $\gamma =1/1\u2212vs2/c2$. We neglect the contribution of the transverse velocity to

*γ*. The values of $\kappa 2(s)$ and $\kappa 3(s)$ depend upon the strength of the quadrupole and sextupole magnets, respectively, and also the velocity of the particle in the axial direction. The angle

*α*is determined by the orientation of the sextupoles with respect to the quadrupoles. It is assumed that

*κ*

_{2}and

*κ*

_{3}are periodic in

*s*with periodicity

*S*, i.e., $\kappa 2(s+S)=\kappa 2(s)$ and $\kappa 3(S+s)=\kappa 3(s)$. and have a zero average over

*S*. That is,

and the same for *κ*_{3}.

The Hamiltonian in Eq. (56) has the form given by Eq. (23) with $A=12(pr2+l2/r2),\u2009U=12\kappa 2(s)r2\u2009cos(2\theta )+13\kappa 3(s)r3cos(3\theta +\alpha ),\u2009\xi =q/m0vs2\gamma 3.$ Applying Eq. (43) to the Hamiltonian given by Eq. (56), we get

As before, the superscript $I$ represents an integration with respect to “time,” *s* with the same conditions as described in Sec. III A.

Equation (58) has only one term which depends on the azimuthal variable Θ. Thus, the transformed Hamiltonian may be rendered azimuthally symmetric if the lattice functions *κ*_{2} and *κ*_{3} satisfy the condition

This condition is satisfied by alternating gradient focusing systems, where *κ*_{2} and *κ*_{3} are step functions alternating in sign along with a phase lag of a quarter lattice period between each other as shown in Fig. 1. Although we use this as a simple example, a method of imposing optimum integrability to more general cases has been described in Ref. 30. In the rest of this paper, we restrict to cases where Eq. (59) is satisfied. Thus, the external potential, which is independent of “*s*” will be

In addition, we define the total potential or the term of the Hamiltonian that depends purely on the position coordinates as

We can now use Eqs. (48), (49), and (50) to obtain the number densities corresponding to a solution $F=F(K)$. Doing so gives us

and

In Eq. (64), we have a term that is proportional to the derivative of the zeroth order number density with respect to the radial coordinate *r*. In highly space charge dominated beams, we know that the density falls off sharply as one approaches the edge of the distribution. This would lead to rather large values for the ratio $n2/n0$, where the perturbation expansion may not converge. In such cases, one should consider the expansion method valid only in a region where the variation of *n*_{0} with respect to *r* is smooth enough. This is not a crucial limitation as most of the charge occurs within this region. The proper shape of the edge of the beam, where the density tapers off rapidly does not affect the distribution in the core of the beam. The particles constituting the edge of the beam, however, are potential candidates for halo formation, which are influenced by the oscillations within in core of the beam.^{43}

The lowest nonzero velocity vector components may be obtained from Eq. (55). They are

and recognizing that $v\theta =r\theta \u0307$ we have

## IV. REQUIRED PARAMETERS FOR CALCULATING DISTRIBUTION FUNCTIONS

In order to compute the solutions for a system, it is first useful to define a certain set of parameters that completely define the system. In this paper, we are correlating a periodic focusing structure with an equivalent constant focusing one. For this reason, it becomes more convenient to express all quantities relative to *S*, the lattice periodicity. All parameters may be scaled accordingly when applying to a real physical system.

We first begin by describing the quantities for a linear focusing channel, which consists of quadrupoles as the only focusing elements. The fractional magnet occupancy $\Delta s/S$ is defined as the length occupied by the magnetic field of the focusing element relative to the length of the lattice period. The phase advance according to the averaged Hamiltonian of an isolated particle per lattice period, in this case would be,

The two quantities, $\Delta s/S$ and *k*_{2}, completely define the motion of a single isolated particle scaled with respect to *S*, the length of a lattice period for a linear focusing system. We shall now proceed to define the parameters that would completely describe the properties of the beam in the averaged frame of reference. Consider a round continuous beam in a smooth, azimuthally symmetric linear focusing channel. Let us consider a beam with a uniform charge distribution, often referred to in this case as a KV beam. The evolution of the beam radius for such a system can be described by the well known envelope equation (see, for example, Ref. 1). We denote the *rms* value of the radial position *r* over the distribution $n0(r)$ as denoted as *R _{rms}*. We define $R=Rrms/S$ as the rescaled beam size with respect to the lattice period, $N=s/S$ as the rescaled axial distance or simply the number of lattice periods and $\epsilon \u2032=\epsilon /S$ as the rescaled beam emittance with $\epsilon =r2\xaf\u2009pr2\xaf\u2212rpr\xaf2$ where the bar indicates an average over the distribution. Using a simple rescaling of variables from the text book version, the envelope equation, in this system of units may be expressed as

where $Q=q\lambda /2\pi \u03f50m0vs2\gamma 3$ is the usual dimensionless space charge perveance, with *λ* being the linear charge density of the distribution. The tune depression ratio *η* is defined as the ratio between the oscillation frequency of an isolated particle and the oscillation frequency of a particle that is always lying within the beam with a uniform distribution. Using the expressions given in Ref. 44, this may be expressed as $\eta =1+u2\u2212u$ with $u=Q/8k2\epsilon \u2032$. The applicability of the envelope equation is far more general than to that of a round KV distribution. As shown in Ref. 45, it can be generalized to any distribution with charge density constant over elliptical surfaces, and where the envelope equation predicts the evolution of the *rms* size along the particular axis, with the same expression for *Q*, the space charge perveance. Thus, it is not uncommon to describe a non-KV distribution with respect to an equivalent KV beam in an rms sense. While a tune depression ratio does not exist for nonuniform distributions, the expression for *η* is frequently used to provide a measure of the space charge intensity for such cases. By drawing an analogy with a constant focusing system, we have shown that for a linear focusing case, the quantities $\Delta s/S$, *k*_{2}, *Q*, and *η* along with a given beam distribution function would fully describe the system scaled with respect to *S*, the lattice period.

With a nonzero value for *κ*_{3}, the focusing becomes nonlinear and one cannot define quantities such as phase advance per lattice period and tune depression ratio. In order to choose a parameter set that is somewhat equivalent to a purely linear focusing case, we prescribe a set of conditions a nonlinear focusing system must satisfy with respect to the equivalent linear focusing one. As assumed in the analysis of Sec. III, the periodicity *S*, of *κ*_{2} and *κ*_{3} will be set to be the same in all cases studied in this paper. This value of *S* will be the same as that of the quadrupoles in the equivalent, purely linear focusing system. We require that the fractional magnet occupancy $\Delta s/S$ for the quadrupoles and the sextupoles in the nonlinear channel, individually be equal to half that of the quadrupoles in the equivalent, purely linear focusing channel. Thus, the net occupancy fraction of both the nonlinear and equivalent linear focusing case would be equal. The space charge perveance *Q* and the *rms* radius *R _{rms}* will be required to be the same in the nonlinear focusing channel and the equivalent linear focusing channel. To represent the focusing strength from the sextupoles, we define a dimensionless quantity

The averaged force on the particle from the sextupoles and quadrupoles may be obtained from the expression of the Hamiltonian given by Eq. (58). Expressing this in terms of the dimensionless system of units defined in this section, which is scaled with respect to *S*, we get

We require that the net force at $r=Rrms$ be the same in both the systems. Suppose the value of the focusing strength in the equivalent linear focusing channel is $k2\u0303$, then using the above equation with $r=Rrms$ this condition leads to the equation

Finally, we require that the averaged force produced by linear component, due to quadrupoles in the nonlinear focusing system is 4 times that produced from the component due to sextupoles, at $r=Rrms$. This would mean

These conditions provide a measure on how much the linear force is to be reduced when the nonlinear component is introduced. A similar set of conditions on the extent of the nonlinearity in the external focusing, in comparison with an equivalent linear focusing system, was used by us in Ref. 9 for studying the control of beam halo formation through nonlinear focusing.

## V. CALCULATIONS OF DISTRIBUTIONS WITH PURELY LINEAR FOCUSING

We now calculate the number density for the case of $\kappa 3=0$, corresponding to an alternate gradient focusing system, provided solely by quadrupoles. The averaged equilibrium number density *n*_{0} is calculated based on Eq. (62) for a given expression for *F*. Additionally, the higher order term *n*_{2} is calculated from Eq. (64), and the lowest order term in the flow velocity **v** is calculated from Eqs. (65) and (66). While *n*_{0} is independent of *s*, *n*_{2}, and **v** are functions of *s*. These are calculated over three points along the lattice period, which are $s/S=0,0.125,0.25$. For a purely linear focusing case, it is clear from Fig. 1 and Eqs. (64), (65), and (66) that at $s/S=0$, corresponding to the longitudinal center of a quadrupole, the flow velocity is zero at all points, while the perturbation in the density is at its maximum. The opposite is true when $s/S=0.25$, corresponding to a point half way between a focusing and defocusing quadrupole, where the perturbation *n*_{2} vanishes at all points, and the flow velocity is at its maximum.

In this paper, explicit calculations will be made for a thermal equilibrium distribution. Some other equilibrium distribution functions are briefly discussed in Appendix B. The thermal distribution, based on the system of units used in this paper, is given by

where $n\u0302$ is the density at *R* = 0 and $T0$ represents the temperature of the beam. The expression for the zeroth order number density may be obtained by integrating over momentum space as given in Eq. (48), which gives

Inserting this into Eq. (51), we see that $\varphi 3$ needs to satisfy

This is a second order ordinary differential equation that needs to be solved numerically. This equation can be solved using a standard finite difference approach, but the solution would not converge as the tune depression ratio approaches 0.1, in which case specialized methods become necessary.^{22}

As an example, we use a case with the following parameters: $\Delta s/S=0.5,\u2009k2=1.04$ (or 60*°* phase advance per lattice period), $Q=3.81\xd710\u22124$ and $\eta =0.196$. It has been shown that beams with intense space charge usually degrade in beam quality when the phase advance per lattice period exceeds about 85°.^{46} Thus, a value of 60° is a realistic choice. Applying this set of dimensionless parameters for a 1 MeV beam of *K*^{+} ions, we obtain a *rms* beam size of 6.77 mm, current of 0.157 A, and a *rms* emittance of 18.61 mm-mr. These conditions compare reasonably well with the high current transport (HCT) experiment for heavy ion inertial fusion.^{47} They will vary for the same parameter set applied to a different species and energy. Yet the distribution function profiles in dimensionless coordinates will not change. To represent the dimensionless number density, we use the quantity $\pi Rrms2n(r)/Ntot$, where *n*(*r*) is the number density expressed up to the desired order, and *N _{tot}* is the total number of particles per unit length, obtained by integrating

*n*(

*r*) over “

*r*.” The radial distance itself is represented by $r/Rrms$. These units are identical to those used in Ref. 28 where is shown that for thermal equilibrium distributions, a beam under constant focusing is specified by a single dimensionless parameter. In the present system of units, this dimensionless parameter is given by

The tune depression ratio has then been shown to be related to *δ* by $\nu =1/(1+\delta )1/2$. Thus, the tune depression alone determines the zeroth order density profile of the beam distribution in thermal equilibrium in the given system of units.

Figure 2 shows the zeroth order and second order number density combined, in the same dimensionless units and at various values of *s*/*S*. Thus, in this case, the density would be $\pi Rrms2(n0(r)+n2(r,\theta ,s))/Ntot$, where we see that when perturbations are included, it becomes a function of *θ* and *s*. The lattice corresponds to that given in Fig. 1 without the sextupoles. The functions are displayed in Cartesian coordinates and based on the symmetry of the quadrupole fields, it is sufficient to calculate the density over one quadrant. At $s/S=0$, in Figure 2(a), the particles are at the longitudinal center of a horizontally focusing quadrupole ($\kappa 2>0$). At this point, the perturbation is at its maximum, as seen by Eq. (64) and the relative value of $\kappa 2II$ corresponding to this point is shown in Fig. 1. Figure 2(b) shows the density at $s/S=0.125$, which is away from the midpoint of the quadrupole, where we see that the perturbation becomes weaker. At $s/S=0.25$, shown in Figure 2(c), the particles are located at the midpoint between two adjacent quadrupoles. At this point, $\kappa 2II=0$, and so $n2=0$ which means the density is independent of *θ* and is identical to the unperturbed distribution obtained from the constant focusing term of the Hamiltonian. It should be noticed that there is a small area with a negative density in Figs. 2(a) and 2(b), which is unphysical. Additionally, there is a region of high density where the perturbation in density itself exceeds the zeroth order density. This occurs due to the steep drop in the number density near the edge of the beam, and because the unperturbed density is this region is very low. While a particle simulation would not produce these effects, this discrepancy may be disregarded as long as the it occurs over the edge, where the density is always very low. The general qualitative feature, corresponding to a depletion of charge near the edge along the defocusing direction, and an accumulation of charge at the edge along the focusing direction is still expected to be take place and could be determined quantitatively only through a high resolution particle simulation. The exact features of the edge may not be easily measurable in experiments.

The flow velocity function of the beam may be calculated from Eqs. (66) and (65), where the relative values of the quantity $\kappa I$ are shown in Fig. 1. In the presence of just quadrupoles, i.e., $\kappa 3=0$, the profile of the velocity distribution remains the same with only the magnitude changing as $\kappa I$ varies along *s*/*S*. We see that at $s/S=0$, the velocity is zero everywhere and this gradually increases to a peak when it exits the quadrupole. On entering a quadrupole of opposite polarity, it begins to decrease, reaching a zero at the center, and then reaching a negative peak upon exiting this quadrupole. Thus, we see that the velocity vanishes everywhere when the density perturbation reaches a maximum, and the velocity is at a maximum when the density perturbation vanishes. Figure 3 shows the velocity profile for $\kappa I>0$. When $\kappa I<0$ the flow direction is reversed. One can deduce from inspection that the direction of the flow shown in Fig. 3 is consistent with the density profile changing from Fig. 2(b) to Fig. 2(c). The densities and velocities at subsequent locations with the same intervals in *s*/*S* along the lattice would only be rotations of the corresponding distributions already shown here.

## VI. CALCULATIONS OF DISTRIBUTIONS WITH NONLINEAR FOCUSING

This section presents similar calculations for a nonlinear system that is equivalent to the one used in Sec. V. All the conditions for equivalence described in Sec. IV were made to be satisfied in these calculations, with the perveance and *rms* beam radius being within 1% of the required values. With the simultaneous presence of quadrupoles and sextupoles, it becomes necessary to calculate the beam properties over all quadrants, because the quadrupoles exhibit a 90° symmetry where as the sextupole exhibits a 60° symmetry. Figure 4 shows the density functions in the same dimensionless units used in Sec. V, for (a) $s/S=0$, (b) $s/S=0.125$, (c) $s/S=0.25$, and (d) $s/S=0.325$. As shown in Fig. 1, *s*/*S* corresponds to a point at the center along the length of a quadrupole, which also corresponds to a point midway between two sextupoles of opposite polarity.

At this point, the effect of the sextupoles cancels, i.e., $\kappa 3II=0$, and as a result the density profile takes a form that bears the symmetry of a quadrupole. The effect of the sextupole on the density shows up at $s/S=0.125$, where the quadrupole symmetry is broken. At *s* = 0.25, the beam is at the center of a sextupole and midway between two quadrpoles of opposite polarity. At this point, the density acquires the symmetry of the sextupole, with the effect of the adjacent quadrupoles canceling each other. At the next point, $s/S=0.325$, the effect of the vertically focussing quadrupole ($\kappa 2<0$) is clearly seen. The point $s/S=0.5$, not shown here would be the center of a vertically focusing quadrupole. The density function here would be a 90° rotation of the density function at $s/S=0$. Distributions at subsequent locations along the lattice with the same intervals in *s*/*S* may be obtained by rotating the corresponding distribution functions already shown here. Similar to the linear focusing case, these calculations also show some exaggeration in the accumulation and depletion of charge near the edge along the direction of focusing and defocusing respectively.

The velocity distribution is shown in Fig. 5 for the same points as in Fig. 4. These were calculated using Eqs. (66) and (65), where the relative values of the quantity $\kappa I$ are shown in Fig. 1. In these figures, we see that in the middle of a quadrupole, the profile has a 60° rotational symmetry, due to the sextupoles in the system. On the other hand, in the middle of a sextupole, we see a 90° symmetry in the velocity profile due to the quadrupoles in the system. A close examination of the direction of the velocity vectors clearly indicate that the density is moving from the one at the corresponding location to the one at the subsequent location shown in Fig. 4. Moreover, the velocity vectors at *s* = 0.625 indicate that the distribution is transforming from that shown in Fig. 4(d) to a 90° rotation of Fig. 4(a), where we see that the density is becoming symmetric about the origin along the horizontal axis.

Figure 6 shows the unperturbed densities of the linear and equivalent nonlinear cases discussed in this paper. We see that there is a slight hollowing of the density in the presence of nonlinear focusing, which is expected. We also see the steep drop in density as the beams passes the edge of the distribution for the nonlinear focusing case. This happens because the contribution from the nonlinear term increases at higher radial distances thus pushing the particles inward. In the linear focusing term, one would expect the particle oscillations to be largely linear in the region where the density is flat, while this will not be the case in the nonlinear focusing case, where neither the density is flat anywhere nor is the external focusing linear. The nonlinear oscillations of such a system are expected to help in mitigating instabilities driven from resonances arising from beam collective effects.

## VII. SUMMARY

In this paper, we have evaluated equilibrium distribution functions for beams with space charge effects in linear and nonlinear alternate gradient focusing systems. We first presented a method to systematically perform an averaging over fast oscillations for beams with space charge effects satisfying the Vlasov-Poisson system of equations. The Lie transform method was used, and it involved a perturbative canonical transformation. The formulation led to an expression for the averaged Hamiltonian along with higher order terms. Phase space distribution functions were obtained describing the unperturbed (or averaged) motion and also higher order terms which include the dynamics associated with the faster time scale, in this case the lattice period. Using the phase space distribution functions, expressions for the density and the flow velocity were obtained, that described the evolution of the beam along the direction of propagation. The functions describing the beam were always for a quasi thermal equilibrium case. The linear focusing channel would always have a zeroth order Hamiltonian that was azimuthally symmetric. In addition, it was shown that the elements in the nonlinear focusing channel can be chosen such that its Hamiltonian had azimuthal symmetry at the zeroth order. This condition of azimuthal symmetry brings the lattice closer to integrability which is known to lead to better confinement of the beam.

The equilibrium density and flow velocity functions were explicitly calculated for a certain class of parameters which involved highly intense space charge dominated beams. The features of these functions were discussed in detail. The calculations were performed along a series of point over a lattice period, showing how the beam properties evolve along this time scale. In some cases, the accumulation and depletion of charge due to focusing and defocusing were exaggerated, and this is a limitation of the analysis method. This effect takes place only near the edge of the beam and is more pronounced in highly space charge dominated beams which have a sharp edge.

The work in this paper can have various applications. For example, stability of a beam with respect to modes within the lattice period may be analyzed. One could also analyze resonances between single particle oscillations and such modes, where the beam distribution would be determined from the analysis of the paper, along with test particles whose trajectories may be studied. In numerical particle or fluid simulations, it is often useful to study the evolution of a beam that is perturbed about an equilibrium distribution. This helps eliminate various other effects such as phase space mixing which in turn helps isolate the specific phenomena being studied in the simulation. Thus, the analysis given in this paper would provide a guidance on setting up the initial conditions of such a simulation.

Nonlinear focusing has the potential to contribute to more stable motion of particles because the focusing force increases as a particle wanders away from the core of the beam. Nonlinear focusing also causes the collective mode oscillations of the beam to damp, such as oscillation of the *rms* radius. These phenomena help in mitigating beam loss due to halo formation and collective instabilities. Thus, it is important to understand the behaviour of beams in the presence of nonlinear focusing at a more fundamental level, which may not be possible from a simulation. This in turn can prove useful in studying the inclusion of nonlinear focusing effects while designing a beam transport system.

The KV distribution and its associated envelope equations, which is valid only for linear focusing systems, has been used extensively to analyze beams with space charge effects. The method used in this paper may serve as an alternate to KV type distributions for two reasons. First, since the formulation of this paper is applicable in more general forms of periodic focusing, it could include nonlinearities. Besides this, certain unphysical effects of the KV beam distribution may be eliminated when the methods of this paper are applied.

In the method described in this paper, the size of the perturbation depends upon the focusing strength of the elements as well as the periodicity of the lattice. Besides, the perturbation near the edge of the beam increases with increased space charge dominance, which would lead to inaccuracies, but may be disregarded since it occurs far away from the core of the distribution. Given that the size of the perturbation needs to be sufficiently small in order to ensure convergence in the series expansion, one needs to further investigate the parameter space in which this method may have limitations. However, we have shown that the method yields reasonable results with little unphysical effects for a realistic set of parameters. A more detailed study of beam evolution through simulations using these distributions as initial conditions under varying focusing strengths and tune depression ratios would also be very useful in this regard and may be the subject of a follow up study.

In conclusion, the procedure described in this paper can have applications for other Vlasov-Poisson systems beyond charged particle beams. In this regard, the treatment of the formulation has been kept as general as possible.

## ACKNOWLEDGMENTS

We are indebted to Steve Lund for his careful and critical evaluation of an early version of the paper, which greatly helped improve and enhance the content of the paper, bringing it to its present form. This work was funded by the U.S. Department of Energy under Grant No. DE-FG03–95ER40926.

### APPENDIX A: PROOF THAT TRUNCATION OF SERIES CONSERVES CHARGE

In Sec. II, it was stated that the perturbed distribution functions did not contribute to the total charge. This is true with every term in *f* higher than the zeroth order term *f*_{0}, and as a result the total charge remains an invariant irrespective of the order at which the series is truncated. To prove this, we first analyze a Poisson bracket expansion of *w _{n}* with

*f*

_{0}

We first prove that this expression must be zero when integrated over all phase space. Integrating the first term of the expansion by parts over $p$, we get

where we have used the condition that $f0\u21920$ as $p\u2192\u221e$.

Integrating the first term in the expansion once again by parts over $q$, and using the condition that $f1\u21920$ as $q\u2192\u221e$, we obtain

It is straight forward to extend this proof to the next order Poisson bracket operation on *f*_{0}. For example, using the fact that ${wn,f0}\u21920$ as $(q,p)\u2192\u221e$ and following the same procedure of integration by parts, it can be shown that

Continuing this extension to higher order operations, it is clear that the result would be the same for any order Poisson bracket operation of *w* on *f*_{0}.

From Eqs. (38) and (20), we know that $fn=Tnf0$. It can be seen from Eqs. (10) and (11) that *T _{n}* can be expressed as a combination of Poisson bracket expressions with respect to different terms of the function

*w*. Thus, based on the analysis in this Appendix, we have shown that

for all *n* > 0. This proves that each term *f _{n}*, for all

*n*> 0 does not contribute to total charge but only redistributes the zeroth order charge distribution in phase space. Consequently, truncating the series expansion in

*f*does not violate conservation of total charge.

### APPENDIX B: EXAMPLES OF BEAM EQUILIBRIA

While explicit calculations were made for a thermal equilibrium distribution, we now discuss examples of some well known forms of alternate equilibrium beam distributions for a constant focusing model, all of which satisfy Eq. (44). In our case, these distributions correspond to the zeroth order term, as shown by Eq. (38), which can then be used to obtain higher order terms using the results of Sec. III. Each of these distributions is expressed in terms of a beam temperature $T0$ and a characteristic beam size *R _{b}*. Depending upon the type of equilibrium distribution, the corresponding expressions for $T0$ and

*R*will be different.

_{b}The simplest example of a beam equilibrium is the KV distribution, which is valid only for a linear external focusing case. In this case, this would correspond to $\kappa 3=0$. The distribution for this is given by

where $T0$ is a constant that represents the temperature of the system expressed in units of the Hamiltonian. The corresponding number density is given by

where the characteristic beam radius *R _{b}* is given such that $V(Rb)=T0$. Using Eq. (51), it is easy to see that the potential $\varphi 3$ may be expressed as $\varphi 3=\u2212(qn0/4\u03f50)r2$. The KV equilibrium beam distribution can be generalized to any linear periodic focusing system with a stable lattice, where

*f*would be a

*delta*function of the Courant-Snyder invariants. This is the only known exact solution for a periodic focusing case, where the width of the beam varies with the lattice beta function. As indicated here, in the constant focusing case, the function is independent of

*s*. It should be noted that when the perturbation method of this paper is applied to an equivalent constant focusing KV distribution, it does not lead to the usual KV distribution corresponding to the full lattice. Additionally, while applying Eq. (64), the term $\u2202n0/\u2202r\u2192\u221e$ at the edge of the beam. However, since calculations need to be performed only within the region of nonzero

*n*

_{0}, one can set $\u2202n0/\u2202r=0$ for $r\u2265Rb$. Overall, the KV distribution has features that are unphysical, and are valid only for purely linear focusing systems. Thus, it is not an ideal case of study for the perturbation scheme of this paper, which is better suited for more physical distributions, and a focusing force more general than a purely linear one.

The rest of the distribution functions mentioned here have exact solutions only for a constant focusing system which may be linear or nonlinear. These include the well known water bag equilibrium distribution, the parabolic distribution,^{22} and the truncated thermal distribution.^{48}

In the system of units used in this paper, the expression for the water bag distribution takes the form

where Θ is the Heaviside step function defined by $\Theta (x)=1$ for $0\u2264x<1$, and $\Theta (x)=0$ for *x* > 1. The corresponding zeroth order number density may be obtained by integrating this in momentum space, which is

As in the KV distribution case, *R _{b}* is such that $V(Rb)=T0$. Substituting this into Eq. (75), the electrostatic potential $\varphi 3$ may be determined. For $0\u2264r<Rb$, we have

For a purely linear focusing case, transforming from $\varphi 3$ to *V*(*r*) as the dependent variable, the above equation leads to a linear inhomogeneous Bessel's equation. Analytic expressions for *n*_{0} have been obtained for this case (see, for example, Ref. 1). The more general case, where $\kappa 3\u22600$ will still need to be solved numerically.

The parabolic function is given as follows:

Integrating over the momentum space as before, we get

with *R _{b}* satisfying $V(Rb)=T0$.

The truncated thermal distribution combines aspects of a thermal distribution and a water-bag distribution and is given as follows:

Integrating over the momentum space as before, we get

In this case *R _{b}* is defined such that $V(Rb)=\chi T0$. The parameter

*χ*specifies the point at which the distribution gets truncated.

Overall, the equilibrium distribution functions given in this section are not as well suited as the thermal distribution for the analysis of this paper. This is because they are all described by either a *delta* function as in the case of a KV distribution, or the Heaviside step function. The analysis of the paper requires the functions to be well behaved. On the other hand, for all the cases described in this Appendix, the errors occur near the edge of the beam and may be considered unimportant based on the nature of the study.