Analytical methods for describing charged particle dynamics in general focusing lattices using generalized Courant-Snyder theory

The dynamics of charged particles in general linear focusing lattices with quadrupole, skew-quadrupole, dipole, and solenoidal components, as well as torsion of the fiducial orbit and variation of beam energy is parameterized using a generalized Courant-Snyder (CS) theory, which extends the original CS theory for one degree of freedom to higher dimensions. The envelope function is generalized into an envelope matrix, and the phase advance is generalized into a 4D symplectic rotation, or an U(2) element. The 1D envelope equation, also known as the Ermakov-Milne-Pinney equation in quantum mechanics, is generalized to an envelope matrix equation in higher dimensions. Other components of the original CS theory, such as the transfer matrix, Twiss functions, and CS invariant (also known as the Lewis invariant) all have their counterparts, with remarkably similar expressions, in the generalized theory. The gauge group structure of the generalized theory is analyzed. By fixing the gauge freedom with a desired symmetry, the generalized CS parameterization assumes the form of the modified Iwasawa decomposition, whose importance in phase space optics and phase space quantum mechanics has been recently realized. This gauge fixing also symmetrizes the generalized envelope equation and express the theory using only the generalized Twiss function \beta.The generalized phase advance completely determines the spectral and structural stability properties of a general focusing lattice. For structural stability, the generalized CS theory enables application of the Krein-Moser theory to greatly simplify the stability analysis. The generalized CS theory provides an effective tool to study coupled dynamics and to discover more optimized lattice design in the larger parameter space of general focusing lattices.


I. INTRODUCTION
In accelerators and storage rings, charged particles are confined transversely by electromagnetic focusing lattices. Many different kinds of focusing lattices have been successfully designed and implemented. The fundamental theoretical tool in designing an uncoupled quadrupole lattice is the Courant-Snyder (CS) theory [1], which can be summarized as follows. For a given set of focusing lattices in the x-and y-directions κ x ðtÞ and κ y ðtÞ, a particle's dynamics is governed by the oscillation equation q þ κ q ðtÞq ¼ 0; (1) where q represents one of the transverse coordinates, either x or y. The solution of Eq. (1) can be expressed as a symplectic linear map MðtÞ that advances the phase space coordinates In CS theory, the linear map MðtÞ is given as where αðtÞ and βðtÞ are two of the so-called Twiss parameters, and ϕðtÞ is the phase advance. They are defined by the envelope function wðtÞ as ϕðtÞ ¼ Z t 0 dt βðtÞ ; (6) and the envelope function wðtÞ is determined by the envelope equation̈w þ κ q ðtÞw ¼ w −3 : In Eq. (3), q 0 ¼ qðt ¼ 0Þ, _ q 0 ¼ _ qðt ¼ 0Þ, β 0 ¼ βðt ¼ 0Þ, and α 0 ¼ αðt ¼ 0Þ are the initial conditions at t ¼ 0.
Associated with the dynamics of Eq. (1), there exists a constant of the motion, I CS , known as the Courant-Synder invariant Here γðtÞ is the third Twiss parameter. It turns out that the transfer matrix MðtÞ can be decomposed into the elegant form [2] MðtÞ which seems to indicate a certain structure for MðtÞ. The CS theory can be viewed as a parametrization method of the time-dependent 2 × 2 symplectic matrix MðtÞ for a standard uncoupled lattice. Not surprisingly, there exist other parametrization schemes mathematically. Why is the CS parametrization preferable? This is because it describes the physics of charged particle dynamics. The main components of the CS theory, i.e., the phase advance, the envelope equation, the transfer matrix, and the CS invariant are physical quantities describing the dynamics of the particles. For example, the CS invariant defines the emittance in phase space, and the envelope function describes the transverse dimensions in configuration space. This theoretical framework also makes it possible to investigate collective effects associated with high-intensity beams, such as in the construction of the Kapchinskij-Vladimirskij distribution [3][4][5].
However, the CS theory can only be applied to the x-or y-dynamics separately for the ideal case of uncoupled quadrupole focusing lattices. In realistic accelerators, there exist bending magnets, torsion of the design orbit (fiducial orbit), and skew-quadrupole components, which are introduced intentionally or by misalignment [6,7]. Solenoidal magnets are also used in certain applications [8]. When these additional components are included, the transverse dynamics in the x-and y-directions are coupled, and the focusing force depends on the transverse momentum as well. In this most general case, the transfer matrix MðtÞ is a time-dependent 4 × 4 symplectic matrix, which has 10 time-dependent parameters and admits many different schemes for parametrization. The first set of parametrization schemes for MðtÞ was developed by Teng and Edwards [9][10][11] and Ripken [12][13][14], some of which have been adopted in lattice design and particle tracking codes, such as the MAD code [15,16]. A class of different parametrizations by directly generalizing the Twiss parameters to higher dimensions has also been developed by Dattoli et al. [17][18][19]. However, in contrast to the original CS theory, these parametrization schemes are designed from mathematical considerations, and fail to connect with physical parameters of the beam. The elegant and muchneeded connection with the physics of beam dynamics in the original CS theory for one degree of freedom is not transparent in these parametrization schemes. This is probably why there is no de facto standard yet adopted by the accelerator community. Another main reason is that, for most present-day accelerators and rings, the transverse dynamics are so nearly decoupled that perturbative treatments often work satisfactorily. Even for lattices with strong coupling, elementary methods can be used to analyze the dynamics, even though the calculation often becomes rather involved and requires diligence and patience.
In a recent Letter [20], we reported on the development of a generalized CS theory for focusing lattices with the most general form in Eq. (11), including bending magnets, torsion of the design orbit, and solenoidal magnets, in addition to quadrupole and skew-quadrupole components. In this generalized theory, the physics elements of the original CS theory, i.e., the phase advance, the envelope equation, the transfer matrix, and the CS invariant are all generalized to the 2D coupled case with identical structure. This new development also generalizes our previous results for coupled dynamics including only a skew-quadrupole lattice component [21][22][23][24]. In this paper, we give a detailed derivation of the generalized CS theory reported in Ref. [20], describe the theoretical structure of the theory in terms of gauge freedoms and group decomposition, and demonstrate the application of the theory in a stability analysis.

II. THEORETICAL MODEL AND SUMMARY OF RESULTS
In this section, we outline the theoretical methods used and summarize the main results obtained in this paper. As discussed in Sec. I, when realistic components such as skew-quadrupoles, bending magnets, torsion of the design orbit, solenoidal magnets are included, in addition to the standard quadrupole components, the transverse dynamics in the x-and y-directions are coupled, and the focusing force depends on the transverse momentum. In this case, the linear dynamics of a charged particle relative to the fiducial orbit are governed by a general time-dependent Hamiltonian [25]  Here, the components of z ¼ ðx; y; p x ; p y Þ T are the transverse phase space coordinates, and κðtÞ, RðtÞ, and m −1 ðtÞ are time-dependent 2 × 2 matrices. The matrices A, κðtÞ, and m −1 ðtÞ are also symmetric. In this most general Hamiltonian, the skew-quadrupole and dipole components are included in the off-diagonal terms of the κðtÞ matrix, and the solenoidal component and the torsion of the fiducial orbit are included in the RðtÞ matrix. There are several different methods to include the effect of torsion, which were reviewed by Hoffstaetter [26]. Typically, Frenet-Serret coordinates along the fiducial orbit are used. When the fiducial orbit is straight, the Frenet-Serret coordinates are not uniquely defined. In this case, we can choose any particular set of Frenet-Serret coordinates in the straight section, as long it is smoothly connected to those in the curved sections. The variation of beam energy along the fiducial orbit is reflected in the mass matrix m −1 ðtÞ, which is allowed to be any real symmetric matrix for complete generality. The transfer matrix MðtÞ corresponding to H is a time-dependent 4 × 4 symplectic matrix, which has 10 time-dependent parameters. Our goal is to develop a generalized Courant-Snyder parametrization method for MðtÞ, which has the same elegant structure and direct connection to beam dynamics as the original Courant-Snyder theory for one degree of freedom. We will use a time-dependent symplectic transformation technique [21,24,27] to analyze the charged particle dynamics governed by the Hamiltonian given in Eq. (11). This technique is described in Sec. III. The concept of scalar envelope function is generalized to a 2 × 2 envelope matrix, and the envelope equation in 2 × 2 matrix form is developed [see Eq. (36)] in Sec. IV. In the original CS theory, the envelope equation (7) is one dimensional and plays a central role. It also has been discovered or rediscovered many times [28][29][30][31][32] in other branches of physics. In quantum physics, it is known as the Ermakov-Milne-Pinney equation [28][29][30], which has been utilized to study 1D time-dependent quantum systems [33,34] and associated nonadiabatic Berry phases [35]. A brief account of the history of the 1D envelope equation can be found in Ref. [36]. We expect the generalization of the envelope equation to higher dimensions for the most general Hamiltonian to have applications in areas other than beam physics as well. The 1D CS invariant given by Eq. (8), also known as the Lewis invariant [31,32] in quantum physics, is generalized to higher dimensions in Eqs. (44) and (45).
Also in Sec. IV, the 1D phase advance is generalized to a time-dependent matrix P, which belongs to the symplectic rotation group Spð4Þ⋂SOð4Þ ¼ Uð2Þ. Here, Spð4Þ, SOð4Þ, and Uð2Þ denote the groups of 4 × 4 symplectic matrices, 4 × 4 rotation matrices, and 2 × 2 unitary matrices, respectively. For dynamics with one degree of freedom, the phase advance is naturally an angle [an element of SOð2Þ] in the 2D phase space. For dynamics with two degrees of freedom, the phase space is 4D, and it is tempting to represent the phase advance by two angles. This is what has been adopted in previous parametrization schemes. From the viewpoint of theoretical physics and geometry, however, it is more natural to represent the phase advance for dynamics with two degrees of freedom by a 4D rotation [an element of SOð4Þ], which is not equivalent to a pair of 2D rotations. Because of the symplectic nature of the Hamiltonian dynamics, generalized phase advance in higher dimensions thus belongs to the symplectic rotation group. Of course, one can adopt different views on this. In the normal-form analysis of accelerator rings, the 4D transfer matrix is reduced to a 2D rotation after block diagonalization. In a sense, we can compare coupled betatron motion to the Dirac equation. In a fully quantum mechanical limit, the only correct approach is to treat electrons and positrons as inextricably coupled; only the 4D approach is permissible in this limit. But in accelerator physics, we routinely treat electrons and positrons as completely separate entities, and two 2D descriptions are adopted without much hesitation.
The generalized decomposition for the symplectic map MðtÞ is given by Eq. (43), which has exactly the same structure as the original 1D CS theory given by Eq. (10). In addition to its aesthetic elegance, the generalized CS theory provides an effective tool to describe the beam dynamics governed by the most general Hamiltonian. The 2 × 2 envelope matrix w defines the transverse dimension of the beam, and the generalized CS invariant defines the emittance. These components of the generalized CS theory are derived in detail in Sec. IV. For the present application to beam transverse dynamics, there are two degrees of freedom. But the theory developed is valid for any degree of freedom. For a system with n degrees of freedom, the time-dependent matrix AðtÞ specifying the Hamiltonian in Eq. (11) will be 2n × 2n, the envelope matrix will be n × n, and the phase advance will belong to Spð2nÞ⋂ SOð2nÞ ¼ UðnÞ.
In Sec. V, we investigate the group structure of the generalized CS theory, which is built on the decomposition of the time-dependent symplectic coordinate transformation G in the form of Eq. (42). There exists a gauge freedom in this decomposition specified by a 2D rotation element c ∈ SOð2Þ for every t. The transfer map MðtÞ is independent of this gauge. By fixing the gauge freedom with a desired symmetry, the decomposition of G as PS assumes the form of the modified Iwasawa decomposition (or pre-Iwasawa decomposition), whose importance in phase space optics [37,38] and phase space quantum mechanics [39] has been recently realized. This specific gauge fixing also symmetrizes the generalized envelope equation and expresses the theory using only the generalized Twiss function β. For a symplectic matrix, the modified Iwasawa decomposition is equivalent to the well-known Iwasawa decomposition for a semisimple Lie group [40]. However, the unique feature of the theory described here is that the decomposition is constructed as a function of time, and from the viewpoint of dynamics using the generalized envelope equation. Nevertheless, it is a pleasant surprise to find the deep connection between the original CS theory for charged particle dynamics [1] and the Iwasawa decomposition for Lie groups [40], two theoretical formalisms developed concurrently. This connection also demonstrates that beam dynamics, phase space optics, and quantum dynamics have a similar theoretical structure at the fundamental level. In order to satisfy the symmetry requirement of the modified Iwasawa decomposition, the gauge freedom needs to be selected locally as a function of time, which is characteristic of gauge theories in theoretical physics. This procedure also results in a symmetrized envelope equation in terms of the generalized Twiss parameter β, which is a symmetric, positive-definite matrix. The beam dimensions and emittance can be expressed using the β matrix only.
We show in Sec. VI how the generalized CS theory can be used to analyze the stability of charge particle dynamics in realistic accelerators with quadrupole, skew-quadrupole, dipole, and solenoidal components, as well as torsion of the fiducial orbit and variation of beam energy. It turns out that the generalized phase advance as a symplectic rotation completely determines the spectral and structural stability properties of the general lattice after a matched solution of the envelope equation is found. For structural stability, the generalized CS theory enables us to apply the Krein-Moser theory [41][42][43][44][45] to greatly simplify the stability analysis. This general result includes the well-known stability criterion for sum/difference resonances for uncoupled quadrupole lattices as a special case.

III. METHOD OF TIME-DEPENDENT CANONICAL TRANSFORMATION
We will construct the generalized Courant-Snyder theory for the general focusing lattice given by Eq. (11) using a time-dependent canonical coordinate transformation technique. Let us consider a linear, time-dependent Hamiltonian system with n degrees of freedom Here, AðtÞ is a 2n × 2n time-dependent, symmetric matrix. The Hamiltonian in Eq. (11) has this form with n ¼ 2.
The basic idea is to introduce a time-dependent linear canonical transformation [27] z ¼ SðtÞz; (13) such that in the new coordinatesz, the transformed Hamiltonian has the desired form whereĀðtÞ is a targeted symmetric matrix. Because the transformation (13) is canonical, it requires that Here, J the 2n × 2n unit symplectic matrix of order 2n, and I is the n × n unit matrix. Equation (15) implies that S is a symplectic matrix. In addition, it needs to satisfy a differential equation, which can be derived as follows.
Hamilton's equation for z is given by Using index notation, Eq. (17) becomes Switching back to matrix notation, Eq. (18) can be expressed as Similarly, Meanwhile, _ z can be directly calculated from Eq. (13) by taking a time-derivative, which gives Combining Eqs. (20) and (21) gives the differential equation for S The remarkable feature of the canonical transformation S is that it is always symplectic if S is symplectic at t ¼ 0. This assertion can be proved by two methods. For the first proof, we follow Leach [27] and consider the dynamics of the matrix K ¼ SJS T , Equation (23) has a fixed point at for all t, and S is symplectic for all t. A more geometric proof can be given from the viewpoint of the flow of S. Because A is symmetric, we have JJĀ −Ā T JJ ¼ 0, which indicates that JĀ belongs to the Lie algebra spð2nÞ. We now show that if S is symplectic at a given t, then JĀS belongs to the tangent space of Spð2nÞ at S, i.e., JĀS ∈ T S SPð2nÞ. Let us examine the Lie group right action: S: a↦aS for any a in Spð2nÞ, and the associated tangent map It is evident that JĀS is the image of the Lie algebra element JĀ under the tangential map T S . This means that JĀS is a vector tangential to Spð2nÞ at S, if S is on Spð2nÞ.
By the same argument SJA ∈ T S SPð2nÞ as well. Thus, the right-hand side of Eq. (22) is a vector in Spð2nÞ, and the S dynamics will stay in Spð2nÞ. We can always choose initial conditions such that S is symplectic at t ¼ 0, and this will guarantee that the time-dependent transformation specified by Eq. (22) is symplectic for all t.

IV. GENERALIZED COURANT-SNYDER THEORY
We now apply the technique developed in Sec. III to the Hamiltonian system in Eq. (11). Our goal is to find a new coordinate system where the transformed Hamiltonian vanishes. This idea is identical to that in Hamilton-Jacobi theory. Applications of Hamilton-Jacobi theory include the construction of action-angle variables for periodic systems [46] and finding geodesic curves on an ellipsoid [47]. It is often required that the variables in the Hamilton-Jacobi equation can be separated in order for the technique to be effective for practical problems. This limits its application. Since our dynamics is linear, the new coordinate system can be more easily constructed using the method developed in Sec. III. We will accomplish this goal in two steps. First, we seek a coordinate transformationz ¼ Sz such that, in thez coordinates, the Hamiltonian assumes the form where μðtÞ is a 2 × 2 matrix to be determined. To write Eq. (22) in the format of 2 × 2 blocks, we let and split the differential equation for S, i.e., Eq. (22), into four matrix equations, Including μðtÞ, we have five 2 × 2 matrices unknown. The extra freedom is introduced by the to-be-determined μðtÞ.
Based on the analogy with Eq. (10), we choose S 2 ≡ 0 to remove the freedom. We rename S 4 to be w, i.e., w ≡ S 4 , because it will be clear later that S 4 is the envelope matrix. Equations (26)-(29) become for matrices S 1 , S 3 , w, and μ. Because ðS 1 ; S 2 ¼ 0; S 3 ; S 4 ¼ wÞ describes a curve in Spð4Þ, they are consistent with the symplectic condition From Eq. (31), we obtain It is straightforward to verify that Eq. (30) is equivalent to another symplectic condition S 3 S T 4 ¼ S 4 S T 3 . Substituting Eqs. (32)-(35) into Eq. (32), we immediately obtain the following matrix differential equation for the envelope matrix w, This is the desired generalized envelope equation. It generalizes the 1D envelope equation (7), or the Ermakov-Milne-Pinney equation [28][29][30], as well as the previous matrix envelope equation for cases with only quadrupole and skew-quadrupole magnets, i.e., R ¼ 0 [21][22][23][24]. For n degrees of freedom, the envelope matrix w will be n × n, and the generalized envelope equation has the same form as Eq. (36).
Once w is solved for from the envelope equation, we can determine S 1 from Eq. (34) and S 3 from Eq. (33). In terms of the envelope matrix w, the symplectic transformation S and its inverse are given by The second step is to use another coordinate transformationz ¼ PðtÞz to transformH into a vanishing HamiltonianH ≡ 0 at all time, thereby rendering the dynamics trivial in the new coordinates. The governing equation for the transformation PðtÞ is As explained in Sec. III, the P matrix satisfying Eq. (39) is symplectic because JĀ ∈ spð4Þ. From μ ¼ μ T , we know that JĀ is also antisymmetric, i.e., JĀ belongs to the Lie algebra soð4Þ of the 4D rotation group SOð4Þ. Thus JĀ ∈ spð4Þ⋂soð4Þ, and PðtÞ is a curve in the group of 4D symplectic rotations, i.e., PðtÞ ∈ Spð4Þ⋂SOð4Þ ¼ Uð2Þ, provided the initial condition for PðtÞ is chosen such that Pð0Þ ∈ Spð4Þ⋂SOð4Þ ¼ Uð2Þ. We call PðtÞ the generalized phase advance, an appropriate descriptor in light of the fact that PðtÞ is a symplectic rotation. The Lie algebra element (infinitesimal generator) is the phase advance rate, and it is determined by the envelope matrix through Eq. (35). Since Spð4Þ⋂ SOð4Þ ¼ Uð2Þ, P and its inverse must have the forms Combining the two symplectic coordinate transformations, we obtain the transformation In thez coordinate representation, becauseH ≡ 0, the dynamics is trivial, i.e.,z ¼ const. This enables us to construct the symplectic matrix specifying the map between z 0 and z ¼ MðtÞz 0 as where subscript "0" denotes initial conditions at t ¼ 0, and P 0 is taken to be I without loss of generality. This expression for MðtÞ generalizes the decomposition of the symplectic map for the original 1D CS theory given by Eq. (10). The first and the third matrices in Eq. (10) obviously have the same construction as their counterparts in Eq. (43). The phase advance, as a 4D symplectic rotation P T in Eq. (43), generalizes the 2D rotation matrix, which is also symplectic, in Eq. (10). The phase advance P is generated by its infinitesimal generator JĀ, which is determined by the envelope matrix through μ ¼ ðwmw T Þ −1 . This mechanism for phase advance in 4D phase space is identical to the original 1D CS theory where the infinitesimal generator of the phase advance is w −2 for a scalar envelope w. The importance of the decomposition in Eqs. (43) and (10) can be appreciated from both physical and mathematical points of view. We explain the physical meaning of the decomposition here, and leave the mathematical analysis to Sec. V. The first matrix from the right is a matching transformation of the initial conditions to an equivalent focusing system where the phase space dynamics can be characterized by a timedependent rotation. The second matrix from the right is a transformation along the time axis in this equivalent focusing system, with the phase advance playing the role of a timelike evolution parameter. The third matrix from the right is a back-transformation to the original coordinate system at t > 0.
The coordinate transformation can also be used to construct invariants of the dynamics. A general description of linear symplectic invariants can be found in Refs. [48,49]. For any constant 4 × 4 positive-definite matrix ξ, the quantity is a constant of the motion, sincez ¼ PSz is a constant of motion. The subscript "ξ" in I ξ is used to indicate that it is an invariant associated with ξ. For the special case of ξ ¼ I, the phase advance P in Eq. (44) drops out, and where α, β, and γ are 2 × 2 matrices defined by β ≡ w T w; Here, we have used I CS to denote this special invariant because it is the invariant that generalizes the CS invariant [1] (or Lewis invariant [31,32]) for one degree of freedom in Eq. (8). The matrices, α, β, and γ are the generalized Twiss parameters in higher dimensions. It is straightforward to verify that they satisfy which is a familiar relationship in the original CS theory between the scalar Twiss parameters defined by Eqs. (4), (5), and (9). The symplectic condition wS T 3 ¼ S 3 w T has been used in obtaining Eq. (49).
It has been demonstrated that the envelope matrix w and the invariant I ξ define the beam dimensions and emittance for both low intensity beams and high intensity beams with strong space-charge potential [4,5].
Note that we have "overloaded" the symbols "M, w, α, β, γ, I CS " to represent the same physical quantities in both the original CS theory for one degree of freedom and the generalized CS theory in higher dimensions without causing any confusion. It is actually more appropriate to do so than not, because the quantities in higher dimensions recover their counterparts for one degree of freedom as special cases, and the correspondence between them is exact.

V. GROUP STRUCTURE OF THE GENERALIZED COURANT-SNYDER THEORY-ROTATION GAUGE AND MODIFIED IWASAWA DECOMPOSITION
In Sec. IV, we noted that initial conditions for the envelope matrix w need to satisfy the symplectic condition; otherwise they can be arbitrary. There are freedoms in the initial conditions and thus the solutions for w. But the transfer matrix M is independent from these freedoms, which are thus gauge freedoms. A subset of the gauge freedoms has the structure of the orthogonal group OðnÞ. For a time-independent element c ∈ OðnÞ, we define the gauge transformation c∶ ðw; PÞ↦ðw;PÞ as w ¼ cw; Let us show that the transformedw andP also satisfy Eqs. (36) and (39) From Eqs. (37) and (38), the S matrix and its inverse transform asS Combining Eqs. (50) and (51), (54) and (55), we conclude that M is invariant under the gauge transformation c∶ ðw; PÞ↦ðw;PÞ, i.e.,M ¼ M.
The OðnÞ gauge group introduces an equivalent class ½ðw; PÞ for the decomposition of M using ðw; PÞ. The dimension of this equivalent class is the dimension of M as a Spð2nÞ group. To specify S by w and _ w, 2n 2 numbers are needed. To specify P ∈ UðnÞ, additional n 2 numbers are needed. The symplectic condition for S, S T 3 S 1 ¼ S T 1 S 3 brings ðn 2 − nÞ=2 constrains on w and _ w, and the OðnÞ gauge freedom for the equivalent class is also ðn 2 − nÞ=2. The dimension of the decomposition is therefore the same as the dimension of Spð2nÞ.
According to the polar decomposition theorem, any nondegenerate square matrix X can be uniquely factored into an orthogonal matrix O and a symmetric, positivedefinite matrix Q, i.e., X ¼ OQ. As a matter of fact, Q ¼ ffiffiffiffiffiffiffiffiffi X T X which is in the form of a modified Iwasawa decomposition (or pre-Iwasawa decomposition), whose importance in phase space optics [37,38] and phase space quantum mechanics [39] has been recently realized. The modified Iwasawa decomposition is the unique decomposition of a 2n × 2n symplectic matrix G into the form where P ∈ Spð2nÞ⋂SOð2nÞ ¼ UðnÞ and Y is symmetric. The matrix Q is also symmetric, which is equivalent to the condition to be symplectic. These facts are also true if the decomposition is alternatively defined to be For a symplectic matrix, the modified Iwasawa decomposition is equivalent to the well-known Iwasawa decomposition for a semisimple Lie group [40]. Makingw symmetric at t ¼ t 1 fixes the OðnÞ gauge becausewðt 1 Þ is unique according to the polar decomposition theorem. However, such a choice only makesw symmetric at t ¼ t 1 . As in general gauge theories, we would like to pick a gauge such that the envelope matrix is symmetric for all t. To accommodate this desired symmetry, we need to modify the governing equations, especially the envelope equation. Let be the time-dependent polar decomposition of wðtÞ. Here we use uðtÞ to denote this specialwðtÞ, which is symmetric and positive-definite for all t. The matrix uðtÞ is the "symmetrized" wðtÞ and equals the square-root of the generalized β function. We will recast the envelope equation (36) in terms of β, as in the original Courant-Snyder theory for one degree of freedom [1]. The difference is that the procedure here has to be carried out in matrix form. We rewrite Eq. (36) as We symmetrize Eq. (60) by taking w T ½Eq:ð60Þ þ ½Eq:ð60Þ T w to obtain a second order ordinary differential equation for β, It is a second-order equation for β because _ w T _ w and w T _ w can be expressed in terms of β and _ β as follows. First note that Both _ u and D in Eqs. (64) and (65) can be expressed as functions of β and _ β. For _ u, from the definition of u we obtain whose left-hand side can be viewed as a linear operator on _ u associated with u, The properties of the linear operator L is discussed in the Appendix. Since u ¼ ffiffiffi β p is symmetric and positivedefinite, L u can be invertible to give where L −1 is the inverse of L defined in Eq. (A3). To express D in terms of β and _ β, we examine the symplectic condition Substituting in the polar decomposition w ¼ c −1 u gives Therefore, Equation (63) is a second equation for β. Its solutions do not uniquely determine the envelope matrix w, which is not surprising considering that β ¼ w T w is a "symmetric" version of w. However, due to the OðnÞ gauge freedom, β contains enough information to determine the transfer map M. In terms of u and c −1 , Even though the rotation matrix cðtÞ here is a function of t, the transformed phase advance is defined the same way as in the case of a global gauge, i.e., What is modified is the governing equation for P u , The second term on the right-hand side of Eq. (77) is due to the dependence on t of the local gauge. At last, the canonical coordinate transformation between z andz is and the transfer map is The symmetric decomposition furnished by Eqs. (75), (76), and (80) is equivalent to the decomposition described in Sec. IV, but it has three desirable features by comparison. The canonical coordinate transformation P u S u in Eq. (79) has the modified Iwasawa format for all t. It comprises a curve of the modified Iwasawa decomposition, developed from a dynamical point of view. The gauge freedom is removed, and the dimension of the symplectic transfer map is directly reflected by the dimension of the decomposition. At every t, MðtÞ is specified by two n × n symmetric matrices β and _ β, and a UðnÞ matrix P u . The dimension of MðtÞ is thus ðn 2 þ nÞ=2 þ ðn 2 þ nÞ=2 þ n 2 ¼ nð2n þ 1Þ.
Before ending this section, we emphasize that the purpose of studying the gauge freedom is to simplify the calculation of the symplectic map M and other lattice functions and beam parameters. By investigating the gauge freedom in the matrix envelope equation for w, we have found that we can actually bypass this gauge freedom and solve for the β matrix instead, which is symmetric and does not have the gauge freedom. From Eqs. (46), (48), and (49), the generalized Twiss parameters α and γ can also be expressed in terms of β and _ β. One important advantage of using the β matrix is that the symmetric matrices β and _ β form a linear space, which makes the numerical algorithms of searching for matched solutions for β much more efficient than for matched solutions for w.

VI. STABILITY ANALYSIS-SPECTRAL STABILITY AND STRUCTURAL STABILITY
The classical analysis by Courant and Snyder [1] on the instability induced by sum resonance for uncoupled transverse dynamics may give a wrong impression that coupling effects are always deleterious. The coupled dynamics can be stable or unstable depending on the specific configuration of the lattice, but certainly not more unstable than the uncoupled dynamics. The parameter space for a stable coupled lattice is probably much larger than that of a stable uncoupled lattice. In the conceptual design of the Möbius accelerator [50] and N-rolling lattice [24,51], it was argued that strongly coupled lattices are more preferable for highintensity beams. Strongly coupled systems have been implemented in the spiral line induction accelerator (SLIA) [52][53][54][55][56][57][58], which reached up to 10KA electron current at 5 MeV beam energy. Our understanding of the stability properties of coupled dynamics has been limited by the theoretical tools available. In this section, we demonstrate how the generalized Courant-Snyder theory can be applied to study the stability of the most general focusing lattice given by Eq. (11) with weak and strong coupling components in realistic accelerators.
For a thorough understanding, it is necessary to distinguish two types of linear stability (or instability). The first type is spectral stability, which means the linear dynamics is stable for all initial perturbations. The system is spectrally unstable if there exists an initial condition that grows without bond. In most contexts, the meaning of stability is that of spectral stability. The second type is the so-called structural stability (or strong stability). It mostly applies to systems that are spectrally stable. A spectrally stable system is structurally unstable if there is a spectrally unstable system infinitesimally close by. Otherwise, the spectrally stable system is also structurally stable. The wellknown result with respect to the stability properties of sum/ difference resonances for uncoupled lattices refers to the structural stability under the influence of an infinitesimal coupling component [1].
The spectral and structural stability of the transverse dynamics in a periodic focusing lattice is determined by its one-turn (or one-period) map MðTÞ. The fact that MðTÞ is a symplectic matrix regulates the stability properties in a significant way [41][42][43][44][45]. We list here the relevant results without presenting details of the proof.
The spectral property is determined by the eigenvalues and their multiplicities. There are four possibilities: (i) All eigenvalues are distinct and on the unit circle of the complex plane. (ii) All eigenvalues are on the unit circle. There are repeated eigenvalues. But the geometric multiplicity for all eigenvalues is the same as the algebraic multiplicity. (iii) All eigenvalues are on the unit circle. There are repeated eigenvalues with algebraic multiplicity greater than the geometric multiplicity. (iv) There exists at least one eigenvalue not on the unit circle. Cases (iii) and (iv) are spectrally unstable, and Cases (i) and (ii) are spectrally stable. For Cases (i) and (ii), we would like to know whether they are also structurally stable. It has been shown that Case (i) is structurally stable using the symplectic nature of MðTÞ [41][42][43][44][45]. Case (ii) needs to be subdivided into two categories: (iia) For all repeated eigenvalues, the corresponding eigenvectors have the same signatures. (iib) There is at least one repeated eigenvalue whose eigenvectors have different signatures. According to the Krein-Moser theorem [41][42][43][44][45], Case (iia) is structurally stable and (iib) is structurally unstable. For an eigenvector ψ of MðTÞ, its signature is defined to be the sign of its self-product hψ; ψi ¼ ψ Ã iJψ. The product between two eigenvectors ψ and ϕ in general is defined to be hψ; ϕi ≡ ψ Ã iJϕ, where ψ Ã denotes the complex conjugate of ψ T .
To design a coupled lattice, it is desirable to be in Case (i), which is both spectrally and structurally stable. As mentioned previously, for the general Hamiltonian given by Eq. (11), the parameter space satisfying this condition is large enough for most applications. Given a periodic lattice, we can search for a matched solution for β, as in the original Courant-Snyder theory for one degree of freedom [1]. After a matched β is found, the one-turn map is which implies that MðTÞ is similar to PðTÞ −1 . Their eigenvalues and multiplicity are identical. Because PðTÞ is a symplectic rotation, all of its eigenvalues are on the unit circle, automatically ruling out the unstable situation in (iv). The phase advance PðTÞ also determines the structural stability of the system. To prove this assertion, let ψ and ϕ be the eigenvectors of MðTÞ. Then S 0 ψ and S 0 ϕ are the eigenvectors of PðTÞ −1 , and where use had been made of the fact S 0 is symmetric, i.e., S T 0 JS 0 ¼ J. Equation (81) states that the signatures of eigenvectors of PðTÞ −1 and thus its structural stability are identical to that of MðTÞ.
These analyses lead to the important conclusion that the phase advance matrix PðTÞ completely determines both the spectral and structural stability of the general focusing lattices. This fact can significantly simplify the stability analysis in lattice design. For example, if the system is in Case (ii), we only need to look at the signatures of the eigenvectors of PðTÞ −1 to know if it is structurally stable. According the Krein-Moser theorem, if the eigenvectors for all repeated eigenvalues of PðTÞ −1 have the same signatures, then the system is structurally stable. Otherwise, it is structurally unstable. Let us show that this conclusion recovers the classical results on the stability properties of sum/difference resonances for uncoupled quadrupole lattices as special cases. In this case, the phase advance matrix is calculated to be [21,22] where ϕ x and ϕ y are the one-turn phase advance in the x− and y−directions. Its four sets of eigenvalues, eigenvectors, and signatures are (82) λ yþ ¼ cosϕ y þ isinϕ y ; ψ yþ ¼ ð0;1;0;iÞ T ; σ yþ ¼ −1; (84) λ y− ¼ cosϕ y − isinϕ y ; ψ y− ¼ ð0; 1; 0; −iÞ T ; σ y− ¼ 1: Resonance occurs when two or more eigenvalues collide, which has four possibilities: (a) Self-resonance in the x-direction. ϕ x ¼ nπ and λ x− ¼ λ y− . Case (a) is structurally unstable because σ xþ and σ x− are different. Case (b) is structurally unstable for the same reason.
For the sum resonance at the repeated eigenvalue λ xþ ¼ λ y− , the signatures σ xþ and σ y− of the corresponding eigenvectors have different signs. The sum resonance is thus structurally unstable. For the difference resonance at the first repeated eigenvalue λ xþ ¼ λ yþ , the corresponding eigenvectors ψ xþ and ψ yþ have the same signature. This is also true at the second repeated eigenvalue λ x− ¼ λ y− . The difference resonance is thus structurally stable. These results are well-known previously [1], but recovered here as a special case of a more general criterion based on the generalized phase advance and the Krein-Moser theorem [41][42][43][44][45]. We expect that the more general stability criterion expressed in terms of the phase advance matrix PðTÞ to be a powerful tool for future lattice design with strong coupling.

VII. CONCLUSIONS AND FUTURE WORK
We have presented in this paper a detailed derivation of the generalized Courant-Snyder theory for the most general linear focusing lattices with quadrupole, skewquadrupole, dipole, and solenoidal components, as well as torsion of the fiducial orbit and variation of beam energy. The theoretical structure of the theory in terms of gauge freedoms and group decomposition were described. We have also demonstrated the application of the theory in stability analysis for strongly and weakly coupled lattices. In addition to being more realistic, the most general Hamiltonian in Eq. (11) enables a much larger parameter space for designing strongly coupled lattices that are spectrally and structurally stable. The generalized Courant-Snyder parametrization scheme developed here provides an effective tool to study the coupled dynamics and to discover more optimized lattice design in the larger parameter space of general focusing lattices. The formalism also sets the theoretical foundation for investigating collective phenomena in high-intensity beams, such as the self-consistent solutions of the Vlasov-Maxwell equations in phase space including strong self-field effects that can couple the transverse dynamics [59][60][61][62][63].
As mentioned in Sec. IV, the theoretical framework developed here is valid for a linear system with any number of degrees of freedom. In particular, we can apply it to the 3D coupled dynamics, which includes the sychrotron oscillation in RF cavities, and the linear coupling between transverse and longitudinal dynamics as in the recent investigations of emittance change [64][65][66]. In this case, n ¼ 3 and the focusing matrix κðtÞ in Eq. (11) becomes a 3 × 3 matrix which describes both sychrotron and betatron oscillations as well as possible coupling between them. The envelope matrix w is 3 × 3 and satisfies Eq. (36) with R, m and κ being 3 × 3 matrices. The Twiss parameters α, β, and γ are 3 × 3 matrices, the symplectic matrix M is 6 × 6, and all the equations they satisfy are the same as in the case of 2-degrees of freedom. Studies in these directions will be reported in future publications.

APPENDIX: L x ðyÞ
In this Appendix, we derive the mathematical properties of the linear transformation used in Eqs. (68) and (72). Let A and X denote n × n matrices. For a symmetric, positivedefinite matrix A, define the linear function associated with A on X as L A ðXÞ ≡ AX þ XA: We prove that L A is invertible. It is enough to show that L A is injective, i.e., L A ðXÞ ¼ 0 only if X ¼ 0. Let X be in the kernel of L A , i.e., Since A is symmetric, the eigenvectors of A form a basis for vectors in R n . Expressed in this basis, Eq. (A1) is where u and v are any pairs of vectors in the basis, and λ u and λ v are the corresponding eigenvalues, respectively. Because the eigenvalues of a positive-definite matrix are positive, Eq. (A2) is possible only when v T Xu ¼ 0. This proves that L A is invertible. In terms of its components in this basis, L −1 A is given as where Y ¼ L A ðXÞ.