Generalized parametrization methods for centroid and envelope dynamics of charged particle beams in coupled lattices

Although there is continued interest in producing an elegant description of coupled linear beam dynamics (see, for example, Refs. 1–4), no single method is yet in general use among the beam physics community in the same way as the well-known Courant-Snyder (CS) theory is for uncoupled linear systems.⁵ Based on a time-dependent canonical transformation, Qin and Davidson proposed a parametrization method that extends the original CS theory for one degree of freedom to higher dimensions and involves remarkably similar mathematical expressions and physical meanings.^6,7 Their parametrization method has been applied to general linear focusing lattices composed of quadrupole, skew-quadrupole, and solenoidal magnets, including the effects of energy variations.^8,9 Furthermore, the Kapchinskij-Vladimirskij (KV) distribution¹⁰ has been generalized to coupled transverse dynamics, based on the generalized CS invariant identified through the Qin-Davidson parametrization.^11–13 

In both the original CS theory and the generalized CS theory introduced by Qin and Davidson, it is assumed that any centroid offset from the reference orbit is negligible, or the coordinates are redefined with respect to the offset centroid.¹⁴ Therefore, most beam dynamics studies have only considered the envelope dynamics. However, in order to design certain sophisticated beam manipulation devices, such as the oscillating wobbler^15,16 and injection/extraction devices,^14,17 a precise description of the beam centroid dynamics is important. Moreover, for the analysis of the centroid oscillations resulting from the imperfections in bending magnets¹⁷ and certain beam instabilities (e.g., two-stream electron cloud instability and beam-beam interactions),^15,16 effective physical models of the centroid motions are required.

In uncoupled lattice systems, it is straightforward to show that the envelope dynamics and centroid dynamics are decoupled when the image charge effects are negligible.^15,17 To the authors' knowledge, there has been little effort to incorporate beam centroid motions into general coupled beam dynamics. On the other hand, from the viewpoint of wave packet motions, the centroid dynamics have been investigated for quantum-mechanical systems by introducing time-dependent linear terms in quadratic Hamiltonians.^18,19 In Refs. 18 and 19, it was found that the centroid motion of any wave packet becomes decoupled from that of any moments relative to the centroid, regardless of what time dependence is present in the quadratic Hamiltonian.

Hence, in this paper we revisit the Qin-Davidson parametrization method and extend it to include the beam centroid motions. In Secs. II and III, the linear terms in the quadratic Hamiltonian for the coupled beam dynamics are introduced and treated by means of a special time-dependent canonical transformation. The generalized CS theory is reviewed in Sec. IV, together with the redefinition of the CS invariant in the presence of a centroid motion. In Sec. V, we derive the governing equations for the beam centroid evolutions in response to coupled lattices and time-dependent driving forces. Our conclusions and directions for future work are summarized in Sec. VI.

As a general remark, the effects of the linear terms in quadratic Hamiltonians have not been fully explored in linear beam dynamics,¹⁵ in particular for coupled lattice systems. To investigate these effects systematically, let us consider a time-dependent quadratic Hamiltonian in a general linear focusing lattice with n degrees of freedom

Here, z = (x₁,…,x_n, p₁,…,p_n)^T denotes the canonical coordinates; s is the path length, which plays the role of a time-like variable; and A(s) is a 2n × 2n time-dependent symmetric matrix. We also note that f(s) = (f₁,…,f_2n)^T is a 2n × 1 time-dependent column vector, and the term f(s)^Tz represents terms that are linear in positions and momenta. For example, a term that is linear in position appears in the Hamiltonian when there is a spatially uniform driving force, such as a uniform electric force.¹⁹ The superscript “T” denotes the transpose operation of a matrix. The canonical momenta are normalized by a fixed reference momentum p₀ = γ₀m_bβ₀c, where m_b is the rest mass of a beam particle, c is the speed of light in vacuo, and β₀ and $γ0=(1−β02)−1/2$ are the relativistic factors.

The Hamiltonian equation of motion for z is given by^9,20

where J is the 2n × 2n unit symplectic matrix

and I is the n × n unit matrix. Making use of index notation, we obtain that

Here, we use the fact that the matrix A is symmetric, that is, A_jk = A_kj. The prime denotes a derivative with respect to s. If we rewrite Eq. (4) in matrix notation, we obtain the expression

We note that the time-dependent function g(s) does not appear in Eqs. (4) and (5). The Hamiltonian equations of motion are not affected by the addition of an arbitrary time-dependent function g(s) in the Hamiltonian, and thus, we set g(s) = 0 for the remainder of this discussion.²¹ 

For the linear coupled dynamics of a charged particle relative to a reference orbit, the matrix A can be expressed in its most general form as

where κ(s), R(s), and m⁻¹(s) are time-dependent n × n matrices. The matrices κ(s) and m⁻¹(s) are also symmetric. The Hamiltonian equations of motion can then be expressed as

Here, $x=(x1,…,xn)T, p=(p1,…,pn)T, fx=(f1,…,fn)T$⁠, and $fp=(fn+1,…,f2n)T$ are n × 1 column vectors.

There are several cases in practical accelerators in which the Hamiltonian describing the charged particle motion contains linear terms in the canonical coordinates.

The wobbler is a set of electrically (or magnetically) biased plates driven by an RF source.¹⁵ The wobbler can actively control the centroid dynamics by means of oscillating deflecting forces, which can be included in the Hamiltonian as linear terms. In ion-beam-driven high energy density physics and heavy ion fusion, the wobbler system has been considered as a possible beam smoothing technique to achieve a uniform illumination on a target.¹⁶ In certain proton and carbon therapy facilities, wobbler systems have been employed to produce a desired dose distribution on patients.²² 

Regarding particle motion in a dipole magnet, a particle with the reference momentum will follow the desired reference orbit.⁴ When the particle has a non-zero energy deviation δ, it will deviate from the reference orbit. If there is no transverse-longitudinal coupling element (e.g., TM₁₁₀ deflecting cavity), then the energy deviation can be treated independently of the transverse motion and gives rise to a term linearly proportional to the transverse coordinate in the Hamiltonian. This linear term produces a transverse deflecting force in the equation of motion, which is proportional to the energy deviation δ.

The Hamiltonian of a particle in an RF cavity is often obtained by averaging the Hamiltonian over the cavity length.⁴ The resulting averaged Hamiltonian contains a term that is linearly proportional to the longitudinal coordinate. This linear term represents acceleration, where the energy gain is proportional to the cavity voltage and also depends on the RF phase.⁴ In the equation of motion, the resulting force term becomes spatially uniform. If the RF phase is set in such a way that the reference particle passes the cavity at a zero-crossing, then there will be no acceleration.

Normally, nonlinear perturbations are assumed to be very small, and so, we can insert the unperturbed solution x₀ and p₀ into the expression for a nonlinear perturbation, making it a time-dependent driving term.²³ Indeed, any driving force term in the equation of motion, which is independent of x and p, can be treated using linear terms in the Hamiltonian. Often, such perturbations arise either externally from the errors in magnetic field strengths and alignments or internally from the nonlinear space-charge forces.

We introduce a time-dependent canonical transformation^8,9,21,24

such that in the new coordinates $z¯$⁠, the Hamiltonian (1) may be transformed into

where $A¯(s)$ is a targeted symmetric matrix. Here, r is a real column vector. Because S(s) is a symplectic matrix, it is required that

Similar to Eq. (5), the Hamiltonian equation of motion for $z¯$ is given by

At the same time, $z¯′$ can be obtained by taking a direct time-derivative of Eq. (9), which yields that

By equating Eqs. (12) and (13), we obtain that

Equation (14) should be valid for any z. Therefore, it is required that

Meanwhile, the remaining terms in Eq. (14) can be expressed as

Because a symplectic matrix S must have an inverse matrix S⁻¹, Eq. (16) can be reduced to

We note that the vector r, which removes the linear terms in the original Hamiltonian, satisfies the same equation of motion for a single charged particle [see Eq. (5)].

In addition, one can check that the linear transformation $z¯=S(s)(z−r)$ is canonical when the matrix S is symplectic by evaluating the square matrix Poisson bracket^20,21

Here, the symplectic condition S^TJS = J is equivalent to Eq. (11).

As a first step, we now apply the coordinate transformation introduced in Sec. III, $z¯=S(s)(z−r)$⁠, such that in the new coordinate system the Hamiltonian takes the following form:^8,9

where μ(s) is an unknown n × n symmetric matrix to be determined together with the symplectic matrix S(s). Because both A(s) and $A¯(s)$ are given as 2 × 2 block matrices, we also let

We now can split the Eq. (15) into four matrix differential equations

From the symplectic condition SJS^T = J, we have three more constraints

In the four matrix differential equations (21)–(24), there are five n × n matrices to be determined. Therefore, this is an overdetermined system. Among the infinite ways to remove the redundant freedom, we choose S₂ ≡ 0, following the notational convention used in the original CS theory. Furthermore, we define $W=S4T$⁠, in order to clearly indicate that the physical meaning of $S4T$ will turn out to be the matrix version of the beam envelope. From the symplectic condition S₁W = I, we obtain that S₁ = W⁻¹. Then, the matrix differential equations (21)–(24) can be expressed as

From Eq. (29), we obtain that

After defining $V=−S3T$⁠, we can express Eqs. (30) and (31) as

where the variable V can be considered to be the matrix associated with the envelope momentum. Unlike in the previous generalized CS theory in Refs. 8 and 9, we express the matrix envelope equation in terms of two first-order equations. In this manner, we observe that Eqs. (33) and (34) have a similar Hamiltonian structure to the single particle equations of motion,²⁵ except for the term (W^TmWW^T)⁻¹ [see Eqs. (7) and (8) for comparison].

So far, we have not explicitly imposed the symplectic condition in Eq. (26), which can be expressed as

By the basic identity for the derivative of an inverse $d(W−1)/ds=−W−1(dW/ds)W−1$⁠, it is straightforward to show that Eq. (28) is indeed equivalent to Eq. (35). This condition is trivial in the original CS theory, where V and W are scalar and commutative. On the other hand, in the generalized CS theory particular initial conditions for V and W are required when solving the matrix envelope Eqs. (33) and (34) in order to ensure that (V^TW − W^TV) = 0 for all s. By taking $[Eq. (33)]TV−VT[Eq. (33)]−[Eq. (34)]TW+WT[Eq. (34)]$⁠, it can be shown that

Therefore, if we set the initial values for W and V at s = 0 in such a way that (V^TW − W^TV)₀ = 0, then the symplectic condition is satisfied for all s > 0, and it is not necessary to additionally solve Eq. (28). Except for this symplectic condition, the initial conditions for W and V can be arbitrary. Once W and V are solved from Eqs. (33) and (34) with an appropriate initial condition satisfying Eq. (35), we can determine the symplectic matrix S(s) and its inverse in the following forms:

In calculating S⁻¹ in Eq. (38), we have used the symplectic condition (V^TW − W^TV) = 0.

The second step is to apply another canonical transformation $z¯¯=P(s)z¯$ to transform $H¯$ into a vanishing Hamiltonian, that is, $H¯¯=0$ for all s. Equation (15) can be expressed as

Following the same procedure as in the first step, the differential equation for P is split into four matrix differential equations

Because Eqs. (40) and (41) have the same form as Eqs. (42) and (43), we can set P₄ = P₁ and P₃ = –P₂. Of course, there are many other possible choices, such as P₃ = P₁ and P₄ = P₂. However, none of these are consistent with the symplectic condition PJP^T = J except for the solutions with P₄ = P₁ and P₃ = –P₂. Then, the symplectic condition can be expressed as

We should solve Eqs. (40) and (41) together with the symplectic conditions (44) and (45). By taking $[Eq. (40)]P2T−[Eq. (41)]P1T−P2[Eq. (40)]T+P1[Eq. (41)]T$⁠, we note that

Furthermore, by taking $[Eq. (40)]P1T+[Eq. (41)]P2T+P1[Eq. (40)]T+P2[Eq. (41)]T$⁠, we note that

Therefore, if we set the initial values for P₁ and P₂ at s = 0 in such a way that $(P1P2T−P2P1T)0=0$ and $(P1P1T+P2P2T)0=I$⁠, then the symplectic conditions are satisfied for all s > 0. The symplectic conditions (44) and (45) are also equivalent to

which indicates that P is not only a symplectic matrix but also a rotation matrix. Therefore, P corresponds to a symplectic rotation in 2n-dimensional phase space. Based on the analogy with the 2D rotation matrix, we rename P₁ = C_o and P₂ = −S_i. Once C_o(s) and S_i(s) are solved from Eqs. (40) and (41) with appropriate initial conditions satisfying Eqs. (44) and (45), we can determine the symplectic rotation matrix P(s) and its inverse in the following forms:

Here, C_o and S_i satisfy

where the term μ = (W^TmW)⁻¹ represents the phase advance rate.

By combining the two symplectic transformations from the first and second steps, we obtain

Because $H¯¯=0$⁠, the dynamics of $z¯¯$ is trivial, that is, $z¯¯=z¯¯0=const$⁠. From $z¯¯=P(s)S(s)(z−r)$ and $z¯¯0=P0S0(z0−r0)$⁠, we have a linear map

where subscript “0” denotes values at s = 0. Because P is a symplectic rotation matrix, P⁻¹ = P^T and P₀ can be taken to be a unit matrix without loss of generality. The transfer matrix M is defined by

Similar to the original CS theory, the linear motion with respect to centroid becomes a time-dependent rotation after the normalization of the phase-space coordinates.²⁶ The first matrix from the left in M(s) is a back-transformation to the original coordinate system.⁹ 

Because $z¯¯$ is a constant column vector, we can construct an invariant of the dynamics that is quadratic in the phase-space coordinates by introducing a constant 2n × 2n positive-definite matrix ξ

The quantity I_ξ is the generalized version of the Courant-Snyder invariant, which can be defined even in the presence of a centroid motion. The matrix ξ acquires a meaning associated with emittance when the beam distribution is defined in terms of the generalized CS invariant I_ξ.¹³ 

In principle, the vector r can evolve according to Eq. (17) with an arbitrary choice of initial conditions r₀ at s = 0. However, in this case the physical interpretation of the second moments of the beam distribution in terms of the CS invariant I_ξ becomes non-trivial. To examine the physical meaning of the vector r, we consider the following form of the beam distribution function:

which is a solution of the Vlasov equation (i.e., df_b/ds = 0, because I_ξ is a constant of motion). We note that the contours of the phase-space density are determined by a single invariant I_ξ, which is the usual assumption in beam physics. This beam distribution can be characterized by a multivariate Gaussian distribution

where

is a covariant matrix, and $〈⋯〉$ denotes the statistical average over the distribution. Then, we can set

By comparing Eqs. (58) and (62), we readily obtain that

The vector r acquires a clear physical meaning when it is set to be the centroid of the beam. Then, from Eq. (17) the equation describing the evolution of the beam centroid can be expressed as

The above equation can be derived more rigorously by using the statistically-averaged rate equations. The statistical average of a phase function χ(x, p, s) over the 2n-dimensional phase space (x, p) is denoted by $〈χ〉(s)$ and is defined by²⁷ 

Here, the beam distribution function f_b(x, p, s) is normalized as $1=∫dxdpfb$ and evolves according to the Vlasov equation²⁷ 

For notational simplicity, we introduce the spatial gradient $∇x=(∂/∂x1,…,∂/∂xn)T$ and momentum gradient $∇p=(∂/∂p1,…,∂/∂pn)T$ and assume that the dot product of two column vectors a and b is equivalent to the matrix multiplication of the row vector representation of a and the column vector representation of b, that is, $a·b=aTb$⁠.²⁸ From the definition of the statistical average in Eq. (66), it follows that

After multiplying the Vlasov equation (67) by χ and operating with $∫dxdp$⁠, we obtain that

By applying integration by parts with respect to x with the appropriate boundary conditions for f_b at the beam boundary,²⁷ the third term in Eq. (69) can be expressed as

Similarly, the fourth term in Eq. (69) can be expressed as

Substituting Eqs. (70) and (71) into Eq. (69) then gives the general rate equation for the coupled lattice systems

Here, the terms involving the trace of the matrix R are cancelled out, that is, $−〈χTr(RT)〉+〈χTr(R)〉=0$⁠. Substituting χ = x₁,…,x_n into Eq. (72) gives us that

Similarly, for χ = p₁,…,p_n, we obtain from Eq. (72) that

Equations (73) and (74) describe the evolution of the beam centroid in phase space $〈z〉=(〈x〉T,〈p〉T)T$ in response to the coupled lattices and uniform driving forces in the phase space. We also note that Eqs. (73) and (74) can be combined to give

which is indeed the same as Eq. (65). Therefore, the centroid motion follows exactly the same equations of motion as the individual beam particle in phase space, even for general coupled lattices.

For charged particle beam dynamics in uncoupled lattices, a fundamental aspect is that the centroid dynamics and envelope dynamics are decoupled in the absence of image charge effects. In this work, we have extended this desirable decoupling feature of the centroid and envelope dynamics to the case of general linear coupled systems. It has been found that the envelope matrices (W and V), the phase advance matrix (P), and the beam matrix (Σ) are defined relative to the beam centroid, and they evolve independently of the linear terms in the quadratic Hamiltonian. We clarified further that all of the formulations established by Qin-Davidson's generalized CS theory are applicable, even in the presence of a centroid motion. The generalized CS theory provides an effective framework for studying beam dynamics in complex coupled lattices, with the freedom to control centroid motions. Hence, it would help to discover a more optimized lattice design and beam manipulation scheme in larger parameter spaces. The addition of space-charge effects to the centroid and envelope motions in coupled beam dynamics is already being investigated, and that work will be published elsewhere, together with numerical examples.

This research was supported by the National Research Foundation of Korea (Grant Nos. NRF-2015R1D1A1A01061074 and NRF-2017M1A7A1A02016413). This work was also supported by the U.S. Department of Energy Grant No. DE-AC02-09CH11466.