Efficient computation of matched solutions of the Kapchinskij-Vladimirskij envelope equations for periodic focusing lattices

A new iterative method is developed to numerically calculate the periodic, matched beam envelope solution of the coupled Kapchinskij-Vladimirskij (KV) equations describing the transverse evolution of a beam in a periodic, linear focusing lattice of arbitrary complexity. Implementation of the method is straightforward. It is highly convergent and can be applied to all usual parameterizations of the matched envelope solutions. The method is applicable to all classes of linear focusing lattices without skew couplings, and also applies to all physically achievable system parameters -- including where the matched beam envelope is strongly unstable. Example applications are presented for periodic solenoidal and quadrupole focusing lattices. Convergence properties are summarized over a wide range of system parameters.


INTRODUCTION
The Kapchinskij-Vladimirskij (KV) envelope equations [1,2,3] are often employed as a simple description of the transverse evolution of intense ion beams. The equations are coupled ordinary differential equations that describe the evolution of the beam edge (or rms radii) in response to applied linear focusing forces of the lattice and defocusing forces resulting from beam space-charge and transverse phase-space area (emittances). Although the KV envelope equations are only fully Vlasov consistent with the singular KV distribution, the equations can be applied to describe the low-order evolution of a real distribution of beam particles when the variation of the statistical beam emittances is negligible or sufficiently slow [2]. Large nonlinear fields that can be produced by non-ideal applied focusing elements, nonuniform beam space-charge, and species contamination (electron cloud effects, etc.) drive deviations from the KV model. Such effects are suppressed to the extent possible in most practical designs, rendering the KV model widely applicable.
The matched solution of the KV envelope equations is the solution with the same periodicity as the focusing lattice [1,2,3]. The matched beam envelope is important because it is believed to be the most radially compact solution supported by a periodic linear focusing channel [3]. Matched envelopes are typically calculated as a first step in the design of practical transport lattices and for use in initializing more detailed beam simulations to evaluate machine performance [2]. The matched envelope solution is typically calculated by numerically integrating trial solutions of the KV equations from assumed initial conditions over one lattice period and searching for the four initial envelope coordinates and angles that generate the solution with the periodicity of the lattice [2,3,4]. An optimal formulation of the conventional root finding procedure for envelope matching has been presented by Ryne [4]. Conventional root finding procedures for matching can be surprisingly problematic even for relatively simple focusing lattices. Variations in initial conditions can lead to many inflection points in the envelope functions at the end of the lattice period. Thus initial guesses close to the actual values corresponding to the periodic solution are often necessary to employ standard root finding techniques. This is especially true for complicated focusing lattices with low degrees of symmetry and where focusing strengths (or equivalently, undepressed single particle phase advances) are large. For large focusing strength and strong space-charge intensity, the matched envelope solution can be unstable over a wide range of system parameters [2,3]. Such instabilities can restrict the basin of attraction when standard numerical root finding methods are used to calculate the needed matching conditions -especially for certain classes of solution parameterizations.
In this article we present a new iterative procedure to numerically calculate matched envelope solutions of the KV equations. The basis of this procedure is that the KV distribution of particle orbits internal to the beam must have a locus of turning points consistent with the beam edge (envelope). In the absence of beam space-charge, betatron amplitudes calculated from the sine-and cosine-like principal orbits describing particles moving in the applied focusing fields of the lattice directly specify the matched beam envelope [5]. For finite beam space-charge, the principal orbits describing the betatron amplitudes and matched beam envelope cannot be a priori calculated because the defocusing forces from beam space-charge uniformly distributed within the (undetermined) beam envelope are unknown. In the iterative matching (IM) method, the relation between the betatron amplitudes and the particle orbits is viewed as a consistency equation. Starting from a simple trial envelope solution that accounts for both space-charge and and applied focusing forces in a general manner, the consistency condition is used to iteratively correct the envelope functions until converged matched envelope solutions are obtained that are consistent with particle orbits internal to the beam.
The IM method offers superior performance and reliability in constructing matched envelopes over conventional root finding because the IM iterations are structured to reflect the periodicity of the actual matched solution rather than searching for parameters that lead to periodicity. The IM method works for all physically achievable system parameters (even in cases of envelope instability) and is most naturally expressed and rapidly convergent when relative beam space-charge strength is expressed in terms of the depressed particle phase advance. All other parameterizations of solutions (specified perveances and emittances, etc.) can also be carried out by simple extensions of the IM method rendering the approach fully general. The natural depressed phase advance parameterization is also useful when carrying out parametric studies because phase advances are the most relevant parameters for analysis of resonancelike effects central to charged particle dynamics in accelerators. The IM method provides a complement to recent analytical perturbation theories developed to construct matched beam envelopes in lattices with certain classes of symmetries [6,7,8,9]. In contrast to these analytical theories, the IM method can be applied to arbitrary linear focusing lattices without skew couplings. The highly convergent iterative corrections of the IM method have the same form for all order iterations after seeding, rendering the method straightforward to code and apply to numerically generate accurate matched envelope solutions.
The organization of this paper is the following. After a review of the KV envelope equations in Sec. II, various properties of matched envelope solutions and the continuous focusing limit are analyzed in Sec. III. These results are used in Sec. IV to formulate the IM method for calculation of matched solutions to the KV envelope equations. Example applications of the IM method are presented in Sec. V to illustrate application and convergence properties of the method over a wide range of system parameters for a variety of systems. Concluding comments in Sec VI summarize the advantages of the IM method over conventional techniques.

II. THEORETICAL MODEL
We consider an unbunched beam of particles of charge q and mass m coasting with axial relativistic factors β b = const and γ b = 1/ 1 − β 2 b . In the KV model, the beam is propagating in a linear focusing lattice without skew couplings and has uniform charge density within an elliptical cross-section with principal radii r x and r y along the (transverse) x-and y-coordinate axes. When self-fields are included and image effects are neglected, the envelope radii consistent with the KV distribution evolve according to the so-called KV envelope equations [1,2,3] Here, primes denote derivatives with respect to the axial machine coordinate s, the subscript j ranges over both transverse coordinates x and y, the functions κ j (s) represent linear applied focusing forces of the transport lattice, Q = const is the dimensionless perveance, and ε j = const are the rms edge emittances. Equations relating the functions κ j to magnetic and/or electric fields of practical focusing elements are presented in Ref. [3]. The perveance provides a dimensionless measure of self-field defocusing forces internal to the beam [2] and is defined as Here, I is the constant beam current, c is the speed of light in vacuo, and ǫ 0 is the permittivity of free space. The perveance Q can be thought of as a scaled measure of space-charge strength [2]. The rms edge emittances ε j provide a statistical measure of beam phase-space area projections in the x-x ′ and y-y ′ planes [2]. When the emittances are constant (ε j = const), the KV envelope equations (1) are consistent with the Vlasov equation only for the KV distribution [1,10], which is a singular function of Courant-Snyder invariants. This singular structure can lead to unphysical instabilities within the Vlasov model [11]. However, the KV envelope equations can be applied to physical (smooth) distributions in an rms equivalent beam sense [2], with the envelope radii and the emittances defined by statistical averages of the physical distribution as and Here, · · · denotes a transverse statistical average over the beam distribution function, and for notational simplicity, we have assumed zero centroid offset (e.g., x = 0). In this rms equivalent sense, the emittances ε j will generally evolve in s. If this variation has negligible effect on the r j , then the KV envelope equations can be applied with ε j = const to reliably model practical machines. This must generally be verified a posteriori with simulations of the full distribution. For appropriate choices of the lattice focusing functions κ j (s), Eq. (1) can be employed to model a wide range of transport channels, including solenoidal and quadrupole transport. For solenoidal transport, the equations must be interpreted in a rotating Larmor frame (see Appendix A of Ref. [3]). In a periodic transport lattice, the κ j are periodic with fundamental lattice period L p , i.e., The beam envelope is said to be matched to the transport lattice when the envelope functions have the same periodicity as the lattice: For specified focusing functions κ j (s), beam perveance Q, and emittances ε j , the matching condition is equivalent to requiring that r j and r ′ j satisfy specific initial conditions at s = s i when the envelope equations (1) are integrated as an initial value problem. The required initial conditions generally vary with the phase of s i in the lattice period (because the conditions vary with the local matched solution). In conventional procedures for envelope matching, needed initial conditions are typically found by numerical root finding starting from guessed seed values [3]. This numerical matching can be especially problematic when: applied focusing strengths are large, the focusing lattice is complicated and devoid of symmetries that can reduce the dimensionality of the root finding, choices of solution parameters require extra constraints to effect, and where the matched beam envelope is unstable.
The undepressed particle phase advance per lattice period σ 0j provides a dimensionless measure of the strength of the applied focusing functions κ j describing the periodic lattice [3,5]. The σ 0j can be calculated from [5] cos σ 0j = 1 2 Tr M 0j (s i + L p |s i ), where M 0j (s|s i ) denotes the 2×2 single particle transfer matrix in the j-plane from axial coordinate s i to s. Explicitly, we have where the C 0j (s|s i ) and S 0j (s|s i ) denote cosine-like and sine-like principal orbit functions satisfying with F representing C or S with C 0j subject to cosine-like initial (s = s i ) conditions C 0j (s i |s i ) = 1 and C ′ 0j (s i |s i ) = 0, and with S 0j subject to sine-like initial conditions S 0j (s i |s i ) = 0 and S ′ 0j (s i |s i ) = 1. Equation (7) can be expressed in terms of C 0j and S ′ 0j as The σ 0j are independent of the particular value of s i used in the calculation of the principal functions. For some particular cases such as piecewise constant κ j the principal functions F 0j can be calculated analytically. But, in general, the F 0j must be calculated numerically. In the absence of space-charge, the particle orbit is stable whenever σ 0j < 180 • and parametric bands of stability can also usually be found for σ 0j > 180 • [2,3,12]. For a stable orbit, the scale of the κ j (i.e., κ j → ακ j with α = const setting the scale of the specified κ j ) can always be regarded as being set by the σ 0j . In this context, Eq. (10) is employed to fix the scale of the κ j in terms of σ 0j and other parameters defining the κ j . Because there appears to be no advantage in using stronger focusing with σ 0j > 180 • in terms of producing more radially compact matched envelopes [3,13], we will assume in all analysis that follows that the κ j are sufficiently weak to satisfy σ 0 < 180 • . The formulation given above for calculation of the undepressed principal orbits C 0j and S 0j and the undepressed particle phase advances σ 0j can also be applied to calculate the depressed principal orbits C j and S j and the depressed phase advances σ j in the presence of uniform beam space-charge density for a particle moving within the matched KV beam envelopes. This is done by replacing in Eqs. (9) and dropping the subscript 0s in Eqs. (7)-(10) for notational clarity (i.e., C 0j → C j and S 0j → S j ). Explicitly, the depressed principal functions satisfy with F representing C or S with C j subject to C j (s i |s i ) = 1 and C ′ j (s i |s i ) = 0, and S j subject to S j (s i |s i ) = 0 and S ′ j (s i |s i ) = 1, and the depressed phase advances satisfy For a stable orbit, it can be shown that the σ j can also be calculated from the matched envelope as [3,5] This formula can also be applied to calculate σ 0j by using the matched envelope functions r j calculated with Q = 0. Matched envelope solutions of Eqs. (1) can be regarded as being determined by the focusing functions κ j , the perveance Q, and the emittances ε j . The lattice period L p is implicitly specified through the κ j . We will always regard the scale of the κ j as being set by the undepressed phase advances σ 0j through Eq. (10). For σ 0j < 180 • there is no ambiguity in scale choice and the use of the σ 0j as parameters allows disparate classes of lattices to be analyzed in a common framework [3]. The depressed phase advances σ x and σ y can be employed to replace up to two of the three parameters Q, ε x , and ε y . Such replacements can be convenient, particularly when carrying out parametric surveys (for example, see Ref. [3]) because σ j /σ 0j is a dimensionless measure of space-charge strength satisfying 0 ≤ σ j /σ 0j ≤ 1 with σ j /σ 0j → 1 representing a warm beam with negligible space-charge (i.e., Q → 0, or ε j → ∞ for finite Q), and σ j /σ 0j → 0 representing a cold beam with maximum space-charge intensity (i.e., ε j → 0). We will discuss calculation of matched beam envelopes for the useful parameterization cases listed in Table I. In cases typical of linear accelerators the focusing functions have equal strength in the x-and y-planes giving σ 0x = σ 0y . In such plane symmetric cases we denote σ 0j ≡ σ 0 . In practical situations where the focusing lattice and emittances are both plane symmetric with σ 0j ≡ σ 0 and ε j ≡ ε, then the depressed phase advance is also plane symmetric with σ j ≡ σ and parameterization cases 2 and 3 are identical. It is assumed that a unique matched envelope solution exists independent of the parameterization when the κ j are fully specified. There is no known proof of this conjecture, but numerical evidence suggests that it is correct for simple focusing lattices (i.e., simple κ j ) when σ 0j < 180 • . In typical experimental situations, note that transport lattices are fixed in geometry and excitations of focusing elements in the lattices can be individually adjusted. In the language adopted here, such lattices with different excitations in focusing elements (both overall scale and otherwise) correspond to different lattices described by different κ j with corresponding different matched envelopes.

III. MATCHED ENVELOPE PROPERTIES
In development of the IM method in Sec. IV, we employ a consistency equation between depressed particle orbits within the beam and the matched envelope functions (III A) and use a continuous focusing description of the matched beam (III B) to construct a seed iteration. Henceforth, we denote lattice period averages with overbars, i.e., for some quantity ζ(s), κj (σ0j ), εj, and one of σj 3 κj (σ0j ), σj, and one of εj A. Consistency condition between particle orbits and the matched envelope We calculate nonlinear consistency conditions for the matched envelope functions r j and the depressed principal orbit functions C j and S j as follows. First, the transfer matrix M j of the depressed particle orbit in the j-plane is expressed in terms of betatron function-like formulation as [5] with, Here, is the change in betatron phase of the particle orbit from s = s i to s and the principal functions C j and S j are calculated including the linear space-charge term of the uniform density elliptical beam from Eq. (12) . Note that r j ≡ ε j β j can be used in Eqs. (17) and (18) to express the results more conventionally in terms of the betatron amplitude functions β j describing linear orbits internal to the beam in the j-plane [5]. These generalized betatron functions are periodic [i.e., β j (s + L p ) = β j (s) and include the transverse defocusing effects of uniformly distributed space-charge within the KV equilibrium envelope. Recognizing that ∆ψ j (s i + L p ) = σ j [see Eq. (14)] and that the matched envelope functions r j have period L p gives Here, [M j ] 12 denotes the 1, 2 component of the 2 × 2 matrix M j and σ j can be equivalently calculated from either Eq. (13) or Eq. (14). Equation (19) can be applied to numerically calculate the consistency conditions for the matched envelope functions r j on a discretized axial grid of s locations. As written, the principal orbit functions employed (i.e., the C j and S j ) need to be independently calculated at each s-location on the grid through one lattice period. The fact that every period is the same can be applied to simplify the calculation. For any initial axial coordinate s i we have M j (s + L p |s) = M j (s + L p |s i + L p ) · M j (s i + L p |s). Multiplying this equation from the right side by the identity matrix Some straightforward algebra employing Eqs. (16), (19), and (20), and the Wronskian (or symplectic) condition on Equation (22) explicitly shows that the linear principal functions C j and S j need only be calculated in s from some arbitrary initial point (s i ) over one lattice period (to s i + L p ) to calculate the consistency condition for the matched envelope functions r j (s), or equivalently, the betatron amplitude functions β j (s) ≡ r 2 j (s)/ε j . Equation (22) can also be derived using Courant-Snyder invariants of particle orbits within the beam.
Equations (13) and (22) form the foundation of an iterative numerical method developed in Sec. IV to calculate the matched beam envelope for any lattice. These equations express the intricate connection between the bundle of depressed particle orbits within the uniform density KV beam and the locus of maximum particle excursions defining the envelope functions r j . The method will be iterative because the consistent matched envelope functions r j are necessary to integrate the linear differential equations for the depressed orbit principal functions C j and S j . However, in the limit Q → 0, the principal functions do not depend on the r j and the matched envelope can be immediately calculated from the equations. Thus, the periodic zero-current matched beam envelope can be directly calculated using Eq. (22) in terms of the two independent, aperiodic linear orbits (i.e., C 0j and S 0j ) integrated over one lattice period.
Additional constraints on the matched envelope functions r j and/or betatron functions β j are necessary to formulate the IM method for parameterizations where one or more of the parameters Q and ε j need to be eliminated (see Table I). Appropriate constraints can be derived by taking the period average of Eq. (1) for a matched envelope, giving

B. Continuous limit
In the continuous focusing approximation, we take the lattice focusing functions κ j as constants set according to with the σ 0j calculated from Eq. (10) consistent with the actual s-varying periodic focusing functions κ j . Then we replace r j → r j in the KV envelope equations (1) and take r j = const to obtain the continuous limit envelope equation Equation (25) provides an estimate of the lattice period average envelope radii r j in response to applied focusing forces and defocusing forces from beam space-charge and thermal (emittance) effects. Solutions for r j will be employed to seed the IM method of constructing matched envelope solutions. In general, the continuous limit approximations tend to be more accurate for weaker applied focusing strengths with σ 0j < ∼ 80 • . However, even for higher values of σ 0j < 180 • , the formulas can still be applied to seed iterative numerical matching methods if the methods have a sufficiently large "basin of attraction" to the desired solution.
For case 0 parameterizations (specified σ 0j , Q, and ε j ) the solutions of Eq. (25) will, in general, need to be calculated numerically from a trial guess. Certain limits are analytically accessible and often relevant. If the beam perveance Q is zero, or equivalently if σ j = σ 0j , then Eqs. (25) decouple and are trivially solved as Alternatively, this result can be obtained using r j = r j in Eq. (14) with σ j → σ 0j . In the case of a symmetric system with σ 0x = σ 0y ≡ σ 0 and ε x = ε y ≡ ε, then r x = r y ≡ r b and the envelope equations decouple and the resulting quadratic equation in r 2 b is solved as In parameterization cases 1-3, the continuous limit solutions r j must be expressed using the depressed phase advances σ j to eliminate one or more of the parameters Q and ε j . In these cases, if the emittances ε j are known, then Eq. (14) can be employed to estimate Alternatively, if the perveance Q is known but one or more of the emittances ε j is unknown, we can use Eq. (28) to eliminate the emittance term(s) in Eq. (25) obtaining (σ 2 0j − σ 2 j )r j = 2QL 2 p /(r x + r y ). Taking the ratio of the x-and y-equations yields Back-substitution of this result in (σ 2 0j − σ 2 j )r j = 2QL 2 p /(r x + r y ) then gives . (30) Alternative smooth-limit formulations in the literature [14,15] can also be employed to estimate the r j for systems with high degrees of symmetry.

IV. NUMERICAL ITERATIVE METHOD FOR MATCHED ENVELOPE CALCULATION
We formulate a numerical iterative matching (IM) method to construct the matched beam envelope functions r j (s) over one lattice period L p using the developments in Sec. III. The IM method is formulated for arbitrary periodic focusing functions κ j . Constraints necessary to apply the IM formalism to all cases of envelope parameterizations listed in Table I are derived. Label all quantities varying with iteration number with a superscript i (i = 0, 1, 2, · · · ) denoting the iteration order. For example, the ith order envelope functions are labeled r i j . The iteration label should not be confused with the initial coordinate s i and the initial "seed" iteration corresponds to i = 0. Parameters such as the perveance Q or emittances ε j will also be superscripted to cover parameterization cases where the quantities are unspecified and are calculated from the envelope functions and other parameters (see Table I). For example, ε i j denotes the j-plane emittance at the ith iteration and for parameterization cases where the value of ε j is specified, then ε i j = ε j = const. For iterations i ≥ 1, we calculate refinements of the principal orbit functions [see Eqs. (9) and (11)] in terms of the envelope calculated at the previous, i − 1 iteration from Here, F i j denotes C i j (s|s i ) or S i j (s|s i ) which are subject to the initial (s = s i ) conditions C i j (s i |s i ) = 1, C i ′ j (s i |s i ) = 0 and S i j (s i |s i ) = 0, S i ′ j (s i |s i ) = 1. Note that the F i j depend on the envelope functions and perveance of the prior, Here, if the parameterization does not specify the depressed phase advances as σ i j = σ j , then they are calculated [see Eq. (13)] for all i from In parameterization cases 0 to 3 (see Table I), one or more of the needed quantities among Q i , ε i j , and σ i j are not specified (e.g., for case 1: ε i j = ε j specified) and must be calculated to apply Eq. (32) and/or to calculate the next (i + 1) iteration principal functions from Eq. (31). Equations (33) and/or the constraint equations (23) with Q → Q i , ε j → ε i j , and r j → r i j (or in some cases r j → ε i j β i j ) can be employed to calculate parameter eliminations necessary to fully realize each iteration as follows for each case: Case 0 (κ j , Q, ε j specified) The σ i j can be calculated from Eq. (33).
Case 1 (κ j , Q, and σ j specified) The ε i j can be calculated using Eq. (23) expressed in betatron form to obtain with the ratio ε i y /ε i x on the right hand side of the equations determined by Note that expressing the constraints in terms of betatron functions β i j is necessary in this case because the envelope functions r i j cannot be calculated from Eq. (32) until the ε i j are known, whereas because of the structure of the envelope equations, the β i j = (r i j ) 2 /ε i j can be calculated from Eq. (32) without a priori knowledge of the values of the ε i j .
Case 2 (κ j , ε j , and σ x specified; or κ j , ε j , and σ y specified) If necessary, either σ i x or σ i y can be calculated from Eq. (33) to enable full specification of the functions β i j or r i j . Then, Q i can be calculated using Eq. (34) and the β i j , or alternatively, using with ε i j = ε j .
Case 3 (κ j , σ j , and ε x specified; or κ j , σ j , and ε y specified) First, Eq. (35) and the β i j functions can be applied to calculate ε i y from specified ε x , or ε i x from specified ε y . Then, Q i can be calculated from the ε i j (if specified, ε i j = ε j ) using Eq. (34) and the β i j , or alternatively, with Eq. (36) and the r i j .
The seed i = 0 iteration is treated as a special case where the continuous limit formulas derived in Sec. III B are applied to estimate the leading order defocusing effect of space-charge on the beam. In this case the principal functions are calculated from Here, F 0 j denotes C 0 j (s|s i ) or S 0 j (s|s i ) subject to the initial (s = s i ) conditions C 0 j (s i |s i ) = 1, C 0 ′ j (s i |s i ) = 0 and S 0 j (s i |s i ) = 0, S 0 ′ j (s i |s i ) = 1, and Q and r j denote the continuous focusing approximation perveance and envelopes calculated from the formulation in Sec. III B with Q → Q and ε j → ε j . The continuous focusing values of Q and ε j used in calculating the r j are set by the parameterization values in cases where they are specified (e.g., Q = Q for Q specified). Otherwise, Q and/or the ε j are calculated in terms of other parameters using the appropriate constraint equations from Eqs. (26)-(30) applied with Q → Q and ε j → ε j .
Note that the seed envelope functions r 0 j calculated under this procedure are not the continuous limit functions (i.e., r 0 j = r j ). Likewise, in parameterizations where they are not held fixed, the seed perveance and emittances will not equal the continuous focusing values (i.e., Q 0 = Q and/or ε 0 j = ε j ). Due to Eq. (32), the seed envelope functions r 0 j will have a (dominant) contribution to the envelope flutter from the applied focusing fields of the lattice with a correction due to space-charge defocusing forces from the continuous limit formulas. This approximation should produce seed envelope functions r 0 j that are significantly closer to the actual periodic envelope functions r j than would be obtained by simply applying continuous limit formulas (i.e., taking r 0 j = r j ) or neglecting the effects of space-charge altogether [i.e., by calculating r 0 j using Eq. (22) with Q = 0]. Generally speaking, a seed iteration closer to the desired solution can reduce the number of iterations required to achieve tolerance, and more importantly, can help ensure a starting point within the basin of attraction of the method, thereby reducing the likelihood of algorithm failure. At the expense of greater complexity and less lattice generality, alternative seed iterations can be generated using low order terms from analytical perturbation theories for matched envelope solutions [6,7,8,9]. In certain cases, these formulations may generate seed iterations closer to the matched solution.
Iterations can be terminated at some value of i where the maximum fractional change between the i and (i − 1) iterations is less than a specified tolerance tol, i.e., where Max denotes the maximum taken over the lattice period L p and the component index j = x, y. Many numerical methods will be adequate for solving the linear ordinary differential equations for the principal functions C i j and S i j of the iteration because they are only required over one lattice period. Generally, the principal functions will be solved at (uniformly spaced) discrete points in s over the lattice period. These discretized solutions can then be employed with quadrature formulas to calculate any needed integral constraints to affect the envelope parameterizations given in Table I. Finally, the convergence criterion (38) can be evaluated at the discrete s-values of the numerical solution for r i j at the ith iteration using saved i − 1 iteration values for r i−1 j that are needed for calculation of the ith iteration. It is usually sufficient to evaluate the convergence criterion at some limited, randomly distributed sample of s-values within the lattice period. Issues of convergence rate and the basin of attraction of the method are parametrically analyzed for examples corresponding to typical classes of transport lattices in Sec. V.

V. EXAMPLE APPLICATIONS
In this section we present examples of the IM method developed in Sec. IV to construct matched envelope solutions and explore parametric convergence properties for example solenoidal and quadrupole periodic focusing lattices. For simplicity, examples are restricted to plane-symmetric focusing lattices with equal undepressed particle phase advances in the x-and y-planes (i.e., σ 0j ≡ σ 0 ) and a symmetric beam with equal emittances in both planes (ε j ≡ ε). Under these assumptions, the depressed phase advances σ j are also equal in both planes (σ j ≡ σ) and parametrization cases 2 and 3 of Table I are identical. First, parameterization cases 1 (specified κ j , Q and σ) and 2 (specified κ j , ε, and σ) are examined in Sec. V A. Both of these cases have specified depressed phase advance σ and represent the most "natural" parameterization of the IM method. Then results in Sec. V A are extended in Sec. V B to illustrate how the IM method can be applied to parameterization case 0 (specified κ j , Q, and ε) with unspecified σ while circumventing practical implementation difficulties.
For simplicity, we further restrict our examples to periodic solenoidal and quadrupole doublet focusing lattices with piecewise constant focusing functions κ j (s) as illustrated in Fig. 1. Solenoidal focusing has κ x (s) = κ y (s), and alternating gradient quadrupole focusing has κ x (s) = −κ y (s). For both the solenoid and quadrupole lattices illustrated, η ∈ (0, 1) is the fractional occupancy of the focusing elements in the lattice period L p . The focusing strength of the elements is taken to be |κ j | =κ = const within the axial extent of the optics and zero outside. For solenoids, κ j =κ > 0 in the focusing element; and for quadrupoles κ x = −κ y =κ > 0 in the focusing-in-x element of the doublet, and κ x = −κ y = −κ < 0 in the defocusing-in-x element. The free drift between solenoids has axial length d = (1 − η)L p . For quadrupole doublet focusing, the two drift distances d 1 = α(1 − η)L p and d 2 = (1 − α)(1 − η)L p separating focusing and defocusing quadrupoles can be unequal (i.e., d 1 = d 2 ). A syncopation parameter α ∈ [0, 1] provides a measure of this asymmetry. Without loss of generality, the lattice can always be relabeled to take α ∈ [0, 1/2], with α = 0 corresponding to the focusing and defocusing lenses touching each other [d 1 = 0 and d 2 = (1 − η)L p ] and α = 1/2 corresponding to a so-called FODO lattice with equally spaced drifts [d 1 = d 2 = (1 − η)L p /2]. These focusing lattices are discussed in more detail in Ref. [3], including a description of how the focusing strength parameterκ is related to magnetic and/or electric fields of physical realizations of the focusing elements.  As mentioned, for general lattices the scale of the focusing functions κ j can be set by the undepressed phase advances σ 0j using Eq. (10). For the piecewise constant κ j defined in Fig. 1, this calculation [3] shows that the focusing strength |κ| is related to σ 0 by the constraint equations: Here, for both solenoidal and quadrupole focusing lattices, Θ = |κ|ηL p /2. In the analysis that follows, Eq. (39) is employed to numerically calculate Θ for a specified value of σ 0 , and thenκ is calculated in terms of other specified lattice parameters as |κ| = 4Θ 2 /(ηL p ) 2 . The undepressed phase advance σ 0 is measured in degrees per lattice period.
Integrations of needed principal orbits to implement the IM method are carried out with the with initial conditions (s = s i ) corresponding to the axial middle of drifts separating focusing elements. Typical matched envelope solutions r j are shown for one lattice period in Fig. 2 for solenoid, FODO (α = 1/2) quadrupole, and syncopated (α = 1/2) quadrupole focusing lattices. Scaled x-plane lattice focusing functions κ x are shown superimposed. Excursions of the matched envelope functions are in-phase for solenoidal focusing (r x = r y ) because the applied focusing is plane symmetric (κ x = κ y ). In contrast, for quadrupole focusing, the anti-symmetric plane focusing (κ x = −κ y ) results in out of phase envelope flutter in each plane (focus-defocus) leading to net focusing over the lattice period in both planes. Expected symmetries of the matched solutions are present for both the solenoidal and quadrupole focusing lattices (see Appendix A). For the quadrupole solutions, note that the FODO case exhibits a higher degree of sub-period symmetry than the syncopated case. Leading order terms of an analytical perturbation theory for the matched beam envelope solution [7] can be applied in the limit σ → 0 to show that the envelope excursions (flutter) scales as Solenoidal Focusing.
Equation (40) shows that for solenoidal focusing the matched envelope flutter increases with decreasing lattice occupancy η and increasing focusing strength σ 0 . In contrast, for FODO quadrupole focusing the flutter depends weakly on η (the variation of Max[r x ]/r x − 1 in η has a maximum range of 0.07) and more strongly on σ 0 (variation of 0.5). Envelope flutter changes only weakly when space-charge strength is reduced (i.e., σ/σ 0 increased). Although the system symmetries assumed simplify interpretation of the matched envelope solutions obtained in the examples, we note that the numerical methods employed in calculation of the particle principal orbits functions and any necessary constraint equations are not structured to take advantage of the symmetries of the matched solutions. Because of this, the examples provide a better guide to the performance of the IM method in situations where there are lesser degrees of system symmetry. Mathematica [16] based programs used to in the examples have been archived [17]. These programs can be easily adapted to more complicated lattices.
A. Case 1 and 2 parameterizations The IM method described in Sec. IV is applied with σ i = σ specified and the unknown parameters of the ith iteration ε i (case 1: Q and σ specified) or Q i (case 2: ε and σ specified) calculated from the constraint equations (34) and (35). The continuous focusing approximation envelope radii r j used in the seed (i = 0) iterations are calculated from Eq. (30) (case 1) and Eq. (28) (case 2). The number of iterations needed to achieve a 10 −6 fractional envelope tolerance [see Eq. (38)] are presented in Fig. 3 as a function of σ 0 and σ/σ 0 for solenoidal, FODO quadrupole, and syncopated quadrupole focusing lattices employing both case 1 and case 2 parameterization methods. In Fig. 4, iterations corresponding to one data point in Fig. 3 (case 2 solenoidal focusing)   The data in Fig. 3 shows that the IM method converges rapidly to small tolerances over a broad range of applied focusing (σ 0 ) and space-charge (σ/σ 0 ) strength. Not surprisingly, stronger focusing strength (i.e., increasing σ 0 ) requires more iterations for both solenoidal and quadrupole focusing at the fixed value of lattice occupancy η employed. Also, lesser degrees of lattice symmetry result in more iterations being necessary for convergence (e.g., lattice convergence rate order: solenoidal, FODO quadrupole, syncopated quadrupole). Iterations required appear to depend only weakly on space-charge strength (σ/σ 0 ) -except for solenoidal focusing lattices with very high σ 0 where required iterations become abruptly larger for weak space-charge with σ/σ 0 close to unity. Even parameters deep within the regime of strong envelope instability [3] converge rapidly. Points for σ/σ 0 = 1 are eliminated in the case 1 examples because the perveance Q is held to a fixed, finite value and this limit would correspond to a matched beam envelope with infinite cross-sectional area. Conversely, for the limit σ/σ 0 = 1 in the case 2 examples, only one iteration is required for convergence to a finite solution because for zero space-charge strength the trial seed iteration generated by Eq. (32) for i = 0 corresponds to the exact matched envelope to numerical error (i.e., when Q 0 = 0 the C 0 j and S 0 j are the principal undepressed particle orbits which generate the matched envelope of the undepressed beam). The IM method applies for extremely strong space charge with σ/σ 0 ≪ 0.1, but probing the limit σ → 0 requires careful analysis of various terms in the formulation presented in Sec. IV.
Complementary to Fig. 3, the decrease in the log of the fractional tolerance [see Eq. (38)] achieved with iteration number is plotted in Fig. 5 for solenoidal and FODO quadrupole focusing lattices for one set of system parameters. The matched envelopes are calculated using the case 2 methods. Case 1 methods and other system parameters yield similar results to those presented. We find that the IM method converges rapidly, with the fractional tolerance achieved increasing by one to two orders of magnitude per iteration till saturating at a value reflecting the precision of numerical calculations employed (∼ 10 −15 fractional accuracy for the examples in Fig. 5).

B. Case 0 parameterization
In parameterization case 0, the matched envelope functions r j are specified by κ j , Q, and ε j . For the ith iteration, the depressed phase advances σ i j needed to calculate the iteration envelope functions r i j are most simply calculated using Eq. (33). Unfortunately, this simple method can fail if space-charge is strong and iterations result in envelope corrections where the radial cross-section of the beam is compressed sufficiently relative to the actual matched envelope solution over the lattice period. Compressive over-corrections can produce iteration principal orbits that are depressed below zero phase advance (i.e., σ i j = 0). In this situation, σ i j becomes complex and we find that Eq. (32) can fail to generate an iteration closer to the desired matched envelope solution.
To better understand where the problem described above can occur, a simple continuous focusing estimate (see Sec. III B) is applied. Taking κ j = (σ 0 /L p ) 2 and ε j = ε, we estimate the envelope compression factor f needed to fully depress particle orbits within the matched envelope. A particle moving within the continuous matched envelope Replacing r b → f r b and σ → 0 in this phase advance formula gives But for continuous focusing, we have [3] σ 0 ε Together, Eqs. (41) and (42) show that f = 0.99, 0.95, and 0.90 (corresponding to ∼ 1%, 5% and 10% compressive over-corrections) will produce fully depressed particle orbits for σ/σ 0 < 0.14, 0.31, and 0.44. Numerically analyzed examples below indicate that this problem can occur in periodic focusing lattices for more moderate space-charge and compression factors than the continuous focusing estimates suggest. The parameter region where the IM method can be applied using the "conventional" case 0 procedure for example periodic solenoid and FODO quadrupole lattices is illustrated in Fig. 6. The region of applicability corresponds to parameters where Eq. (33) can be employed to calculate the iteration depressed phases advances σ i j without obtaining complex values. Iterations necessary to achieve tolerance are plotted as a function of σ 0 and σ/σ 0 . Rather than plotting results in terms of the perveance Q, Eq. (14) was used to calculate σ/σ 0 from the matched envelope functions and system parameters to better quantify the relative space-charge strength where the method fails. Values of Q were chosen to uniformly distribute points in σ/σ 0 . Note that the IM method works with the simple initial seed iteration when space-charge is moderate to weak (0.6 < σ/σ 0 ≤ 1) but abruptly fails with increasing space charge (σ/σ 0 < 0.6). Near the point of failure, convergence becomes slow (iteration counts for the example in Fig. 6 can become thousands if points are chosen sufficiently close to the start of the failure region). Several alternative methods were attempted to render the IM method applicable to all case 0 parameters with arbitrary space-charge strength. We describe these methods without assuming σ 0x = σ 0y and ε x = ε y to better reflect general case 0 applications. First, rather than employing Eq. (33) to calculate the depressed phase advance σ i j of the iteration, the integral formula (14) is applied with the envelope functions of the previous i − 1 iteration with The anticipation is that σ i j calculated from Eq. (43) should be sufficiently close to the actual depressed phase advance σ j of the converged solution to correct the problem. Unfortunately, this method, when applied to example solenoid and quadrupole lattices, results in systematic convergence to unphysical solutions. Replacing Eq. (43) with an "underrelaxed" average over previous iterations might address this problem but was not analyzed. In cases where complex phase advances resulted, various other simple replacements of Eq. (33) were attempted without obtaining satisfactory results.
Several alternative procedures extend applicability to general case 0 parameters. First, slowly increasing the perveance Q from some sufficiently small (or zero) value while implementing the conventional case 0 iteration method using Eq. (33) proves workable in our tests. In this scheme, if Eq. (33) fails (i.e., produces unphysical complex values for σ i j ) then Q is adaptively decreased while iterating until the formula becomes valid before increasing Q again toward the target value. For strong space-charge this procedure can result in many iterations being necessary for convergence because small increases in Q were required in various test cases examined. It is also difficult to determine optimal increments to increase the perveance -which complicates practical code development and can limit the range of method applicability.
Another, simpler to implement, alternative procedure is formulated by combining the Sec. V A method for solving case 1 parameterizations with numerical root finding. In this "hybrid" procedure, the emittances ε j calculated from the x-and y-plane constraint equations (23) are regarded as an undetermined function of the σ j [i.e., ε j | specified = ε j (σ x , σ y )] and trial matched envelope solutions r j are rapidly calculated to tolerance using matched envelopes obtained with case 1 methods for specified (guessed) values of the σ j . Numerical root finding can be employed to refine the guessed values for the σ j to obtain the values of σ j consistent with the target values of ε j . Because the ε j (σ x , σ y ) are smooth, monotonic functions of the σ j for 0 < σ j < σ 0j , the consistent values of the σ j can be found with relatively small numbers of root finding iterations. This is particularly true for plane-symmetric systems (σ 0j = σ 0 and ε j = ε) because one-dimensional root finding can be employed.
The total number of two-dimensional (i.e, the calculations do not assume plane symmetry) iterations needed to implement this hybrid method for case 0 is shown in Fig. 7 for example periodic solenoid and FODO quadrupole lattices. Here, the total iteration number represents the sum of all iterations needed to calculate the emittances to a specified fractional tolerance while calculating all trial matched envelope solutions to a separate specified tolerance over all two-dimensional root finding steps. The same lattices and presentation methods used in Fig. 6 are employed to aid comparisons. Note that the full case 0 parameter space is accessible in this procedure with only relatively modest total iteration counts in spite of the additional numerical work resulting from the root finding. A secant-like multi-dimensional root finding method is employed [16]. Note that only two-dimensional root finding is necessary in contrast to four-dimensional root finding associated with conventional procedures for constructing matched envelope solutions by finding appropriate initial envelope coordinates and angles. Initial root finding iterations are seeded using continuous focusing model estimates for σ j calculated from Eq. (28) using the seed values of r j . Subsequent root finding steps in σ j employ the previous step matched envelopes as a seed envelope in the case 1 iterations. For small root finding steps in σ j this previous step seeding saves considerable numerical work. Only one iteration is necessary for the limit points with σ/σ 0 = 1 because for zero space-charge strength the trial seed iteration is exact to numerical error. Iteration counts at fixed σ 0 likely increase and decrease in σ/σ 0 due to approximate iteration seed guesses being (accidentally) farther and closer to the actual root than in other cases. If the plane symmetries are employed (i.e., using ε x = ε y and σ x = σ y ), then total iterations required can be further reduced. Matched envelopes for general case 0 parameters can also be calculated in similar number of total iterations by analogously combining case 2 methods with numerical root finding. In this case values of σ j consistent with specified values of Q are calculated using the two components of Eq. (23) [i.e., set Q → Q j in the j = x and y components of Eq. (23) and then root solve for σ j consistent with Q = Q j (σ x , σ y )]. (α = 1/2) quadrupole lattices. The same system parameters and presentation format is employed as is used in Fig. 6.

VI. CONCLUSIONS
An iterative matching (IM) method for numerical calculation of the matched beam envelope solutions to the KV equations has been developed. The method is based on orbit consistency conditions between depressed particle orbits within a KV beam distribution and the envelope of orbits making up the distribution. Application of the IM method in simplest form requires numerical solution of linear ordinary differential equations describing principal particle orbits over one lattice period and the calculation of a few axillary integrals over the lattice period. A large basin of convergence enables seeding of the iterations with a simple trial solution that takes into account both the envelope flutter driven by the applied focusing lattice and leading-order space-charge defocusing forces. All cases of envelope parameterizations can be employed, but the method is most naturally expressed, and highly convergent, when employing the depressed particle phase advances σ j as parameters -which also corresponds to a natural choice of parameters to employ for enhanced physics understanding. Virtues of the IM method are: it is straightforward to code and applicable to periodic focusing lattices of arbitrary complexity; it is efficient for arbitrary space-charge intensity; and it works for all physically achievable system parameters -even in bands of parametric envelope instability where conventional matching procedures can fail.