Stability properties of the transverse envelope equations describing intense ion beam transport

The transverse evolution of the envelope of an intense, unbunched ion beam in a linear transport channel can be modeled for the approximation of linear self-fields by the Kapchinskij-Vladimirskij (KV) envelope equations. Here we employ the KV envelope equations to analyze the linear stability properties of so-called mismatch perturbations about the matched (i.e., periodic) beam envelope in continuous focusing, periodic solenoidal, and periodic quadrupole transport lattices for a coasting beam. The formulation is analyzed and explicit self-consistent KV distributions are derived for an elliptical beam envelope in a periodic solenoidal transport channel. This derivation extends previous work to identify emittance measures and Larmor-frame transformations to allow application of standard form envelope equations to solenoidal focusing channels. Perturbed envelope equations are derived that include driving sources of mismatch excitation resulting from focusing errors, particle loss, and beam emittance growth. These equations are solved analytically for continuous focusing and demonstrate a factor of 2 increase in maximum mismatch excursions resulting from sudden driving perturbations relative to adiabatic driving perturbations. Numerical and analytical studies are carried out to explore properties of normal mode envelope oscillations without driving excitations in periodic solenoidal and quadrupole focusing lattices. Previous work on this topic by Struckmeier and Reiser [Part. Accel. 14, 227 (1984)] is extended and clarified. Regions of parametric instability are mapped, new classes of envelope instabilities are found, parametric sensitivities are explored, general limits and mode invariants are derived, and analytically accessible limits are checked. Important, and previously unexplored, launching conditions are described for pure envelope modes in periodic quadrupole focusing channels.


INTRODUCTION
Intense ion beams are usually transported in a periodic lattice of linear focusing elements that radially confine the beam against defocusing space-charge and thermal forces arising from the distribution of beam particles. Envelope equations are frequently used to model the low-order, statistical evolution of the beam edge. The KV envelope equations are coupled ordinary differential equations that describe the transverse evolution of the edge of an unbunched elliptical beam in a linear focusing lattice consistent with linear space-charge defocusing forces and constant emittances (i.e, invariant phase space projections) [1][2][3]. Even though the KV equations are self-consistent only for the singular and unphysical KV distribution function which has higher-order collective instabilities internal to the beam [4], the KV equations are valid for general distributions with evolving emittances in a rms equivalent beam sense [2,5,6]. Moreover, because space-charge profiles on intense beams tend to be fairly uniform in the core of the beam, the linear self-fields in the KV model should be accurate at high space-charge intensities except in a thin layer near the beam edge where the charge-density drops rapidly to zero. Consequently, for intense beams propagating in linear transport channels without significant particle loss or emittance growth, the KV envelope equations can provide a good estimate of the actual (rms statistical measure) envelope of the beam. Because of this and the fact that the ordinary differential equations describing the KV envelope evolution are straightforward to solve numerically, the KV envelope equations are extensively employed in the design of practical linear transport lattices.
The matched beam envelope is the solution to the KV envelope equations with the periodicity of the focusing lattice. The matched solution is expected to have smallest maximum radial excursions relative to other possible envelope evolutions in the lattice and it requires particular initial conditions in the envelope of beam particles. There will always be some finite mismatch error or deviation of the beam envelope from this matched evolution.
These mismatch errors can be analyzed in terms of the small-amplitude modes about the matched solution that are supported by the KV envelope equations [2,7]. Parametric instabilities of the mismatch modes must be avoided for practical machine operating points to maintain beam control. Even in cases where the envelope is stable, the structure of the mismatch modes is important. The frequencies of the envelope modes can resonate with lattice structures, particle orbits, and collective space-charge waves causing deleterious beam effects. For example, the high frequency breathing envelope mode is well-known to drive large-amplitude resonant beam halo causing loss of beam particles [8]. Therefore the control of this mode, even when stable, can be of critical importance.
In this paper we carry out a systematic parametric analysis of mismatch modes supported by the KV envelope equations for continuous focusing, periodic solenoid, and periodic quadrupole doublet linear transport lattices. Regions of strong parametric (band) instability for the periodic lattices are numerically mapped in a scaled manner rendering results immediately applicable to a large range of beam parameters and lattices. New classes of envelope instability are found and results are checked against analytical calculations of the envelope modes in the zero space-charge limit and in the thin-lens limit with a beam of maximal space-charge intensity. Important and often overlooked mode launching conditions and sources of envelope mismatch (driving terms) due to focusing errors, particle loss, and emittance evolution are analyzed. These launching conditions are of practical importance in understanding how to launch pure mode oscillations in experiments and simulations which can aid in understanding of complex effects.
This paper is organized as follows. In Sec. II, the envelope model employed is developed. The basic equations are reviewed, classes of transport lattices are defined, characteristic envelope responses are calculated, convenient parameterizations are identified, equations describing small amplitude driven perturbations about the matched beam envelope are derived, and general properties of envelope mode solutions are analyzed. Self-consistent KV distributions for solenoidal focusing channels are developed in Appendix B. These results are important to establish the validity of models employed. In Sec. III, the linear mode structure is completely solved in the continuous focusing limit employing a Green's function approach and the structure of the formal solutions are illustrated with specific examples to understand the differences between adiabatic, sudden, and harmonic driven perturbations. Next, normal envelope modes of periodic solenoidal and quadrupole doublet lattices are analyzed in Secs. IV and V. In both sections numerical results are presented that illustrate properties of: the matched envelope solutions, mode properties and symmetries, regions of parametric (band) instability, and pure mode launching conditions.

II. MODEL EQUATIONS
We consider an unbunched beam of ions of charge q and mass m coasting with axial relativistic factors β b = const and γ b = 1/ 1 − β 2 b . The ions are confined transversely by linear applied focusing fields that may result from a variety of focusing systems. Linear space-charge forces internal to the beam envelope are assumed, corresponding to a spatially uniform distribution of charge within an elliptical beam envelope with beam image charges induced on any aperture structures neglected.

A. KV envelope equations
As illustrated in Fig. 1, the transverse beam cross section at axial coordinate s is taken to have an s-varying elliptical cross-section with envelope (i.e., edge) radii r x (s) and r y (s) along the principal x-and y-coordinate axes. The charge distribution is centered at x = y = 0 with uniform uniform space-charge within the ellipse (x/r x ) 2 + (y/r y ) 2 = 1, and zero outside. The x-and y-equations of motion for the orbit of a single-particle within this uniform density elliptical beam are given by [1,2] x (s) + κ x (s)x(s) − 2Qx(s) [r x (s) + r y (s)]r x (s) = 0, y (s) + κ y (s)y(s) − 2Qy(s) [r x (s) + r y (s)]r y (s) = 0, where κ x (s) and κ y (s) represent linear applied forces of focusing optics in the transport lattice, and is the dimensionless beam perveance representing self-field defocusing forces internal to the beam. Here, λ = const is the beam line-charge density, c is the speed of light in vacuuo, and 0 is the permitivity of free-space. In the limit of negligible beam space-charge, Q → 0 in Eq. (1), and the particle moves solely in response to the applied focusing forces of the lattice. For a uniform density elliptical beam, the envelope radii r x and r y are connected to statistical moments of the distribution by where · · · ⊥ denotes a transverse statistical average over the beam phase-space. Differentiating Eq. (3) with respect to s and employing the equations of motion (1) yields coupled ordinary differential equations describing the evolution of the envelope radii r x (s) and r y (s) [1,2,7] r j (s) + κ j (s)r j (s) − 2Q r x (s) + r y (s) − ε 2 j r 3 j (s) = 0.
Here and henceforth j ranges over both transverse coordinates x and y, and are the rms emittances which provide statistical measures of the beam phase-space area projections in the x-x and y-y phase-spaces. For charge distributions that are constant on nested elliptical surfaces, the envelope equations (4) are valid with evolving emittances and can be solved for the evolution of r j (s) provided consistent emittance evolutions are employed [6]. In this situation the envelope equations are often referred to as the rms envelope equations and can be applied to a wide variety of nonuniform laboratory beams through an "equivalence" to a uniform density beam with equal second-order statistical moments [2,5,6]. For the case of a KV distribution with uniform space-charge, the ε j are constants and the envelope equations are referred to as the KV envelope equations [1][2][3]. The corresponding self-consistent KV distribution under a Vlasov model is singular and unphysical and has numerous higher-order collective instabilities internal to the beam [4,9]. However, in cases of reasonably smooth initial beam distributions and a transport channel tuned for applied focusing linearity and a lack of instabilities, the emittance evolutions will be small, the linear space-charge model will be accurate, and the KV envelope model with constant emittances can be reliably applied. For a particular focusing lattice [i.e., specified focusing functions κ j (s)] and beam parameters (Q and ε j ), the KV envelope equations (4) are typically integrated from an initial condition (i.e., specified r j and r j at s = s i ) to solve for the envelope evolution. In many cases we will take ε x = ε y ≡ ε = const corresponding to a beam with isotropic transverse temperature (on average) and negligible emittance evolution propagating in a linear transport channel without bends, dispersion, and edge-focusing.

B. Periodic transport channels and beam matching
In a periodic transport lattice, applied focusing elements are arranged in a periodic sequence with lattice period L p and lattice focusing functions κ j (s) satisfying κ j (s + L p ) = κ j (s). (6) The matched beam envelope is the solution to the KV envelope equations (4) in a periodic transport where the envelope radii r j (s) are also periodic functions with period of the lattice, i.e., r j (s) = r jm (s) with r jm (s + L p ) = r jm (s).
For specified focusing functions κ j (s), beam perveance Q, and emittances ε j , this is equivalent to requiring that r j and r j satisfy specific initial conditions at s = s i when the envelope equations (4) are integrated as an initial value problem. Required initial conditions will change with the phase of s i within the lattice period. The matched envelope solution is important because it is believed to have minimum radial excursions in the transport channel relative to all other possible envelope evolutions [10], thereby allowing use of more radially compact transport channels. Although formulations in the following parts of this section apply for arbitrary periodic focusing functions κ j (s), in Secs. III-V we analyze special classes of (a) continuous, (b) periodic solenoidal, and (c) periodic quadrupole doublet focusing lattices described below with κ j (s) that are piecewise constant, corresponding to so-called hard-edge or square-edge models. These three choices are representative of broad classes of lattices used in practical applications and can guide analysis of more complicated lattices. Examples of matched envelope solutions will be presented for each class of lattice in Secs. III-V. These hard-edge models can be applied to a wide range of periodic solenoid and quadrupole lattices with fringe fields where the κ j (s) vary smoothly in s by using equivalent hard-edge replacement prescriptions (see Appendix A).
(a) Continuous focusing, with equal and constant focusing in each plane [see Fig. 2(a)] [2]. In this case, κ x (s) = κ y (s) = k 2 β0 = const > 0. (8) Continuous focusing is equivalent to a partially neutralizing, non-interacting background of charges and is commonly used to simply estimate the average focusing properties of periodic lattices in rapid design estimates.
Here, k β0 is the wavenumber of particle oscillations in the absence of space-charge (see Sec. II D). The choice of "lattice period" L p for continuous focusing is arbitrary.
Here, the envelope equations (4) must be interpreted in a local rotating "Larmor" frame in the sense shown in Appendix B. In analyses of solenoidal focused systems, we will assume that all calculations carried out in the Larmor frame unless otherwise noted. We define a hard-edge lattice where the solenoids have axial length ηL p with κ =κ = const > 0 within the solenoid and zero outside and the solenoids are separated by axial drifts of length d = (1 − η)L p . Here, η ∈ (0, 1] is the occupancy factor of the solenoid in the lattice period. In the Larmor frame, solenoidal focusing is equivalent to continuous focusing in the limit η → 1. (c) Quadrupole doublet focusing in an alternating gradient lattice with period L p [see Fig. 2(c)] [2]. In this case, The quadrupole focusing fields are commonly derived either from electrostatic or magnetostatic lenses. In terms of a simple 2D field model, the focusing strength for electric lenses can be expressed as where [Bρ] ≡ γ b β b mc/q is the particle rigidity and dE x /dx is the (generally s-varying) linear quadrupole field gradient of the electric field. Similarly, the focusing strength for magnet lenses can be expressed as where dB x /dy is the (generally s-varying) linear quadrupole field gradient of the magnetic field. We define a hard-edge doublet lattice where the focusing (κ q > 0 for focusing in x) and defocusing (κ q < 0) quadrupoles have axial length ηL p /2 with κ q = ± κ q with κ q = const within the quadrupoles and zero outside, and are separated by (generally unequal) axial drifts of length d 1 = α(1 − η)L p and d 2 = (1 − α)(1 − η)L p . Here, η ∈ (0, 1] is the quadrupole occupancy factor in the lattice period, and α ∈ [0, 1] is a syncopation factor that measures how close the focusing and defocusing quadrupoles are to each other in the lattice period. For the special case of α = 1/2 corresponding to a symmetric FODO lattice, the drifts are equal [d 1 = d 2 = (1 − η)L p /2] and the lenses are equally spaced. For α = 1/2, the lattice is called syncopated. Without loss of generality, we restrict analysis to 0 ≤ α ≤ 1/2 because values of α with 1/2 ≤ α ≤ 1 can be mapped to the range 0 ≤ α ≤ 1/2 by relabeling lattice quantities.

C. Characteristic Envelope Responses
The envelope evolution in the periodic, hard-edge focusing lattices described in Sec. II B can be interpreted in terms of the sequential response of the envelope to the focusing action of optical elements followed by free-drift expansions without focusing (κ j = 0).
No general analytical solutions to the envelope evolution within hard-edge solenoidal or quadrupole optics are known for arbitrary values of Q and ε j . Some simplifications can be made using the envelope Hamiltonian as a constant of the motion within a particular optical element, but full integrations of the envelope trajectories r j (s) are not known. However, the essential effect of the optics can be simply understood in the thin-lens limit where η → 0 while adjusting the focusing strength such that the desired net particle focusing (measured by σ 0j , see Sec. II D) is produced. This is effected by taking where δ(s) is the Dirac delta-function, f = const is the thin-lens focal length, and s = s o is the axial location of the equivalent model optical impulse. In this thin-lens limit, it follows from Eqs. (4) and (13) that the action of the optic is to transform the envelope as where | s ± o denotes the limit lim s→so from above and below s = s o . For solenoidal and quadrupoles that are focusing in the j-plane, f > 0, whereas for defocusing quadrupoles in the j-plane, f < 0.
General analytical solutions for the envelope evolution in drift regions are also unknown, though some simplifications can be made using envelope Hamiltonian invariants in particular drift regions. However, the range of drift envelope evolution can be analytically understood in the extremes of an emittance-dominated beam with Q = 0 and a spacecharge dominated beam with ε j = 0. First, in the limit of a emittance-dominated beam, Q = 0, and the j-plane envelope equation (4) reduces to This equation is straightforward to integrate using the constancy of the envelope Hamiltonian, r 2 j /2+ε 2 j /(2r 2 j ) = const. After some manipulations we obtain Here, s = s i is the initial condition where r j = r j (s i ) and r j = r j (s i ). In the opposite limit of a space-charge dominated beam, ε j = 0, the envelope equations can be cast in decoupled form using sum and difference coordinates defined as r ± = (r x ± r y )/2 to give The constancy of the envelope Hamiltonian, r 2 + /2 − Q ln r + = const can be employed to simplify integration of the equation for r + and the equation for r − is trivially integrable. We obtain (see Appendix C) Here, erfi(z) = erf(iz)/i = (2/ √ π) z 0 dt exp(t 2 ) is the imaginary error function, i ≡ √ −1, and notation on initial conditions is analogous to that employed in Eq. (16). An alternative form of the solution for r + has been presented by Humphries [11].
To contrast these two limiting free-drift solutions, the envelope expansions r x (s)/r x (s i ) and r + (s)/r ± (s i ) given by Eqs. (16) and (18) are plotted in Fig. 3 for equal initial forces Q/r + (s i ) = ε x /r 3 x (s i ) and zero initial angles. The space-charge dominated expansion is more rapid due to the ∼ 1/r + expansion force being larger than the ∼ 1/r 3 x expansion force at large radii.

D. Single-particle phase advances
Because it is useful to employ single-particle phase advances to parameterize matched solutions to the KV envelope equations (4), we review these phase advances here in order to define notations employed and summarize results needed later.
where M x (s|s i ) denotes the 2 × 2 transfer matrix from the initial x-plane particle phase-space coordinates at axial coordinate s = s i to the phase-space coordinates at s. An analogous equation holds for the y-plane orbit. An eigenvalue analysis of M j (s|s i ) shows that orbit stability corresponds to and that Tr M j (s i + L p |s i ) is independent of the initial axial coordinate s i . A stable orbit of the particle satisfying Eq. (1) can be expressed in terms of the phase-amplitude formulation as [12,13] x where A j = const and the amplitude and phase functions w j and ψ j satisfy For a periodic lattice with κ j (s + L p ) = κ j (s), without loss of generality, w j can be taken to be the positive, real solution to the amplitude equation satisfying w j (s + L p ) = w j (s). Comparing the amplitude equation (22) with the envelope equation (4) for a matched beam, we we observe that r jm and w j are related as The phase advance σ j of the orbit through one lattice period can then unambiguously defined by integrating ψ j = 1/w 2 j to yield independent of s i . Analysis of the transfer matrix elements then shows that Direct application of this formula to calculate σ j requires proper branch identification. For the case of an isotropic beam with ε x = ε y ≡ ε and the three classes of lattices outlined in Sec. II B, the phase advances of the x-and y-orbits will be equal, i.e., σ x = σ y . In this symmetric situation, we will not distinguish between σ x and σ y , and we denote σ x = σ y ≡ σ.
In the limit of zero space-charge (Q → 0), we denote the phase advance σ j calculated from Eq. (24) as The so-called undepressed phase-advance σ 0j measures the phase advance of particles in response to the linear applied focusing fields and provides a measure of the focusing strength of the lattice. For periodic focusing lattices with κ x = κ y or κ x = −κ y , σ 0x = σ 0y , and we denote σ 0x = σ 0y ≡ σ 0 . Equations (1) and (25) can be employed to relate σ 0 = lim Q→0 σ to the lattice focusing functions κ j . For the three classes of lattices described in Sec. II B, we summarize expressions for σ 0 . First, for the continuous focusing lattice, κ x = κ y = k 2 β0 = const and there is no strict periodicity. However, for purposes of comparisons to periodic systems (see Fig. 2), we can define the undepressed particle phase advance σ 0 over an arbitrary "period" L p by For the hard-edge solenoidal and quadrupole doublet focusing channels illustrated in Fig. 2, it can be shown that for the solenoidal channel with Θ ≡ √κ L p /2, and for the quadrupole channel with Θ ≡ | κ q |L p /2. For any periodic channel with specified focusing functions κ j (s) (possibly including fringe variations), particular values of σ 0 can be produced by scaling the amplitude of the focusing functions, i.e., κ x → kκ x and κ y → kκ y with k = const constrained by σ 0 . For example, consider the hard-edge quadrupole lattice. At specified lattice period L p , occupancy η, and syncopation factor α, Eq. (29) constrains the value of the focusing strength | κ q | in terms of σ 0 . The thin-lens limit of Eqs. (28) and (29) can be explored by settingκ = η/(f L p ) for solenoids and κ q = η/(f L p ) for quadrupoles and then taking the limit η → 0. This gives These results can also be derived directly from a transfer matrix approach of the undepressed orbits using Eq. (25) in the limit Q → 0 to show that for solenoidal and quadrupole channels, respectively. When σ 0x = σ 0y and σ x = σ y , the ratio of space-charge depressed to undepressed phase advances σ/σ 0 provides a convenient, normalized measure of space-charge intensity which we call the space-charge depression. For a fixed value of σ 0 , the limit σ/σ 0 → 0 corresponds to a cold-beam at the space-charge limit with zero emittance (i.e., Q = 0 and ε = 0, or Q → ∞ for finite ε), and the limit σ/σ 0 → 1 corresponds to an emittance-dominated beam with zero space-charge effects (i.e., ε = 0 with Q = 0, or ε → ∞ for finite Q). Stability properties of the envelope equations (4) will be characterized in terms of σ 0 and σ/σ 0 for stable undepressed orbits with σ 0 ∈ (0, π) and the full range of space-charge intensity σ/σ 0 ∈ [0, 1].
Matched envelopes can also exist in periodic transport channels for ranges (bands) of σ 0j with σ 0j > π. These bands will correspond to higher focusing strength (e.g., largerκ for the hard-edge solenoidal lattice) and will be separated from the stable range σ 0j ∈ (0, π) by so-called stop bands where (1/2) |Tr M j (s i + L p |s i )| > 1 and the single-particle orbits are unstable. Successive stable bands will have narrower ranges of σ 0j and the first stop band begins at σ 0j = π. Higher stable bands tend to be more sparse and narrower for quadrupole focusing as opposed to solenoidal focusing, as can be seen from analyzing regions where |cos σ 0 | < 1 from Eqs. (28) and (29). This occurs because stronger quadrupole focusing in alternating gradient lattices tends to produce large excursions within the defocusing quadrupoles, leading to narrow parameter bands for stable orbit bundles. Higher stable bands in solenoidal channels are wider than for quadrupole channels because such systems are focusing in all planes. However, contrary to assertions by some authors [14][15][16], we find that machine operation in higher bands of stability for solenoidal focusing will not, in general, prove advantageous. This follows because the higher bands require significantly larger focusing strengths, while for fixed beam perveance Q, these stronger optics also result in matched beams with larger envelope excursions than for σ 0j < π, leading to larger apertures needed for beam confinement [17]. Moreover, matching the beam envelope to transport at high phase advances from lower phase advances can prove problematic. Because of these problems we only consider σ 0j ∈ (0, π) in this paper. Contrary to σ 0j > π, we find that for fixed Q in periodic focusing lattices with high degrees of symmetry, that the matched beam envelope generally is radially more compact for higher values of σ 0j in the range σ 0j ∈ (0, π) [18,19].

E. Scaled envelope equations and system parameters
Solutions for the envelope equations (4) can be parameterized in terms of the lattice focusing functions κ j , the perveance Q, and the emittances ε x = ε y = ε. For example, in the hard-edge solenoid and quadrupole doublet lattices defined in Sec. II B, the functions κ j are described in terms of the lattice period L p , the occupancy η, the focusing strengthκ (solenoids) or κ q (quadrupoles), and the syncopation parameter α (quadrupoles only). Thus for hard-edge solenoids, five "direct" parameters (L p , η,κ, Q, and ε) can be employed to specify the matched beam solution, while for the hard-edge quadrupole doublet six direct parameters (L p , η, α, κ q , Q, and ε) can be employed. Other choices of hard-edge solenoidal and quadrupole lattices yield similar direct parameters. Focusing lattices including axial fringe field models of the elements introduce additional direct parameters that replace the hard-edge occupancy η and focusing strengthsκ and κ q .
Rather than employing the direct parameters, matched beam solutions of the envelope equations can be described with a reduced number of dimensionless parameters [20,21]. The particle phase advance parameters σ 0 and σ (or σ 0 and σ/σ 0 ). This procedure can be carried out by introducing scaled coordinates for continuous focusing, or general solenoidal and quadrupole transport lattices as follows. For equal emittances and finite perveance (i.e., ε x = ε y ≡ ε and Q > 0), the envelope equations (4) can be expressed for the three classes of lattices as Here, S = s/L p is a dimensionless axial coordinate that increases by unity as the beam advances over a lattice period L p , are dimensionless scaled envelope radii, for continuous or solenoidal focusing ξ j = 1, for quadrupole focusing ξ x = +1 and ξ y = −1, F (S; σ 0 , {S}) is a scaled focusing function that depends on S, σ 0 , and a set of "shape" parameters {S} of the focusing function, and is a dimensionless parameter measuring relative space-charge intensity. For a matched beam envelope (i.e., r j = r jm ), the parameter Γ can be related to the depressed particle phase advance σ through Eq. (24) as The parameters of the scaled envelope equations (31) are σ 0 , the shape parameters {S} which describe the action of the applied focusing forces on a single-particle, and the dimensionless parameter Γ defined in Eq. (34), or alternatively to Γ, σ or σ/σ 0 .
Specific examples of the scaled focusing function F (S; σ 0 , {S}) for the hard-edge solenoidal and quadrupole lattices described in Sec. II B clarify the forms given in Eqs. (31)- (34). For the solenoidal lattice, we take F (S; σ 0 , {S}) = F (S; σ 0 , η) = [L 2 pκ ][κ x (s)/κ] where [κ x (s)/κ] is a function of S and η, and [L 2 pκ ] is a function of σ 0 and η through an application of the phase advance equation of constraint (28). Similarly, for the quadrupole lattice, we take is a function of S, η, and α, and [L 2 p | κ q |] is a function of σ 0 , η, and α through an application of the phase advance equation of constraint (29). For the special case of continuous focusing, there are no shape parameters (i.e., {S} = ∅), and the choice of lattice period L p entering in σ 0 and σ is arbitrary. Hence, matched solutions in continuous focusing can be regarded as a function solely of σ/σ 0 (see also Sec. III A).
From the results above, it follows that the matched beam envelope in the hard-edge solenoidal transport lattice considered can be described in terms of three parameters: the magnet occupancy η, the undepressed particle phase advance σ 0 , and the space-charge depression σ/σ 0 rather than five direct parameters [20]; whereas matched solutions to the quadrupole doublet lattice can be described in terms of four parameters: η, σ 0 , σ/σ 0 , plus the syncopation factor α (α ≡ 1/2 for the special case of a symmetric FODO lattice) rather than the five direct parameters [21]. In subsequent analyses we will exploit this reduction by specifying matched beam solutions in terms of these reduced parameter sets. However, to maintain direct connection the usual form of the envelope equations we will carry out most analyses in terms of conventional dimensioned envelope coordinates and direct parameters.

F. Small amplitude perturbations about the matched beam envelope
To analyze the evolution of perturbations about the matched beam envelope, we resolve the envelope coordinates r j evolving according to Eq. (4) from general initial conditions as r j (s) = r jm (s) + δr j (s).
Here, r jm is a matched solution to the envelope equation and δr j are perturbations about the matched beam envelope. In the small-amplitude limit we require that |δr j | r jm over the range of s analyzed. The envelope perturbations δr j are typically referred to as mismatch oscillations about the matched beam envelope.
In addition to the perturbations in the envelope coordinates δr j , we consider additional driving perturbations that can excite the envelope perturbations from the matched beam condition. Specifically, we allow the lattice focusing functions κ j , the beam perveance Q, and the rms emittances ε j to vary by making the following substitutions in the envelope equations (4) Here, δκ j (s) represent perturbations to the focusing functions κ j (s) used in generating the matched beam envelope, δQ(s) with −Q ≤ δQ ≤ 0 represents the effect of particles lost outside the beam envelope, and δε j (s) with ε j +δε j > 0 model evolutions in the rms emittances. The focusing perturbations δκ j (s) need not be periodic and can represent construction and excitation errors in the focusing optics. The long-path trend of the perturbed perveance δQ(s) ≤ 0 will be negative due to halo particles being scraped, but it can locally increase if any halo particles re-enter the core distribution. Similarly, over a long path the emittance perturbations δε j (s) will tend to be positive and increasing due to the relaxation of the beam in response to nonlinear applied focusing fields, space-charge nonlinearities, etc., but can locally decrease because rms emittances are statistical measures of phase-space area (e.g., momentary unwrapping of strong "S"-shaped distortions in the particle phase-space, etc.). In the small amplitude regime we require that δQ Q and that δκ j and δε j be small in the sense that they induce small corrections on the particle orbits relative to κ j and ε, respectively. Higher-order theories outside the scope of the envelope model are generally necessary to calculate δQ(s) and δε j (s) for use in this formulation.
Inserting Eqs. (35) and (36) into the envelope equation (4) and neglecting all terms of order δ 2 and higher yields the linearized perturbation equations δr j + κ j δr j + 2Q (r xm + r ym ) 2 (δr x + δr y ) + These equations can be equivalently expressed in first-order matrix form as where is a vector of perturbed envelope phase-space coordinates, is a 4 × 4 matrix with periodic s-dependence of period L p , and is a vector of driving perturbations. The solution to Eq. (38) can be expressed as where δR h is the general solution to the homogeneous equation (38) with δP = 0, and δR p is any particular solution for the full equation with general δP. Henceforth, the homogeneous and particular components of the solution will be analyzed in turn without explicitly distinguishing them to avoid cumbersome notation. For the special cases of continuous (κ x = κ y = k 2 β0 = const) and periodic solenoidal focusing (κ x = κ y = κ) channels with a round beam equilibrium (ε x = ε y = ε and r xm = r ym = r m ), considerable simplification results by resolving the envelope perturbations into sum and difference variables as In these cases, the perturbed envelope equations (37) can be expressed as and are decoupled. Here, κ = k 2 β0 for continuous focusing. Each of the decoupled equations for δr ± in Eq. (45) can be cast into matrix form analogous to Eqs. (38)-(43). However, in this case the vectors and matrices will be of dimension 2 rather than dimension 4. From the structure of Eq. (45), note that changes in beam perveance (i.e., δQ = 0) drive only δr + and not δr − , whereas changes in focusing strength (δκ j = 0) and emittance (δε j = 0) can project on δr + for in-phase error components in x and y, and δr − for out of phase errors. Also, due to the additional factor of (2Q/r 2 m )δr + , the restoring force coefficient ∝ δr ± is larger in the equation for δr + than in the equation for δr − , and oscillations in δr + will be more rapid than oscillations in δr − for Q = 0.
The homogeneous solution to the linear equations (38) describes the evolution of normal mode perturbations launched by mismatch in the initial conditions at s = s i (i.e., δR = 0 at s = s i ), whereas the particular solution describes how mismatch oscillations can be driven by changes in beam and focusing parameters. The normal modes supported by the homogeneous equations have been previously analyzed for both continuous and periodic focusing channels by Struckmeier and Reiser [7]. In this paper, this previous work is extended to carefully identify classes of modes supported over a wide range of system parameters.
The homogeneous solution to Eq. (38) can be expressed in terms of a transfer map as where M e (s|s i ) is a 4 × 4 matrix that maps the perturbed envelope coordinates δR from an initial axial coordinate s = s i to position s. Stability and launching conditions of normal modes can be analyzed in terms of eigenvalues and eigenvectors of the transfer map M e [12,13,22]. The four complex eigenvectors E n and eigenvalues λ n (not to be confused with the line-charge λ) of the envelope transfer matrix through one lattice period M e (s i + L p |s i ) from an initial condition s i are defined according to where n = 1, 2, 3, and 4. The overall scale of the eigenvectors E n is arbitrary. Without loss of generality, we take all eigenvectors to be normalized according to E n · E * n = 1, where * s denote complex conjugation. The eigenvalues λ n are independent of the location of the initial axial coordinate s i , whereas the eigenvectors E n (s i ) vary with the relative location of s i in the lattice period [22]. Stability in the sense of bounded excursions of perturbations in phase-space correspond to λ n on the complex unit circle with |λ n | = 1 for all n, whereas instability results from any λ n lying off the complex unit circle (i.e., at least one unstable mode).
Because the envelope perturbation equations are real and Hamiltonian, M e is a real symplectic matrix. Consequently, the four eigenvalues λ n of M e occur both as complex conjugate pairs (λ n , λ * n ) and reciprocal pairs (λ n , 1/λ n ) and can be categorized according to four symmetry classes (not counting possible degeneracies) described by Dragt [12,22]. These classes with associated eigenvalue and eigenvector symmetries are illustrated in Fig. 4. Eigenvalues are expressed in polar form as λ n = γ n exp(iσ n ) with real γ n and σ n and the specific labellings in n indicated. Standard analyses [12,13,22] show that γ n is the growth factor of the mode per lattice period with |γ n | = 1 corresponding to stable oscillations, and |γ n | > 1 and |γ n | < 1 corresponding to oscillations with exponential growth and damping, respectively. The number of e-folds that the amplitude of mode oscillations grows per lattice period is ln |γ n |. It is shown later that σ n is related to the phase advance of the mode oscillations per lattice period. For each of the four classes of eigenvalue symmetries in Fig. 4, two real numbers are needed to describe the four eigenvalues [e.g., σ 1 and σ 2 in case (a)]. In case (a), |λ n | = 1 and the eigenvectors corresponding to eigenvalues λ n and 1/λ n are simply related as complex conjugates with E 3 = E * 1 and E 4 = E * 2 . For |λ n | = 1, some relation between eigenvectors associated with the eigenvalues λ n and 1/λ n (i.e., through operations with M e , M −1 e , symplectic group generators, conjugation, etc.) would be useful rather than calculating them both directly from the characteristic equation defining the eigenvalues from Eq. (47). Unfortunately, it was shown that no such relationship exists in cases (c) and (d) [e.g., E 1 and E 3 are not simply related in case (c)]. Although such a relationship may in principle exist in case (b), efforts to derive it failed.
In Fig. 4, eigenvalues in case (a) correspond to stability (γ n = 1), whereas cases (b)-(d) are unstable. Case (b) is called a "confluent resonance" instability because it represents a parametric resonance between both envelope mode oscillation frequencies and the lattice as evident by the eigenvalue symmetry. Cases (c) and (d) are called "lattice resonance" and "double lattice resonance" instabilities because they represent a half-integer parametric resonance between the focusing structure and one or both mode oscillation frequencies. Transitions between cases lead to degenerate eigenvalues which can result in either linear growth in oscillation amplitude with periods traversed or no instability [22]. Which case the degeneracy falls into requires further analysis and is not addressed in this paper (other than a special case in the thin-lens limit) because the degenerate cases correspond to lower-dimensional surfaces in accessible parameter-space that are launched with probability zero. For the types of lattices in Sec. II B, we find in Secs. III-V that the eigenvalues λ n fall into cases (a)-(c) and case (d) does not occur.
Any initial envelope perturbation at s = s i can be expanded in terms of the eigenvectors as where the α n are complex constants subject to δR(s i ) being real. To launch pure envelope modes with specific phase advances and growth factors for each of the eigenvalue classes listed in Fig. 4, it is convenient to resolve the initial Eigenvalues Eigenvectors condition δR(s i ) = δR (s i ) into pure normal modes as summarized in Table I. Mode structures are illustrated in the table by mappings of the pure mode initial condition through one lattice. The eigenvectors E n , and consequently the pure mode launching conditions δR (s i ), vary with the relative location of the initial condition s i in the lattice period. Note that cases (a) and (b) in Fig. 4 have two distinct modes, while cases (c) and (d) have three and four distinct modes. The evolution of the beam mismatch δR(s) can also be interpreted as evolving linear projections of envelope  Fig. 4. Here, A and ψ are real-valued amplitudes (unnormalized) and phases of mode ( = 1, 2, and possibly 3 and 4 under conventions taken), and σ and γ > 1 are the phase advances and growth factors of the mode over one lattice period. Symbols are abbreviated with δR ≡ δR (si), E ≡ E (si), and Me ≡ Me(si + Lp|si).
mismatch onto pure normal modes according to:

cases (a) and (b),
where A and ψ are real-valued. This interpretation allows us to unambiguously define that pure mode phase advance σ as the change in ψ as s advances over one lattice period L p , thereby providing a concrete branch selection criteria to relate the complex eigenvalues λ n to the mode phase advances σ . Alternatively, the correct branch for the mode phase advance advance σ can be selected by numerically integrating Eq. (38) with δP = 0 from the pure mode initial conditions in Table I to calculate the orbit δR (s) , and then Fourier transforming δr x (s) or δr y (s) to identify frequency components. A sufficient propagation distance should be taken to allow clear branch identification from discretized transforms. A specific initial envelope mismatch δR(s i ) can be resolved into normal mode projections according to Eq. (49) by solving for the four (in total) real-valued amplitude A and phase ψ (s i ) parameters needed. Values of A cos ψ (s i ) and A sin ψ (s i ) needed [ψ = π for = 2, 3 in case (c) and for all in case (d)] can be found by matrix inversion of Eq. (49). A pure mode of a given amplitude can be launched by setting one A in Eq. (49) to the desired amplitude and setting all other A = 0 . The pure mode phase ψ can be varied over 2π [except modes that do not have phases: 2 and 3 in case (c) and all in case (d)] to sweep through all possible phase launches. Individual mode launching conditions are not, in general, orthogonal [e.g., . However, in the stable case (a) the two modes are almost orthogonal in a sense that lim s→∞ (1/s) For the special cases of continuous and periodic solenoidal focusing channels with a round beam equilibrium (ε x = ε y = ε and r xm = r ym = r m ), the mode structure is simpler to analyze in terms of the decoupled equations (45) for δr ± = (δr x ± δr y )/2 rather than the coupled equations (37) for δr x and δr y . The transfer matrix M e that advances δR = (δr + , δr + , δr − , δr − ) through one lattice period is of block diagonal form with where M ± (s i + L p |s i ) are 2 × 2 symplectic matrices that can be independently analyzed for the stability of δR ± = (δr ± , δr ± ). The analysis for the normal modes and launching conditions supported by the uncoupled homogeneous equations is analogous to the coupled case above. For each uncoupled equation for δr ± we need only analyze a 2 × 2 eigenvalue equation Each of these reduced equations defines a pair of eigenvalues: λ + and 1/λ + for the δr + equation, and λ − and 1/λ − for the δr − equation. Stable modes with λ ± on the unit circle correspond to [13] 1 2 |Tr M ± (s i + L p |s i )| < 1.
Possible reduced symmetry classes of the eigenvalues are illustrated in Fig. 5 [22]. It follows that for a continuous or periodic solenoidal focusing channel with a round beam equilibrium there can only be a single stable mode (phase advance σ ± ) or two modes: one with pure exponential growth, and one with pure damping (both of the lattice resonance class with phase advance σ ± = π). The situation in Fig. 4(b) corresponding to a confluent resonance is excluded when the 4 × 4 coupled envelope equations reduce to two 2 × 2 decoupled equations. Initial envelope perturbations in x and y can be projected onto the sum and difference perturbations at axial coordinate s = s i as For pure mode launches, one takes This launching condition is schematically plotted in Fig. 6. Note that the pure mode in δr + with δr − = 0 has in-phase perturbations with δr x = δr y (breathing distortion), whereas a pure mode in δr − with δr + = 0 has π out of phase perturbations with δr x = −δr y (quadrupole distortion). For this reason, pure decoupled modes in δr + and δr − are called "breathing" and "quadrupole" modes, respectively.

b) Unstable, Lattice Resonance
Re Im 1 Re In the limit of zero space-charge (Q → 0), the perturbed envelope equations (37) decouple and the envelope modes can be analyzed independently in the x-and y-planes. Therefore, eigenvalue cases for undriven modes reduce to decoupled 2 × 2 problems in the δr j and δr j variables that are directly analogous to the cases discussed above for the δr ± and δr ± variables. A phase-amplitude analysis is employed in Appendix D to show that in this zero space-charge limit the phase advance σ ej of oscillations in δr j over one lattice period is This limit holds for any periodic focusing functions κ j (s). An immediate consequence of Eq. (55) is that all coordinates connected by linear transformation to the δr j and δr j coordinates must have oscillations with the same phase advance. Normal mode coordinates δR of an envelope mode are related to the δr j and δr j coordinates through linear transformations specified by the eigenvectors as given in Table I. Therefore, the phase advance σ of any normal envelope mode for any transport channel with σ 0j = σ 0 and a symmetric beam with ε x = ε y must satisfy the zero space-charge limit regardless of the details of the focusing lattice. This limit provides a consistency check useful in practical calculations. The result can be understood qualitatively as follows. All j-plane particle oscillations in the zero space-charge limit beam have phase advance σ 0j . Particles contributing to the perturbed beam edge in phase-space must have contributions from positive and negative coordinates, leading to a envelope locus of particle orbits that oscillates at double the particle oscillation frequency. A consequence of the limit (56) and the decoupling of the x− and ycoordinates for zero space-charge is that envelope instabilities for zero space-charge can only occur for σ 0 = π/2 or σ 0 = π. This can be understood as follows. Because the coordinates are decoupled, only lattice resonance instabilities [see Fig. 5(b)] are possible in the x-and y-variables or in any other set of variables connected to these by linear transformations. Therefore, for instability, the eigenvalues must lie on the real axis and lim Q→0 σ must be an integer multiple of π. Equation (56) then establishes a necessary (but not sufficient) condition for instability is that σ 0 be an integer multiple of π/2, of which only π/2 and π are relevant to the analysis in this paper because we consider only stable undepressed particle orbits with σ 0 < π. Another feature of the envelope normal modes is that they evolve with quadratic Courant-Snyder invariants in the normal coordinate phase-space variables of the modes. These invariants are valid regardless of the details of the lattice focusing functions κ j (s) and can be interpreted as phase-space area measures in the normal coordinates of the envelope modes. In Appendix E, a complete derivation of the envelope Courant-Snyder invariants is presented for the case of envelope perturbations in continuous and solenoidal focusing channels with κ x = κ y ≡ κ and a symmetric matched beam with ε x = ε y . In this case the normal coordinates of the envelope mode perturbations are the simple sum and difference coordinates δr ± = (r x ± r y )/2 which evolve according to Eq. (45) with zero driving terms. This situation renders the derivation of the invariants directly analogous to standard treatments for decoupled single-particle orbits and shows that where the w ± are positive mode amplitude functions satisfying For stable envelope oscillations, the amplitudes w ± can be interpreted as scaled maximum excursions of pure envelope modes. Analogous Courant-Snyder invariants exist also for the more complicated case of lattices with κ x = κ y , where the coupled envelope perturbations δr j evolve according to Eq. (37) with zero driving terms. In this case, calculation of normal mode coordinates is much more complicated, and is treated in Appendix E by connecting the formulation to a reference where the general coupling problem is analyzed. Although envelope mode Courant-Snyder invariants increase our formal understanding of the mode structure, in practical applications the utility of these invariants is limited because they require calculation of periodic amplitude functions in the lattice that are not, in general, analytically tractable. However, for a given lattice and beam parameters, these amplitude functions can be numerically calculated and the mode invariants can then be employed to improve understanding of mode evolution and launching conditions in an analogous manner to standard treatments of single-particle dynamics. Finally, the particular solution to Eq. (38) for driving perturbations with δP = 0 is complicated for periodic focusing lattices. In Sec. III we construct the particular solution for the special case of continuous focusing where it can be carried out directly. This special case may provide a partial guide to the much more complicated situation for periodic focusing lattices which is not addressed in this paper. However, future studies of the particular solution will be aided from the analyses of the homogeneous solutions for periodic solenoidal and quadrupole focusing lattices presented in Secs. IV and V, because such solutions can aid in formulation of Green's function type methods to construct the particular solution.

III. ENVELOPE MODES IN CONTINUOUS FOCUSING CHANNELS
The continuous focusing model described in Sec. II B is straightforward to analyze and can serve as a good approximation to more realistic periodic focusing lattices when effects do not depend strongly on the periodic nature of the applied focusing forces. In this regard, we find in Secs. IV and V that the continuous focusing model can provide good estimates of normal mode oscillation frequencies in periodic lattices when system parameters are far from bands of parametric instability. Moreover, the continuous focusing model can be solved exactly for driven envelope perturbations (δP = 0), thereby providing new insight into this previously unexplored situation.

A. Matched envelope solution
In a continuous focusing channel, In this situation, if ε x = ε y = ε, then the matched beam envelope is round with r xm = r ym ≡ r m = const, where r m satisfies the reduced envelope equation of constraint or equivalently, From Eqs. (24) and (60), the depressed phase advance of a particle over an axial distance L p is given by Equations (60) and (61) can be used to relate the space-charge depression σ/σ 0 for the continuous focused beam, which is independent of the specific choice of L p , to the beam parameters as (62)

B. Envelope modes
Using the equilibrium condition in Eq. (61), the linearized equations (45) describing small amplitude elliptical beam perturbations δr ± = (δr x ± δr y )/2 in the continuous focusing model can be expressed in scaled form as where The homogeneous solutions to Eq. (63) are trivially expressible as where r ± (s i ) and r ± (s i ) denote initial (s = s i ) values of the perturbed envelope coordinates. The ± normal mode solutions to the homogeneous equation correspond to the well-known breathing (+) and quadrupole (−) modes with phase advances σ + and σ − [7]. Projections of the modes onto general initial envelope perturbations in x and y are given by Eq. (53) and the pure mode launching conditions are presented in Sec. II F and are illustrated in Fig. 6. The stability of these undriven oscillations is expected because the applied forces of the continuous focusing lattice cannot change the energy of the oscillations. Variations of the normalized phase advances σ ± /σ 0 of the breathing and quadrupole modes calculated from Eq. (64) are plotted in Fig. 7 as a function of σ/σ 0 over the full range of spacecharge intensity 0 ≤ σ/σ 0 ≤ 1. As expected for finite space-charge intensity with σ/σ 0 < 1, the phase advance of the breathing mode is always greater than the phase advance of the quadrupole mode. Also, consistent with Eq. (56), the zero space-charge limit phase advances satisfy lim σ→σ0 σ ± = 2σ 0 . The homogeneous solutions for the normal envelope modes given in Eq. (65) is also straightforward to derive using the matrix formulation in Sec. II F. This derivation is carried out in Appendix F [23]. Results presented in Appendix F provide insight on the matrix formulation and may prove useful as a conceptual guide for calculations with periodic lattices where explicit transformations to decoupled variables are not possible. A Green's Function method [13] based on the the homogeneous solutions (65) can be employed to calculate particular solutions of Eq. (63) in terms of an integral as where are Green's functions for the breathing (+) and quadrupole (−) terms. The particular solutions in Eq. (66) satisfies the initial conditions δr ± (s = s i ) = 0 and δr ± (s = s i ) = 0, and the solution can be readily interpreted as a scaled convolution of the driving perturbations onto the natural harmonic response of the normal modes. This is expected because continuous focusing is spatially homogeneous. Given specific initial mismatches (i.e., δr j and δr j at s = s i for j = x, y) and driving terms (i.e., δκ j , δQ, and δε j ), Eqs. (65) and (66) for the homogeneous (denoted | h here) and particular (denoted | p here) solution components can be applied to calculate the general evolution of the envelope perturbations δr To better understand the range of possible envelope evolutions the driving terms induce, it is instructive to calculate the response to specific classes of driving perturbations that are adiabatic, sudden, ramped, and harmonic. First, in the adiabatic limit, variations of the driving terms δp ± (s) are taken to be slow on the scale of the quadrupole mode wavelength 2πL p /σ − , which is longer than the breathing mode wavelength 2πL p /σ + for σ/σ 0 < 1. In this situation, Eqs. (63) and (64) can be employed to express the particular solutions as The particular solutions (68) can also be derived by partial integration of Eq. (66) and neglecting the terms containing δp ± . Coefficients of the driving terms (δκ x /k 2 β0 ± δκ y /k 2 β0 )/2, δQ/Q, and (δε x /ε ± δε y /ε)/2 in the adiabatic solution (68) (i.e., the terms in the square brackets) are plotted as a function of space-charge depression σ/σ 0 in Fig. 8. These curves describe the relative strength of the response in δr ± to different classes of adiabatic driving perturbations as a function of space-charge intensity.
Next, to analyze the opposite limit of driven perturbations δp ± (s) with sudden changes, we take Here, Θ(x) is the Heaviside step function defined by Θ(x) = 0 for x < 0 and Θ(x) = 1 for x > 0, s = s p is the axial coordinate where the perturbations jump, and  σ /σ 0 a) Adiabatic Solution Coefficients for δr + = (δr x + δr y )/2 b) Adiabatic Solution Coefficients for δr -= (δr x -δr y )/2 are constant amplitudes of the sudden, step-perturbations defined according to Eq. (66) with δQ → δQ = const for the perveance perturbation, etc. For these step-function perturbations, Eq. (66) can be integrated, and the response in the particular solutions for δr ± can be expressed in the same form as the adiabatic solutions with the following substitutions in Eq. (68): The response to other sudden perturbations δp ± (s) can be calculated by superimposing perturbations of the form in Eqs. (69) and (70).
Comparing the adiabatic solutions in Eq. (68) with the sudden solutions given by Eqs. (68) and (71), note that the sudden response is oscillatory with periods 2πL p /σ ± in the solutions for δr ± rather than steady as in the adiabatic case. For the same driving perturbation amplitude [e.g., δQ(s) = δQ for s > ∼ s p for the adiabatic perturbation], the excursions of the envelope oscillations for sudden perturbations oscillate between zero and twice the adiabatic excursion. This characteristic response to a step-function perturbation is illustrated in Fig. 9 for the evolution of δr + /r m resulting from a pure emittance perturbation with ( δε x /ε+ δε y /ε)/2 > 0. The dashed line shows the adiabatic response for a slowly varied perturbation that achieves the same emittance change as the sudden perturbation applied at s = s p {i.e., an adiabatic perturbation with [δε x (s)/ε + δε y (s)/ε]/2 = ( δε x /ε + δε y /ε)/2 for s > ∼ s p }. Another instructive example is the envelope response to linear ramp driving perturbations. In this case we take Here, δp ± is a constant that describes the ramp perturbations applied at axial coordinate s = s p and is defined by replacing δQ → δQ with δQ = const in Eq. (70) 9: Contrast of the envelope response δr+(s) to sudden and adiabatic driving emittance perturbations with ( δεx/ε + δεy/ε)/2 > 0 that achieve the same value of emittance change.
following substitutions in Eq. (68): Comparing the solutions for adiabatic, sudden, and linear ramped driving perturbations given by Eqs. (68)-(73), the response to the ramped perturbations can be interpreted as having both adiabatic and "sudden" oscillatory components given by the first and second terms in Eq. (73), respectively. The particular solution for the ramped perturbation is illustrated schematically in Fig. 10 for the evolution of δr + /r m in response to a pure emittance Finally, we analyze harmonic driving perturbations in δp ± (s) by taking Here, δp + are complex constants defined from δp ± (s) that are of the same form as given in Eq.  (68): For k pure imaginary the particular solutions given by Eqs. (68) and (75) lead to oscillatory solutions with damping (Im k > 0) or growth terms (Im k < 0). For |Re k| sufficiently removed from σ ± /L p , the response is pure oscillatory (Im k = 0), or oscillatory with growth or damping terms (Im k = 0). When |Re k| σ ± /L p , the response has classic linear resonance growth with lim k→σ± /Lp Linear superposition of the specific examples above provides insight into the range of envelope responses to smallamplitude driving perturbations formally given by Eq. (66). Several conclusions can be drawn from these results. First, driving perturbations that are fast on the scale of the breathing and quadrupole envelope mode phase advances σ ± are more dangerous than slow changes because they can result in twice the beam envelope excursion for equal amplitude changes in the driving terms. Second, oscillatory driving perturbations of small amplitude are not problematic unless the driving terms have frequency components with phase advances close to the breathing or quadrupole mode phase advances σ ± . This simple continuous focusing analysis of driven perturbations should help guide the more complicated situation in periodic focusing lattices. However, caution must be employed in application of these results, particularly when making higher-order extrapolations. For example, if a periodic focusing perturbation with δκ x = −δκ y is taken, one might expect to recover net focusing because this perturbation corresponds to a superimposed alternating gradient quadrupole lattice. However, from Eq. (66), such a choice has no projection on the breathing mode to linear order. A recent analytical study of matched beam envelopes in periodic quadrupole focusing channels suggest how such higher-order effects might be more accurately modeled but are beyond the scope of the present analysis [18].

IV. ENVELOPE MODES IN SOLENOIDAL FOCUSING CHANNELS
A beam with an elliptical cross-section envelope and zero total canonical angular momentum {i.e.,  (4), provided that we take the focusing strength to be κ x = κ y ≡ κ = [qB z /(2mγ b β b c)] 2 and interpret the analysis as applying in a local (s-varying) rotating Larmor frame. This correspondence is developed in Appendix B by deriving a self-consistent KV distribution describing elliptical beams with P θ = 0 in a solenoidal focusing channel. The derivation identifies correspondences needed to apply Eq. (4) to this situation including replacements for the laboratory-frame rms emittances ε x and ε y in terms of Larmor-frame invariant emittance measures. In the context of interpreting the envelope equations (4) with appropriate Larmor-frame replacements, the formalism in Sec. II F can be applied to analyze elliptical perturbations in the beam envelope for solenoidal transport channels. For simplicity of presentation, we carry out analysis in this section without making explicit reference to the Larmor frame. However, the Larmor frame transformations presented in Appendix B must be applied to project nonaxisymmetric normal modes (δr x = δr y ) into the stationary laboratory frame to view the evolution there -a distinction not made in previous work on envelope modes in solenoidal transport channels [2,7].
In this section we explore normal mode envelope perturbations (δP = 0) in a periodic solenoidal transport channel with piecewise constant κ(s). Physical solenoidal focusing lattices with fringe fields associated with continuously varying κ(s) can be approximated by the model with piecewise constant κ(s) using the equivalency procedures developed in Appendix A. Mode structures and launching conditions are classified and instabilities are parametrically mapped using the decoupled mode formulation in Sec. II F. In addition to elliptical beam perturbations with zero beam canonical angular momentum (P θ = 0), results for axisymmetric perturbations (δr x = δr y ) can also be applied to beams with nonzero canonical angular momentum using the correspondences derived in Appendix B. Results in this section can also be applied without transformation to idealized beam models where a stationary, partially neutralizing charge distribution produces the specified focusing function κ(s). Driving perturbations with δP = 0 are not analyzed.

A. Matched envelope solution
We consider the case of a periodic, matched solution (r x = r xm and r y = r ym ) to the the envelope equations (4) with equal emittances ε = ε x = ε y and a round beam envelope r xm = r ym = r m . A hard-edge focusing lattice is taken with κ x (s) = κ y (s) ≡ κ(s) piecewise constant as illustrated in Fig. 2(b) with lattice period L p , solenoid strengthκ = const, and occupancy η. The matched envelope solutions are parameterized by L p , σ 0 , σ/σ 0 , and η. Numerical solutions are plotted for one lattice period in Fig. 11 in terms of the scaled envelope radius [r m /( √ 2QL p )] versus normalized axial coordinate (s/L p ) for fixed undepressed single-particle phase advance σ 0 , and values of σ/σ 0 . Solutions for large and small occupancies are contrasted in (a) and (b) and symmetry points are indicated on the plot. Comparing (a) and (b), observe that the radial extent of the periodic flutter of the beam envelope varies strongly with the solenoid occupancy η. This strong flutter variation is expected because the limit η → 1 corresponds to continuous focusing in the Larmor frame with r m = const, whereas the limit η → 0 corresponds to focusing with thin-lens kicks (see Appendix G) and maximal flutter motion. The single-particle phase advances σ 0 and σ are interpreted in the Larmor frame. The analysis in Appendix B shows that the matched, round beam envelope with r xm = r ym = r m described above can be applied to model beams with finite canonical angular momentum P θ = xy −yx ⊥ +[qB z /(2mγ b β b c)] x 2 +y 2 ⊥ in the range −1 < 2P θ /ε < 1 if the emittance ε is replaced by the Larmor invariant emittance according to ε 2 → ε 2 x − 4 r 2 θ ⊥ + 4P 2 θ = const to include the defocusing effect of P θ = 0. The matched beam parameterization employed with L p , σ 0 , σ/σ 0 , and η is convenient to understand envelope stability properties in solenoidal focusing and place results in a form directly comparable to those obtained for quadrupole focusing (see Sec. V). However, in calculating perveance (or current) limits, the periodicity of the lattice plays a less direct role for solenoidal focusing than for quadrupole focusing [17,18]. For a quadrupole lattice alternating focusing and defocusing optics in the period provide net focusing that cannot exist without the regular interchange of focusing and defocusing within the period, whereas in solenoid lattices all optics are focusing and net focusing exists even in the limit of full occupancy (η = 1) where there is no lattice period.

B. Normal modes
Normal mode solutions for elliptical perturbations (δr x = δr y ) about the round, matched envelope solutions described in Sec. IV A are numerically analyzed using the formulation in Sec. II F. For the special case considered of an axisymmetric matched envelope the reduced eigenvalue formulation with decoupled perturbed envelopes δr ± = (δr x ± δr y )/2 is employed to resolve the perturbations into breathing (+) and quadrupole (−) modes. This reduced problem is equivalent to numerically calculating eigenvalues λ n of the coupled 4-dimensional matrix equation (47) to characterize the stability properties of the homogeneous equation (38) with κ x = κ y = κ, ε x = ε y = ε, and zero driving perturbations (δP = 0). Eigenvalues are calculated for specified values of η and ranges of σ 0 and σ/σ 0 to map out mode properties including mode phase advances σ ± and growth factors γ ± associated with broad bands of parametric instability. Results we obtained are displayed in Figs. 12-14. We find that bands of parametric instability with γ + > 1 and γ − > 1 for the breathing and quadrupole modes do not overlap. Therefore, in terms of the symmetry classes presented in Fig. 4 Table I]. The growth factor contours are restricted to σ 0 ∈ (90 • , 180 • ) because no instability is found for σ 0 < 90 • . These plots help clarify regions of parametric envelope instability observed in previous studies [4,7,20]. In the present study, modes are unambiguously identified over the practical range of machine parameters and previous errors in branch calculation are corrected.
To better illustrate the mode structures presented in Figs. 12, in Fig. 13 we plot the mode phase advances and growth factors for fixed σ 0 and η as a function of σ/σ 0 . Curves shown correspond to the vertical dashed lines on the contour plots in Fig. 12. Complex plots of the eigenvalues λ ± relative to the unit circle are superimposed at several locations to further illustrate the eigenvalue structure of the modes and the symmetry class of the instability bands. Individual mode branches for the breathing and quadrupole modes are identified outside of regions of instability using the formulation described in Sec. II F and branches are continuously connected across bands of parametric instability. As expected for finite space-charge (σ/σ 0 < 1), the in-phase oscillations of the breathing mode with δr x = δr y have higher phase advance σ + than the phase advance σ − of the out-of-phase oscillations of the quadruple mode with δr x = −δr y . Continuous model results superimposed are calculated from Eq. (64) using the single-particle phase advances σ 0 and σ derived from the periodically focused, matched beam. Far from instability bands, the continuous formulas provide accurate estimates of the mode phase advances [7]. Curves for the continuous and periodic focused calculations of σ ± almost overlay for σ 0 = 80 • and agreement further improves for lower values of σ 0 . Analysis in Sec. II F [see Eq. (56)] shows that the zero space-charge limit Q → 0 corresponds to stable breathing and quadrupole mode oscillations with These limiting values are indicated on the figure to provide a partial check of results. For another check of results, mode eigenvalues are analytically calculated in Appendix G in the thin-lens limit (η → 0 at fixed σ 0 ) for a beam with full space-charge depression (σ = 0). Phase advances and growth rates calculated from this analysis are indicated in Fig. 13(b). The small discrepancies evident are traceable to sensitivities of results on the solenoid occupancy η.
To better understand sensitivities on the solenoid occupancy, variations of the breathing and quadrupole mode phase advances and growth factors with η are illustrated in Fig. 14. From this plot and Fig. 12, note that instability bands get wider and mode growth rates stronger with decreasing η. Limit points indicated on the plots come from analytical calculations of the band structure in the thin-lens limit (η → 0) at full space-charge depression presented in Appendix G and Eq. (77).
The lattice resonance instability bands have been identified as a half-integer parametric resonance between the envelope mode oscillation frequencies and the lattice [7]. This suggests that the location of the envelope bands can be crudely estimated using resonance conditions based on the breathing and quadrupole mode frequencies σ ± derived from the continuous model with σ + = 2σ 2 0 + 2σ 2 and σ − = σ 2 0 + 3σ 2 [see Eq. (64)] as σ ± = π, or equivalently with σ 0 and σ measured in degrees, for the breathing mode, and   for the quadrupole mode. Comparing to Fig. 12, these resonance formulas roughly predict band locations, but correspond more closely to the lower thresholds in σ 0 than the maximum growth factor locations as might be anticipated. Moreover, the formulas cannot capture the broad parametric nature of the instability or changes in band structure with solenoid occupancy η. In order to provide better guidance for the band locations in practical applications, parametric data was employed to calculate curve fits for the instability thresholds. We find that instability threshold data can be roughly fit by curves in σ 0 and σ corresponding to: centered elliptical curves with for the left-and right-boundary curves of the breathing mode band, a curve with form for the left-boundary of the quadrupole mode band, and a line with where average and maximum deviations between the threshold data and the fit formula are ∼ 3 • and ∼ 6 • for the left-edge and ∼ 4 • and ∼ 6 • for the right-edge. Similarly, using the threshold data over the same range for the quadrupole band with the formulas in Eqs. (81) and (82), we obtain the fits g(η) = 1, left-edge, 0.227 + 0.173η, right-edge.
where average and maximum deviations between the threshold data and the fit formulas are ∼ 5 • and ∼ 8 • for the left-edge and ∼ 2 • and ∼ 3 • for the right-edge. Here, we have neglected weak η and σ 0 variation in the formula for the left-edge of the band and we have neglected variation in σ 0 on the right-edge to produce formulas of similar accuracy to the breathing band. Methods employed by Lee and Briggs [17] to analytically estimate breathing mode band widths for a focusing function κ(s) with a constant and sinusoidally varying terms may be generalizable using Fourier expansions of both the piecewise constant lattice focusing function κ(s) employed here and the matched envelope to calculate thresholds of the decoupled breathing and quadrupole bands. Future research on this topic is desirable to obtain more reliable band threshold estimates and to better understand the structure of the envelope instabilities.
As discussed in Sec. II F, the breathing and quadrupole modes are decoupled, and consequently pure modes can be simply launched at axial coordinate s = s i by taking initial conditions as given in Eqs. (53) and (54). This launching condition is illustrated in Fig. 6 Table I for stability and instability, respectively. For the case of instability, Table I and Figs. 4 and 5 show that certain subsets of initial conditions with either breathing or quadrupole symmetry can project onto pure exponentially damping or growing modes. Mode phase-space invariant (see Sec. II F and Appendix E) can be employed analogously to Courant-Snyder invariants in single-particle dynamics to efficiently launch stable breathing or quadrupole modes with specific amplitudes and phases and calculate maximum mismatch excursions at any location in the lattice period.

and corresponds to cases (a) and (c) in
As a practical matter, note in the parametric plots of the bands of mode instability presented in Fig. 12 that the breathing mode band is encountered before the quadrupole band when approaching from lower values of σ 0 . Therefore the breathing mode instability is generally of greater importance in the selection of machine operating points. Finally, it should be stressed that the breathing mode results presented can be interpreted directly as being laboratory frame, but the quadrupole mode rotates according to the local Larmor frame transformations presented in Appendix B. Breathing mode stability results presented are also valid for beams with finite canonical angular momentum P θ = 0, whereas quadrupole mode results are valid only for P θ = 0.

V. ENVELOPE MODES IN QUADRUPOLE DOUBLET FOCUSING CHANNELS
In this section, we explore normal mode envelope perturbations (δP = 0) in alternating gradient quadrupole focusing channels. Driving perturbations with δP = 0 are not analyzed. Results are presented in analogous manner to Sec. IV for solenoidal focusing. However, in contrast to Sec. IV, all results are laboratory frame and can be interpreted directly without transformation. The beam considered in this section has zero total canonical angular momentum P θ = xy − yx ⊥ = 0 and the focusing channel has no skew couplings so that the envelope equations (4) apply (see Appendix B). Piecewise constant focusing strength κ x (s) = −κ y (s) ≡ κ q (s) is taken and results can be applied to realistic focusing channels with continuously varying κ q (s) using the equivalency procedures developed in Appendix A.

A. Matched envelope solution
We consider the case of a periodic, matched solution (r j = r jm for j = x, y) to the the envelope equations (4) with equal emittances ε x = ε y = ε. A hard-edge focusing lattice is assumed with piecewise constant κ x (s) = −κ y (s) ≡ κ q (s) as illustrated in Fig. 2(c) with lattice period L p , quadrupole strength κ q = const, occupancy η, and syncopation factor α. The matched envelope solutions are parameterized by L p , σ 0 , σ/σ 0 , η, and α. Numerical solutions are plotted for one lattice period in Fig. 15 in terms of scaled envelope radii [r j /( √ 2QL p )] versus normalized axial coordinate (s/L p ) for α = 1/2 [FODO channel, (a)] and α = 0.1 [strong syncopation, (b)]. In these plots, the undepressed singleparticle phase advance σ 0 is fixed, and the solution is shown for values of σ/σ 0 . Approximate, analytical solutions can also be constructed using the perturbative formulation of Lee [18]. Symmetry points are indicated on the plot. Contrasting (a) and (b), observe that syncopation with α = 1/2 reduces the symmetry of the matched envelope but that maximum excursions change little. In contrast to solenoidal transport channels, at fixed σ 0 , little variation in matched envelope structure is observed with changes in quadrupole occupancy η. The value of occupancy employed η = 0.6949 corresponds to the hard-edge equivalent value of a FODO electric quadrupole lattice in the High Current Transport Experiment (HCX) at the Lawrence Berkeley National Laboratory [24,25]. This equivalence prescription is derived as an example application in the general discussion of fringe field equivalent models presented in Appendix A. Negligible difference in the structure of the matched beam envelope is observed between the actual focusing function κ q (s) derived from the linear applied field components of the physical lattice and the hard-edge equivalent model with piecewise constant κ q (s).

B. Normal modes
Normal mode solutions are numerically analyzed using the formulation in Sec. II F. For the matched beam solutions described in Sec. V A, mode perturbations are not simply decoupled by sum (+) and difference (−) coordinates δr ± = (r x ± r y )/2 as in the case of continuous and solenoidal focusing channels analyzed in Secs. III and IV, and we employ the full coupled eigenvalue formulation based on Eq. (38) with κ x = −κ y = κ q , ε x = ε y = ε, and zero driving perturbations (δP = 0). Eigenvalues are calculated for specified values of η and α and ranges of σ 0 and σ/σ 0 to map out mode properties. Results we obtained are displayed in Figs. 16 are not decoupled in quadrupole transport channels, we find that the pure modes have properties analogous to the sum (i.e, "breathing") and difference (i.e., "quadrupole") modes studied in continuous and solenoidal channels. We label envelope modes in the quadrupole focusing channel that in stable regions degenerate into quadrupole-and breathing-like perturbations with B and Q, respectively, (e.g. σ with = B, Q for the mode phase advances) to avoid confusion with the purely decoupled sum and difference modes analyzed in Secs. III and IV. The "breathing" and "quadrupole" mode phase advances σ B and σ Q and growth factors γ B and γ Q are contoured in Fig. 16 as a function of σ 0 and σ/σ 0 for σ 0 ∈ (0, 180 • ) and σ/σ 0 ∈ (0, 1]. Results are shown for quadrupole occupancy η = 0.6949 and syncopation factors α = 1/2 (FODO lattice) and α = 0.1 in columns (a) and (b) to contrast a symmetrical FODO lattice with the case of strong syncopation. Growth factors of both the breathing (γ B > 1) and quadrupole (γ Q > 1) modes are superimposed on the same plots because the bands do not overlap, and growth factors of the corresponding damped modes are not plotted because damping factors are given by 1/γ of the growing mode [see Fig. 4 and Table I]. No instability bands are found for σ 0 < 90 • so the growth contours are restricted to σ 0 ∈ (90 • , 180 • ). In terms of the symmetry classes presented in Fig. 4, eigenvalues fall into class (a) when there is stability and class (b) [confluent resonance] or (c) [lattice resonance] when there is instability. The confluent resonance involves locking of the breathing and quadrupole modes and occurs for both FODO and syncopated (α = 1/2) lattices, while the lattice resonance is observed only for syncopated lattices (see the thin band to the left of confluent resonance band on the right column of the figure) which have a lesser degree of symmetry. The lattice resonance also only appears for the breathing mode. As for solenoidal transport, the breathing mode instabilities appear to be more dangerous for practical, intense beam transport lines because it is encountered first when increasing σ 0 from small values. These plots help clarify regions of parametric envelope instability observed in previous studies of FODO lattices [4,21] and a syncopated quadrupole lattice [7]. The thin band associated with the lattice resonance for syncopated quadrupole lattices was previously overlooked. Also, mode branches are unambiguously identified, errors in branch calculations are corrected, and parametric sensitivities and mode properties are throughly characterized over the range of practical machine parameters.
To better illustrate the mode structures presented in Fig. 16, in Fig. 17 we plot the mode phase advances and  growth factors for fixed σ 0 , α, and η as a function of σ/σ 0 . Curves shown correspond to the vertical dashed lines on the contour plots in Fig. 16. Complex plots of the eigenvalues λ n relative to the unit circle are superimposed at several locations to further illustrate the eigenvalue structure of the modes. Individual mode branches for the breathing and quadrupole modes are identified outside of regions of instability using the formulation described in Sec. II F and branches are continuously connected across bands of parametric instability. The phase advance of the breathing mode σ B is larger than the quadrupole mode σ Q for finite space-charge (σ/σ 0 < 1). Continuous focusing model predictions shown for the sum (i.e., breathing) and difference (i.e., quadrupole) mode phase advances σ + and σ − are derived from Eq. (64) using the single-particle phase advances σ 0 and σ calculated from the periodically focused, matched beam. Far from bands of instability, σ + and σ − provide accurate estimates for the numerically calculated breathing and quadrupole mode phase advances σ B and σ Q [7]. This correspondence helps justify the identification of the mode branches being breathing-and quadrupole-like. The continuous focusing and alternating gradient model phase advances are in better agreement for lower values of σ 0 and weaker syncopation (α → 1/2). Analysis in Sec. II F [see Eq. (56)] shows that the zero space-charge limit Q → 0 corresponds to stable breathing and quadrupole mode oscillations with   For FODO lattices with the same values of the particle phase advances σ 0 and σ, little variation is observed in envelope mode structure with changes in quadrupole occupancy η. On the other hand, significant variations can occur in lattices with strong syncopation (α far from 1/2). To illustrate this point, the mode phase advances σ B and σ Q and growth factors γ B and γ Q are plotted in Fig. 18 for several values of η. Results are contrasted in columns (a) and (b) for a FODO lattice (α = 1/2) and a strongly syncopated lattice with α = 0.1. The small variations presented in Fig. 18(a) illustrate the insensitivity of the mode structure to changes in η for a FODO lattice. The location and strength of the confluent resonance instability band changes only weakly with η. Even smaller variations are observed in FODO lattices when making comparisons for σ 0 farther from the first instability band or when comparing results for the hard-edge lattice taken here with full fringe field models calculated from physical lattices (see Appendix A). Results from eigenvalues analytically calculated in Appendix G for the thin-lens limit (η → 0 at fixed σ 0 ) of a FODO lattice with a fully space-charge depressed beam (σ = 0) are indicated in Fig. 18(a). The good agreement with the thin-lens results is traceable to mode structure in a FODO lattice being relatively insensitive to the solenoid occupancy η. For FODO lattices, σ 0 and σ are the most important parameters, whereas shape parameters characterizing the occupancy and fringe components of the focusing field are of much lesser importance. In contrast, for strong syncopation, it is evident from Fig. 18(b) that the mode structure can vary significantly with occupancy η. This strong variation is correlated to the existence of the second band of parametric instability (lattice resonance) for the breathing mode, which was shown in Sec. IV to vary strongly with η in a solenoidal focusing lattice. Observe from the α = 0.1 plots in Fig. 18(b) that the instability band of the lattice resonance is weakly expressed compared to the confluent resonance band for higher lattice occupancies η, whereas it is strongly expressed for η small. High resolution numerical studies suggest that the lattice resonance instability band disappears (becoming both thinner and weaker) as α → 1/2. The confluent resonance instability band has been identified as a half-integer parametric resonance between both envelope mode oscillation frequencies and the lattice, and the lattice resonance instability band has been identified as a half-integer parametric resonance between one envelope mode and the lattice [7]. This suggests that the location of the envelope bands can be crudely estimated using resonance conditions based on the breathing and quadrupole mode frequencies σ ± derived from the continuous model with σ + = 2σ 2 0 + 2σ 2 and σ − = σ 2 0 + 3σ 2 [see Eq. (64)] as (σ + + σ − )/2 = π for the confluent resonance band and σ + = π for the lattice resonance band, or equivalently with σ 0 and σ measured in degrees, for the confluent resonance band (both breathing and quadrupole modes), and for the lattice resonance band (breathing mode). Comparing to Fig. 16, Eq. (86) roughly predicts the maximum growth factor in the confluent resonance band for the breathing and quadrupole modes, and Eq. (87) roughly predicts the location of the thin lattice resonance band associated with the breathing mode in syncopated lattices. However, the formulas cannot capture broad parametric nature of the instability or changes in band structure with solenoid occupancy η or syncopation factor α. Moreover, Eq. (87) fails to suggest the non-existence of the band for FODO lattices with α = 1/2. In order to provide some guidance for band locations in practical applications, parametric data was employed to calculate curve fits for the instability thresholds for a FODO lattice. A least-squares fit was carried out using a large amount of threshold data for the confluent resonance band. Elliptic and linear boundary constraint equations in σ 0 and σ measured in degrees were applied with forms: for the left-edge of the FODO confluent resonance band, and for the right-edge of the FODO confluent resonance band. Here, f (η) and g(η) are functions to be determined by fitting the data, and the forms taken correspond to the most general centered ellipse and line containing the limit point σ = σ 0 = 90 • . We find very little variation in η when fitting to these forms (almost to the accuracy of the numerical data) and errors are minimized by taking Maximum errors (all η values) in threshold prediction from Eqs. (88)-(90) are ∼ 5 • and ∼ 2 • for the left-and rightedge of the band, respectively. No two-dimensional curve fits were attempted to understand and characterize the more intricate band changes observed in a syncopated quadrupole lattice with α = 1/2. Generalizations of methods employed by Lee and Briggs [17] using Fourier expansions of the lattice focusing function κ q (s) might be employed to derive analytical formulas to better understand the instability thresholds. This topic is left for future research. To further illustrate the mode structure, pure mode launching conditions in δR = (δr x , δr x , δr y , δr y ) are plotted in Fig. 19 for both the breathing and quadrupole modes. Launching conditions are plotted at axial locations corresponding to matched beam symmetry points in a FODO lattice: the axial mid-drift with focusing-and defocusing-in-x quadrupoles at smaller and larger s values, and mid-quadrupole within a focusing-in-x quadrupole. Launching conditions at the other axial mid-drift and mid-defocusing-quadrupole locations are simply derivable by symmetry from the plots shown. The launching conditions employed are described in Table I  Lp 0 ds √ r xm r ym } and launching conditions are plotted versus mode phase ψ for −π ≤ ψ ≤ π (the absolute value of the phases are arbitrary). Parameters chosen for Fig. 19 correspond to stable modes with eigenvalues and launching conditions falling into class (a) in Fig. 4 and Table I. To launch a stable mode with a prescribed maximal excursion in δr j (s), it suffices to select a particular axial location s = s i of interest, calculate needed eigenvectors for the full lattice period as a function of s, calculate the maximum value of |δr j (s i )| with respect to both s and the mode phase ψ , and then scale the mode amplitude A such that this maximum value equals the prescribed excursion. Outside of exceptional cases where the mode phase advance σ is a rational fraction of 2π (i.e., the lattice phase advance), this maximal excursion will be approached arbitrarily close as the beam propagates through the periodic lattice.
Comparing the launching conditions in Fig. 19, we observe that initial conditions in δr j and δr j for pure quadrupole or breathing mode launches vary significantly in form with the axial location within the lattice period. This contrasts the situation for the decoupled sum and difference modes analyzed in Secs. III and IV for continuous and solenoidal   focusing channels. In these decoupled cases the launching conditions are the same form independent of axial location [see Eqs. 53 and (54)] and, to launch a given mode amplitude, only the relative projections on the normal mode coordinates and angles vary with axial location in the lattice period. The decoupled mode launch conditions in Eqs. (53) and (54) are a poor approximation to the pure breathing and quadrupole modes launching conditions for an alternating gradient channel even though the modes have rough breathing and quadruple symmetry for specific launch phases and axial coordinate locations. Thus the use of the continuous model as a guide for pure mode launchings in alternating gradient focusing lattices (common in simulations studies etc.) can result in significant launching errors with unwanted projections on the mode desired to be suppressed. A useful reduced prescription for stable pure mode launching (σ 0 < 90 • ) at the mid axial location of a focusing-in-x quadrupole in a FODO lattice (where the x excursion of r xm is maximum and the y excursion in r ym is minimum) is to take δr x to be a specified fraction of δr y and zero perturbed angles (δr j = 0). By symmetry, analogous launching conditions with the same fractions apply at the mid axial location of a focusing-in-y quadrupole. The numerically calculated excursion ratio δr x /δr y needed for this pure mode launch is illustrated in Fig. 20. The excursion ratio is plotted separately for pure breathing and quadrupole mode launches as a function of σ 0 for several values of quadrupole occupancy η and σ/σ 0 . Observe that the excursion ratio varies weakly in both occupancy η and space-charge depression σ/σ 0 relative to variations in σ 0 . For approximate launches these variations in η and σ/σ 0 can sometimes be neglected and the curves shown can then be applied for values of these parameters not shown. Observe that for decreasing values of σ 0 the excursion ratio approaches unity. This is expected because alternating gradient focusing is better approximated by continuous focusing for lower values of σ 0 [3]. The launching conditions illustrated in Fig. 20 hold only for a lattice with FODO symmetry (α = 1/2) and are violated if the lattice is syncopated with α = 1/2.

VI. CONCLUSIONS
This study extends the present understanding of the stability properties of mismatch modes supported by the KV envelope equations describing the transverse evolution of a coasting, unbunched beam envelope. New, generalized linear perturbation equations were derived that describe mismatch evolution resulting from both the usual mismatch source of errors in envelope coordinates from the periodic matched-beam condition as well as driving errors in focusing, perveance, and emittance. Important general properties of undriven envelope modes were analyzed including symmetry classes, zero space-charge limits, and phase-space invariants. The linear perturbation equations were solved in general for continuous focusing channels, and example solutions presented illustrate important differences in the action of adiabatic and sudden driving errors as well as the scaling of contributions to the breathing and quadrupole mismatch modes from the driving errors. The normal modes in the absence of driving errors were also parametrically analyzed in hard-edge, periodic solenoidal and syncopated, alternating gradient quadrupole focusing channels. This analysis extended earlier work by Struckmeier and Reiser [7] by systematically identifying possible envelope instabilities over a wide range of system parameters, thereby providing valuable data for practical applications. For solenoidal focusing, a systematic derivation of the consistent KV distribution identified invariant emittances and the appropriate Larmor frame transformation to interpret the analysis under -both of which are necessary for elliptical beam perturbations and were overlooked in earlier studies. It was shown that the solenoid occupancy strongly influences stability properties and that "quadrupole" symmetry envelope perturbations are generally less critical than "breathing" symmetry envelope perturbations for intense beam transport, thereby simplifying design considerations. For a quadrupole focusing, a previously overlooked mode was found associated with syncopated doublet lattices where one drift between quadrupoles is shorter than the other. For both solenoidal and quadrupole focusing channels, pure mismatch mode launching conditions were constructed. Exact launching conditions were specified for pure breathing or quadrupole envelope mode evolutions in solenoidal channels and approximate launch conditions were formulated for quadrupole focusing channels at lattice symmetry points. These new results provide insight on beam envelope control and have important applications in envelope matching including possible halo mitigation through breathing mode suppression, and in simulation studies.
The present study is limited in that only specific (though representative) cases of linear solenoidal and quadrupole focusing lattices have been analyzed. Also, analysis of properties of mode evolutions due to driving term errors in focusing, perveance, and emittance in periodic focusing channels is left for future study.

APPENDIX A: FRINGE FIELD EQUIVALENT MODELS FOR APPLIED FOCUSING ELEMENTS
Physical electric or magnetic focusing elements in the laboratory have material structures that produce applied fields consistent with the three-dimensional Maxwell equations. The resulting fields will, in general, contain both linear components that contribute to the applied focusing functions κ j (s) (see Sec. II B) as well as nonlinear error terms outside the scope of this paper. Typically, the material extent of a focusing element is axially long relative to the transverse extent of it's clear-bore aperture. Near the axial mid-plane, the elements in this case can often be regarded as two-dimensional axially extruded structures, and the fields within this region can then be well-approximated by lowerdimensional models. Near the ends, this lower-dimensional approximation breaks down, and there are significant field variations as well as unwanted nonlinear field terms. In many cases the ends have approximate left-right symmetries about the axial mid-plane of the element. Consistent with these field variations in the element, the linear focusing functions κ j (s) of the lattice will generally vary in s: typically with near constant value in the neighborhood of the axial mid-plane till the end regions are approached, then smoothly dropping to exponentially small values outside of the ends. This s-variation near the ends is called the fringe field and is associated with all physical focusing elements -leading to the obvious question on how well so-called hard-edge lattice models with piecewise constant κ j (s) such as introduced in Sec. II B apply to physical lattices.
Differences between physical and hard-edge focusing functions κ j (s) can be regarded as driving perturbations [i.e., the δκ j (s) in Sec. II F] about a hard-edge lattice. Alternatively, the physical and hard-edge focusing functions can be thought of as producing different matched beam envelopes for the same beam parameters with deviations between these matched envelopes leading to possible variations in envelope mode structure about the matched beam. In either case, prescriptions must be given for best possible choices of "equivalent" hard-edge focusing parameters. In this Appendix we take the latter interpretation of the physical and model (hard-edge or otherwise) focusing functions giving different matched envelopes and develop procedures for fixing "equivalent" model focusing functions to the physical lattice that minimize deviations between physical and model matched beam envelopes. No fully general optimal procedure can be given, because deviations depend on the specific lattice considered and the structure of the matched envelope. With regards to low-order beam evolutions modeled by the KV envelope equations (4), we generally find that hard-edge models produce results that are close to physical models for a periodic lattice if appropriate, "equivalent" hard-edge equivalent parameters (e.g., for solenoids the occupancy η and focusing strengthκ) are employed. The level of agreement improves for simpler lattices with higher degrees of symmetry [e.g., the solenoidal and quadrupole hardedge lattices introduced in Sec. II B]. This agreement is likely a consequence of compensating errors and that periodic focusing "kicks" in a repeating solution can be well modeled by kicks of varying shape that impart the same total impulse. On the other hand for more general lattices and evolutions (e.g., a matching section where beam envelope parameters are transformed and there is no periodicity), uncompensated accumulated phase errors in the envelope excursions can result in significant errors if accurate fringe-field models are not employed.
No unique prescription exists on how to best fix parameters employed to model κ j (s) for either hard-edge models or for more detailed fringe-field models. In both cases we refer best fit model parameters as "equivalent fringe" parameters. Here we outline several reasonable procedures that can be used to fix equivalent fringe parameters.
One way to fix equivalent fringe parameters is through the dynamics of the particles. Equating matrix elements of the 2 × 2 transfer maps M j (s|s i ) [see Eq. (19)] connecting one side of the optic before the fringe entering the focusing element to the end of the fringe leaving the focusing element will result in equivalent evolutions in drift regions. It is usually sufficient to carry out such procedures only for the undepressed orbits (Q = 0 maps) which provide a direct measure of the applied focusing functions. This will result in three, rather than four, independent constraints in each of the x-and y-planes (the x-and y-planes can yield identical or related constraints for some optics) because M j (s|s i ) is symplectic and has a unit determinant. Unfortunately, these kinematic methods of fixing equivalent fringe parameters prove less than satisfactory in practical applications because they do not result in scale invariant models and must be repeated with changes in particle energy. Moreover, if there are more equivalent fringe parameters than kinematic constraints employed, the method is ambiguous and must be supplemented by additional constraints.
Because of these limitations, we employ an alternative, non-kinematic method of fixing equivalent fringe field parameters defined as follows. Denote the physical and equivalent fringe model focusing functions by κ j (s) and κ f j (s). Then we specify κ f j (s) by requiring: (1) Equal focusing strength at the mid-planes of the elements, i.e., where s = {s m } are the axial coordinates of the mid-planes of the physical focusing elements.
(2) Correct integrated focusing strength in each focusing element, i.e., where the integration range contains the full fringe field of the physical and equivalent elements.
(3) Minimum least square error between the physical and fringe equivalent lattices. Here, we take the error measure to be given by The constraint (1) allows easy understanding of the overall "excitation" needed for the focusing element. The physical element focusing strength κ j (s m ) can also be replaced in (1) by lower-dimensional optical models (e.g. a 2D analytical quadrupole model) to allow analytical estimates of needed excitations. The constraint (2) maintains a good low-order equivalency of focusing strength. To obtain higher-precision equivalences, (2) can be replaced by a constraint for equal values of the single particle phase advance σ 0j in the physical and equivalent lattices to obtain partial kinematic equivalency. Finally, (3) constrains any additional parameters after (1) and (2) establish low-order equivalency. Note that κ f j (s) = κ j (s) corresponds to zero error and an identical kinematic transfer matrix M j (s|s i ) between all points in the lattice with or without space-charge. Constraints (2) and (3) both recover the correct thin-lens limit. Additional criteria can also be added to the error function in (3) to further constrain the equivalent fringe model for minimum envelope excursions and other criteria [26]. For the special case of solenoidal and FODO quadrupole lattices with piecewise constant κ j (s) as defined in Sec. II B, only two equivalent lattice parameters need be set (e.g.,κ and η for the solenoidal lattice) for which the first two constraints are sufficient and the third error constraint is not needed. We illustrate the results of this equivalent fringe procedure with the FODO electric quadrupole lattice of the High Current Transport Experiment (HCX) at the Lawrence Berkeley National Laboratory [24,25]. As shown in crosssection in Fig. 21(a), the quadrupoles are formed from left and right conducting end-plates held at equal magnitude positive and negative biases and conducting rods attached to each end-plate along the principal x-and y-axes as indicated. The structure is mechanically supported at the axial mid-plane by a plate connected to the end-plates by insulators, end-plates are 12.7mm thick with clear-bore aperture holes of radii r p = 23.0mm, rods are of radii ∼ (8/7)r p = 26.3mm (ratio chosen to minimize nonlinear fields in 2D axial extruded geometry of the mid-plane) and axial length 154.5mm, and gaps between rods and end-plates are 17.7mm. The quadrupoles are arraigned with an overall lattice period of L p = 435.2mm and adjacent end-plates are biased to the same potential. This detailed lattice structure was loaded into a 3D electrostatic field solver in the WARP code [27] and the linear applied field moments were extracted from the full 3D electric field solution. These applied field moments can then be linearly scaled for the quadrupole bias (excitation) employed and then Eq. (11) is used to calculate the physical lattice focusing functions κ x (s) = −κ y (s) ≡ κ q (s) shown in Fig. 21. Two equivalent fringe focusing functions derived by the procedure above are also shown in Fig. 21 superimposed with the physical focusing function. The equivalent model in Fig. 21(b) is the hard-edge model in Sec. II B with square-edges [i.e., piecewise constant κ q (s)], and the equivalent model in Fig. 21(c) has κ q (s) with constant value near the mid-plane and linear ramp variations in s (zero to mid-plane value) entering and exiting the quadrupole. Only one quadrupole of the symmetric FODO lattice is shown for clarity. The procedure obtains for the square-edge model a best fit occupancy of η = 0.6949 (corresponding to axial length ηL p /2 = 151.5mm), and for the ramped-edge model a central flat region of axial length 117.1mm with ramped ends of length 34.5mm. The total error according to criteria (3) above is about 74 times larger for the square-edge model than for the ramped-edge model. Matched beam envelopes for both the physical focusing function and the square-edge and linear-ramp equivalent fringe models shown in Fig. 21 are contrasted in Fig. 22. In this procedure the physical quadrupole excitation was chosen such that σ 0 = 80 • . For both solutions values of σ 0 and σ/σ 0 are listed for the same values of perveance and emittance corresponding to a usual operating point in the HCX experiment. Little difference in matched envelopes are observed between the physical and either equivalent fringe models with the linear ramped equivalent model giving smaller deviations than the square-edge model (matched envelopes almost overlay). Deviations between the matched envelope structure of the physical and equivalent lattices are even smaller than those shown in the figure if the excitations of the physical and equivalent lattices are each separately tuned to obtain the same values of σ 0 after the equivalency procedure is applied to set geometric parameters of the equivalent lattice. These results illustrate the general feature that simple equivalent fringe models result in little change in matched envelope structure for lattices with a high degree of symmetry. In situations where there is an s-varying solenoidal magnetic field B(s) = −[B z (s)/2](xx + yŷ) + B z (s)ẑ in addition to (or in place of) other physical laboratory-frame focusing forces described by κ x (s) and κ y (s) (e.g., electric or magnetic quadrupoles or a uniform, partially neutralizing background), the transverse equations of motion of a beam particle can be expressed as [3,28] x Here, is the Larmor wavenumber, ω c = qB z /m is the cyclotron frequency, and −∂φ/∂x and −∂φ/∂y are the x-and ycomponents of the self-electric field of the beam in the electrostatic approximation. For solenoidal focusing κ x = κ y = 0 and and (q/mγ 3 b β 2 b c 2 )∂φ/∂y = −Qy/[(r x + r y )r y ], corresponding to the interior fields of a uniform density elliptical beam in free-space with principal axis radii r x and r y oriented along the x and y coordinate axes.
For solenoidal focusing, the cross-coupled form of Eq. (B1) results in a macroscopic rotation of the beam about the longitudinal axis. Because of this, we analyze the particle equations in ax −ỹ coordinate system rotated by an anglẽ ψ(s) relative to the x − y laboratory-frame as indicated in Fig. 23 with x =x cosψ + y sinψ, Henceforth, we use tildes (˜) to denote quantities measured in the rotating-frame. For a uniform density elliptical beam with (possibly evolving) principal axesrx = 2 x 2 ⊥ andrỹ = 2 ỹ 2 ⊥ along the rotatingx andỹ coordinate axes (see Fig. 23), the self-field forces within the beam are given by Paralleling the analysis of Wiedemann [28], we find thatψ can be selected such that the principal axes of the uniform density elliptical beam remain aligned with thex-andỹ-axes as the elliptical beam evolves. This occurs whenψ satisfies the "Larmor-frame" conditionψ = −k L with the initial conditionψ(s = s i ) = 0, or equivalently that For this choice ofψ, the particle equations of motion (B1) can be simply expressed in the rotating "Larmor-frame" as Consistent with this Larmor transform, the particle angles x and y transform as x =x cosψ + y sinψ + k L (x sinψ − y cosψ), y = − x sinψ + y cosψ + k L (x cosψ + y sinψ). (B7) The particle equations of motion (B6) in the rotating Larmor-frame variables are the same form as the laboratoryframe equations of motion (1) with κ x = κ y → k 2 L . Therefore, the usual KV analysis [29] involving transformations to phase-space invariants that are quadratic forms in conjugate phase-space variables can be employed in the rotatingframe to construct a self-consistent, KV equilibrium solution to the Vlasov equation. In the Larmor-frame, the Vlasov equation for the single-particle distribution functionf(x,x ,ỹ,ỹ ) consistent with the particle equations of motion (B6) can be expressed as The KV equilibrium distributionf satisfying the Vlasov equation can be written as where δ(x) is the Dirac delta-function,εx andεỹ are conserved rms emittances of the uniform density beam in the rotating-frame defined byεx and the envelope radiirx = 2 x ⊥ andrỹ = 2 ỹ ⊥ obey the envelope equations If the˜s are dropped, Eq. (B9) can also be applied to describe the laboratory-frame KV distribution appropriate for quadrupole or continuous focusing channels that is consistent with ε x = const, ε y = const, and the envelope equation (4) [1,3] for the evolution of r x = 2 x 2 ⊥ and r y = 2 y 2 ⊥ . The laboratory-frame envelope equations (4) are identical to the Larmor-frame envelope equations (B11) under the substitutions Therefore, we can apply the laboratory-frame envelope equations (4) to analyze the evolution of an elliptical beam envelope in a solenoidally focused transport channel provided we take κ x = κ y = k 2 L and interpret the results as applying in the rotating Larmor-frame defined by Eqs. (B3) and (B5). This correspondence also requires that the stationary laboratory-frame x-y coordinate system be chosen to lie along the principal axes of the initial (s = s i ) elliptical perturbations with initial values r x (s i ) = 2 x 2 (s i ) ⊥ =rx(s i ) and r y (s i ) = 2 y 2 (s i ) ⊥ =rỹ(s i ) and the particle phase advances (i.e., σ 0 , σ x and σ y ) must also be interpreted as being defined in the rotating Larmor-frame.
Moment calculations verify that the KV distribution (B9) is self-consistent and help elucidate the structure of the rotating equilibrium. The Larmor-frame transformation defined by Eqs. (B3) and (B7) is canonical and phase-space area preserving, with dxdy = dxdỹ and dx dy = dx dỹ and hence the distribution transforms to the lab-frame as f (x,x ,ỹ,ỹ ) = f (x, x , y, y ). (B13) The rotating-frame densityñ = d 2x f is equal to the lab frame density, i.e., n(x,ỹ) = n(x, y) = d 2 x f (B14) and can be calculated asñ verifying consistency. Laboratory-frame statistical averages are simply expressed in terms of rotating-frame distribution as Various moments of the Larmor-frame KV distribution (B9) are summarized in Table II. The moments in Table II can also be applied to the laboratory-frame KV distribution for quadrupole and continuous focusing channels by simply dropping all˜s. It should be stressed that the rotating-frame emittancesεx andεỹ are not, in general, equal to the laboratory-frame emittances ε x and ε y , which evolve in s for an elliptical KV beam in a solenoidal focusing channel. Laboratory-frame projections of the moments are straightforward to calculate from the transformations (B3) and (B7). The angular momentum of the KV distribution (B9) is important in understanding the structure of the equilibrium. The vector potential of the axial magnetic field can be expressed as A = (B z /2)(yx − xŷ), and canonical momenta conjugate to the coordinates Q x = x and Q y = y are P x = x − k L y and P y = y + k L x. Therefore, the single-particle canonical angular momentum p θ ≡ Q×P ·ẑ can be expressed as For any distribution of particles, it follows from the equations of motion (B1) that the system canonical angular For a solenoidal focused beam, κ x = κ y = 0 and the first moment term in Eq. (B19) vanishes. Using Green's functions for the Poisson equation ∇ 2 φ ⊥ = −qn/ 0 , it can be shown that the second moment term in Eq. (B19) also vanishes provided any conducting boundaries are azimuthally invariant (∂/∂θ = 0). Thus, in a solenoidal channel, P θ = const for any distribution function focused within an axisymmetric beam pipe. For the specific KV distribution in Eq. (B9), we calculate This zero value of P θ is appropriate for a beam launched from a field-free source (i.e., B z = 0 at the source location), and results in maximum focusability for given Larmor-frame emittancesεx andεỹ. Beams injected from a source immersed in a finite axial magnetic field B z have nonzero P θ and higher effective emittances [30,31]. Also, for a KV equilibrium evolving in an alternating gradient focusing channel, Eq. (B9) can be applied in the laboratory frame (drop all˜s) to show that P θ = xy − yx ⊥ = 0 in this case too. Elliptical KV distributions analogous to Eq. (B9) can also be constructed for nonzero values of P θ = const for both solenoidally focused beams in the Larmor-frame and quadrupole focused (possibly with skew-couplings) beams in the laboratory-frame. However, these distributions will result in envelope equations that differ in form from Eq. (4) and have more complicated invariant emittances explored by Dragt, Neri, and Rangarajan [32] and Barnard (see Ref. [33]) because the elliptical beam axes rotate within the Larmor-frame. Coupled moment formulations have been employed to analyze the special case of elliptical envelope perturbations about a round-beam equilibrium for a continuously focused beam with P θ = 0 [34]. In the present paper, we explore only cases of envelope evolution that we can readily map to Eq. (4) while concretely connecting the model to physically relevant focusing systems and quantities. Elliptical beams with P θ = 0 in solenoidal or alternating gradient quadrupole focusing channels [35] are beyond the scope of the present analysis. For the special case of a round beam envelope (i.e.,rx =rỹ ≡r b ) with P θ = 0 in a solenoidal focusing channel, the envelope envelope equations can also be mapped to the form of Eq. (B11) with a KV distribution previously analyzed by Chen, Pakter, and Davidson [36]. (In Ref. [36] the authors do not refer to the distribution analyzed as being KV, but it is appropriate to label any distribution as being of KV form if it is a delta-function of a Courant-Snyder invariant that is expressible as a quadratic form in canonically transformed phase-space variables and has a uniform density projection. An infinity of KV distributions can be constructed through linear canonical transformations of the quadratic standard form of the KV distribution [29].) For this special axisymmetric case the self-consistent KV distribution satisfying the Vlasov equation can be expressed as and the corresponding emittances and envelope equation are given by Eqs. (B10) and (B11) withεx =εỹ ≡ε and rx =rỹ =r b . The self-consistent density can be calculated from Eq. (B21) for −1 < 2P θ /ε < 1 (corresponding to sufficiently small angular momentum to be focusable) as and additional moments are summarized in Table II for −1 < 2P θ /ε < 1. Note that the nonzero P θ term results in directed azimuthal flow contributions to various moments. For an axisymmetric beam, it is instructive to employ Eqs. (B3) and (B7) to express the constant Larmor-frame emittanceε in terms of laboratory-frame moments. This givesε where ε x = ε y is the ordinary laboratory-frame emittance given by Eq. (5) and r = x 2 + y 2 and θ = tan −1 (y, x) are the usual cylindrical particle coordinates in the laboratory-frame. Because P θ = const, it follows from Eq. (B23) that can be interpreted as a conserved radial emittance for an axisymmetric beam. In summary, round beam envelopes with r x = r y and P θ = 0 can also be analyzed with the envelope equations (4) provided we interpret the results in the rotating Larmor-frame and substitute in Eq. (4). In this axisymmetric case, the initial phase of orientation of the stationary laboratory-frame is arbitrary. Various authors have derived equivalent axisymmetric beam envelope equations valid for finite P θ including: Lee and Cooper [31], Miller [30], and Chen, Pakter and Davidson [36].
APPENDIX C: DERIVATION OF EQUATION (18) The free-drift differential equation r + = Q/r + [Eq. (17)] can be integrated to give where C 1 = const is the scaled envelope Hamiltonian. We temporarily assume that r + ≥ 0. Then r + = 2Q ln r + + 2C 1 , and integrating once again we have where C 2 = const and erfi(z) = erf(iz)/i = (2/ √ π) z 0 dt exp(t 2 ) is the imaginary error function. Eliminating C 1 and C 2 in terms of the initial conditions r + = r + (s i ) and R + = R + (s i ) at s = s i yields Because r + = Q/r + , r + is an analytic function. Therefore, even though Eq. (C3) was derived under assumption that r + ≥ 0, Eq. (C3) is also valid for r + < 0. Solving Eq. (C3) for r + (s) yields the free-drift solution given in Eq. (18).

APPENDIX D: PHASE ADVANCE OF ENVELOPE MODES IN THE ZERO SPACE-CHARGE LIMIT
In the limit of zero space-charge (Q = 0), the equations of motion [see Eq. (1)] for the particle x-and y-orbits are decoupled and of the same form. Likewise, the resulting x-and y-envelope equations [see Eq. (4)] are decoupled and of the same form. Therefore, we analyze only the x-equation here. The x equation of motion of the particle is given by Here, κ x (s + L p ) = κ x (s) describes the linear applied focusing function of a lattice with periodicity L p that is of otherwise unspecified form. To define the phase advance of the x-orbit, without loss in generality, we express the orbit in phase-amplitude form as [12,13] x(s) Here, A x = const, ψ x (s) = 1/w 2 x (s), and w x (s) is the positive solution to the equation satisfying the periodicity requirement w x (s+L p ) = w x (s). This equation is recognized as the scaled envelope equation with zero space-charge [see Eq. (4)] and therefore the matched beam envelope r xm is related to the amplitude function w x as Consistent with Eqs. (24) and (26), the single-particle phase advance σ 0x of oscillations in x over one lattice period is independent of the initial axial coordinate s i . Consistent with the x-particle equation of motion (D1), the beam envelope equation for r x = 2 x 2 ⊥ is We expand r x as r x = r xm + δr x where r xm = √ x w x is the matched beam radius, and δr x is a small-amplitude perturbation about the matched beam. Then δr x satisfies the linearized envelope equation [see also Eq. (37)] Paralleling the analysis above for the particle orbits, the solution to the linearized envelope equation can be expressed in phase-amplitude form as Here, A e = const, ψ e (s) = 1/w 2 e (s), and w e (s) the positive solution to satisfying the periodicity requirement w e (s + L p ) = w e (s). The phase advance σ ex of oscillations in δr x over one lattice period is given by Comparison of Eqs. (D3) and (D9) shows that the solution for w e is connected to w x as Thus, from Eqs. (D5) and (D10), the phase advances of the single-particle x-orbit and the envelope δr x are related as Equation (D12) proves that in the limit of zero space-charge that the phase advance of oscillations in δr x are two times the phase advance of undepressed particle oscillations in x. An analogous result to Eq. (D12) holds in the yplane. Equation (D12) is straightforward to generalize to show that for zero space-charge and arbitrary, non-periodic lattices with general κ x (s) that the rate of phase accumulation of a small-amplitude envelope perturbation about a non-periodic reference envelope evolution is twice the rate of phase accumulation of the single-particle orbit.
Here, κ(s + L p ) = κ(s) describes the linear applied focusing forces of the lattice with periodicity L p and is otherwise unspecified (κ = k 2 β0 = const for continuous focusing). As in standard treatments of single-particle orbits [12,13], without loss in generality, we express the envelope perturbations in phase-amplitude form as (E2) Here, A ± = const, ψ ± = 1/w 2 ± , and w ± are positive solutions to the equations satisfying the periodicity requirements w ± (s + L p ) = w ± (s). Employing scaling procedures analogous to those presented in Sec. II E, it can be shown that the envelope amplitude functions w ± can be regarded as being parameterized by L p , σ 0 , σ/σ 0 , and parameters that describe the shape of the lattice focusing function κ(s). In contrast to the usual single-particle amplitude functions of undepressed particle orbits [12,13], w ± will depend on matched beam space-charge depression parameter σ/σ 0 in addition to the lattice parameters σ 0 , L p , etc. The phase advance σ ± of oscillations in δr ± over one lattice period is given by independent of the initial coordinate axial s i . Differentiating Eq. (E2) and employing ψ ± = 1/w 2 ± to simplify the result gives Squaring and adding these equations, we obtain quadratic Courant-Snyder invariants in δr ± -δr ± given by These invariants of the envelope perturbations are equations of ellipses in δr ± -δr ± phase-space where the axes of the ellipses do not, in general, align with the δr ± and δr ± axes. They can be interpreted as area measures because as the perturbations evolve the set of all initial conditions satisfying Eq. (E6) (corresponding to all pure modes with the same oscillation amplitude) will trace out an ellipse in δr ± -δr ± phase-space with area πA ± = const. Analysis of the evolution of this constant area ellipse (in terms of elongation and rotation) can be facilitated by introducing Twiss parameters as in standard treatments of uncoupled single-particle dynamics [13]. It is straightforward to show that Courant-Snyder invariants directly analogous to those in Eq. (E6) apply to non-periodic lattices with general κ(s) for any small-amplitude envelope perturbations about a reference envelope evolution that need not be matched. In this non-periodic case, w ± satisfies Eq. (E3) for a general reference envelope solution r m and w ± need not be periodic. Unfortunately, the analysis of Courant-Snyder invariants for the more general case of undriven coupled envelope modes evolving according to [see Eq. (37)] δr j + κ j δr j + 2Q (r xm + r ym ) 2 (δr x + δr y ) + with j = x, y, and κ x = κ y is considerably more complicated than in the continuous or solenoidal focusing case above. This complication arises because when κ x = κ y the normal coordinates of the coupled modes are not simply expressible. In the previous case, the decoupled sum and difference coordinates δr ± = (r x ± r y )/2 provided simple expressions of the normal coordinates. For κ x = κ y , the normal coordinates are expressible in terms of the s-varying eigenvectors of the transfer map in Eq. (47). Fortunately, Edwards and Teng [37] have presented a general formulation for analysis of two-dimensional coupled linear motion that can be directly applied to coupled envelope modes. To make a connection to this work, the equations of motion (E7) of the envelope perturbations can be expressed in where is the envelope Hamiltonian. From this Hamiltonian formulation, the correspondence of variables and coefficients employed by Edwards and Teng to those used here are readily identified. The conjugate variables x, p x and y, p y of Edwards and Teng are replaced by and the coefficients L, F , G, and K are identified as The quadratic invariants calculated by Edwards and Teng can then be directly applied and are generalized Courant-Snyder invariants of the coupled envelope modes. Application and interpretation of these coupled motion invariants are considerably more complicated than in the case presented above for solenoidal and continuous focusing. However, the structure and consequences of the generalized invariants are loosely analogous to the uncoupled motion Courant-Snyder invariants.
is the transfer matrix of the envelope perturbations. This solution is equivalent to the one presented in Eq. (65).
Any linear transform of the envelope coordinates δR can be expressed as where T = const is some invertible 4 × 4 matrix andδ R is the transformed (tilde variables) envelope coordinate vector. The equation of motion in transformed coordinates becomes The solution to Eq. (F5) is given byδ where the transformed transfer mapM is related to the untransformed map M e (s|s i ) bỹ The matrix M e (s i + L p |s i ) has four eigenvalues λ n with n = 1, 2, 3, and 4 defined by M e (s i + L p |s i ) · E n = λ n E n [see Eq. (47)]. Transforming this equation givesM e (s i + L p |s i ) ·Ẽ n = λ nẼn , whereẼ n = T · E n , thereby showing that the eigenvalues of M e are invariant under the transform T. Similarly, the four eigenvalues of the matrix K are equal to the eigenvalues of the transformed matrixK = T · K · T −1 . Consider a transform T to fully decoupled coordinates δR whereM e (s|s i ) is diagonal. Then from Eq. (F8)K is also diagonal. In this representation, the diagonal elements of bothM e (s i + L p |s i ) andK are their respective eigenvalues. The eigenvalues of K can be calculated using Eq. (F1) as ±iσ + /L p and ±iσ − /L p . Therefore, and Eq. (F8) givesM Form Eq. (F10) the eigenvalues of M e (s i + L p |s i ) are identified as λ n = e ±iσ+ , e ±iσ− . To complete the analysis, the fully decoupled solution forδ R(s) given by Eqs. (F7) and Eq. (F10) is connected to the untransformed coordinates δR = (δr x , δr x , δr y , δr y ) byδ R = T · δR through the transformation , It is instructive to examine the connection of this matrix transformation approach to the decoupled sum and difference coordinates δr ± = (δr x ±δr y )/2 employed in Sec. III. Takingδ R = (δr + , δr + , δr − , δr − ) = T·δR corresponds to the linear transformation Under this transformation,K is block diagonal with The corresponding transfer mapM e (s|s i ) = exp[−K(s − s i )] can be directly calculated using the series definitions of the matrix exponential, recognizing the recursion among terms, and summing resulting series to givẽ The block diagonal form of Eq. (F13) corresponds to rotations in the decoupled r + −r + and r − −r − phase-spaces. This solution also follows immediately from Eq. (65). Finally, it is instructive to point out that the transformation to fully decoupled coordinates presented in Eq. (F11) can be derived in terms of two successive linear transformations, first a transformation to decoupled δr ± coordinates as given by Eq. (F12)  Normal envelope modes of thin-lens (η → 0 with σ 0 fixed) solenoidal and FODO (α = 1/2) quadrupole focusing channels can be analyzed analytically in the limit of zero (σ → σ 0 ) and full (σ → 0) space-charge depression [38]. In the zero space-charge limit, the x-and y-envelope coordinates decouple and the general analysis in II F applies to show that all envelope modes have phase advance 2σ 0 . As a consequence of this phase advance limit and the decoupling, instabilities are possible only for σ 0 = 90 • and 180 • . A straightforward transfer matrix analysis of single-particle orbits can then be used to show that there is no envelope instability for σ 0 = 90 • , and there is linear growth in oscillation amplitude with lattice periods traversed for σ 0 = 180 • .
To analyze the opposite limit of full space-charge depression, we proceed as follows. Let L be the length of the free-drift interval between the two subsequent thin-lenses, solenoids or quadrupoles (L p = L for solenoids and L p = 2L for quadrupoles) as illustrated in Fig. 24. By symmetry we only need to consider the envelope evolution of the beam between two lenses only. We take the first lens to be at axial location s = −L/2 and the second one to be at s = L/2.
The envelope equations (4) for the case of a fully depressed beam (with σ = 0), which corresponds to zero beam emittance (ε x = ε y = 0), can be written in terms of scaled sum and difference coordinates R ± = (r x ± r y )/(2 √ 2Q) as: for quadrupole focusing. Here, we have employed κ x = κ y for solenoidal focusing and κ x = −κ y for quadrupole focusing. For the solenoidal channel all thin-lens focusing kicks are focusing in x. Without loss of generality, we assume that in the alternating gradient channel that the lens at s = L/2 is focusing in x. Then for both thin-lens solenoids and quadrupoles we take near s = L/2 [see Eq. (13)] where f = const is the thin-lens focal length and δ(s) is the Dirac delta-function. The focal length f can be related to the undepressed particle phase advance over one lattice period σ 0 as [see Eq. (30)] L f = 2 − 2 cos σ 0 , solenoidal focusing, √ 2 − 2 cos σ 0 , quadrupole focusing.
To analyze the envelope stability of the fully space-charge depressed beam in solenoidal and quadrupole focusing channels we examine the change in the envelope coordinate vector R(s) = (R + (s), R + (s), ζ(s)R − (s), ζ(s)R − (s)) from the mid-drift at s = 0 to the next mid-drift at s = L. Here, ζ(s) = 1 when the next lens to be traversed is focusing, and ζ(s) = −1 when the next lens is defocusing. For solenoidal focusing ζ(s) = 1, and for alternating gradient quadrupole focusing ζ(s) = −1 for 0 ≤ s ≤ L/2 and ζ(s) = 1 for L/2 < s ≤ L. The effect of ζ(s) is to map the matched beam envelope R − in the second drift region of the quadrupole (L/2 < s ≤ L) into the first drift region (0 ≤ s < L/2) [see Fig. 24(c)]. This mapping expedites the use of symmetry in the matched envelope to simplify calculations for the FODO lattice.
To analyze the first-order perturbations in the coordinate vector R(s) we compute the Jacobian matrix M(0, L) where M(s 1 |s 2 ) = ∂R(s 2 )/∂R(s 1 ) and derivatives are evaluated for a matched beam envelope. The Jacobian matrix M(0|L) is identical to the symplectic transfer map in Eq. (46). Analogously to the treatment in Sec. II F, the stability of the envelope to linear perturbations about the matched envelope solution corresponds to M(0|L) having all eigenvalues on the complex unit circle. To characterize mode properties, eigenvalues will be calculated in terms of σ 0 .
To evaluate M s , we consider the action of the thin-lens according to Eq. (14). We obtain for the solenoidal channel and for the quadrupole channel keeping in mind that the ζ(s)R − changes sign from one free-drift region to the next, To complete the evaluation of M f (L/2 − 0), we find relations of the elements to σ 0 and L by deriving equations connecting R + (L/2 − 0) ≡ R + (L/2), R + (L/2 − 0), and R + (0) to these quantities for the matched beam envelope. By symmetry, for a periodic matched envelope Combining these constraints with the matching conditions (G8), we get Similarly, using Eqs. (G1b) and (G2) for alternating gradient focusing and matched beam symmetries (G8), we obtain For quadrupole focusing, Eqs. (G9b) and (14) can be combined to yield an additional constraint Both the solenoidal and quadrupole matching conditions in Eqs. (G9) for R + can be expressed as wherek = L 2f = 1 − cos σ 0 , solenoidal focusing, where M ± (0|L) are 2 × 2 symplectic matrices that can be independently analyzed for the stability of perturbations in R + (breathing mode) and R − (quadrupole mode). Mode phase advances (σ ± ) and growth factors (γ ± ) can be calculated from the eigenvalues λ ± of M ± (0|L). Stable modes with λ ± on the complex unit circle correspond to (1/2) |Tr M ± (0|L)| < 1 [see Eq.