Accurate iterative analytic solution of the Kapchinskij-Vladimirskij equations for the case of a matched beam

The well-known Kapchinskij-Vladimirskij (KV) equations are di ﬃ cult to solve in general, but the problem is simpliﬁed for the matched-beam case with su ﬃ cient symmetry. We show that the interdependence of the two KV equations is eliminated, so that only one needs to be solved—a great simpliﬁcation. We present an iterative method of solution which can potentially yield any desired level of accuracy. The lowest level, the well-known smooth approximation, yields simple, explicit results with good accuracy for weak or moderate focusing ﬁelds. The next level improves the accuracy for high ﬁelds; we previously showed [Part. Accel. 52 , 133 (1996)] how to maintain a simple explicit format for the results. That paper used expansion in a small parameter to obtain results of second-level accuracy. The present paper, using straightforward iteration, obtains equations of ﬁrst, second, and third levels of accuracy. For a periodic lattice with beam matched to lattice, we use the lattice and beam parameters as input and solve for phase advances and envelope functions. We ﬁnd excellent agreement with numerical solutions over a wide range of beam emittances and intensities.


I. INTRODUCTION
This paper analyzes matched beams in alternatinggradient (A-G) focusing systems. Such beams have smaller excursions than mismatched beams, thus requiring smaller transport apertures, and are the starting point for analysis of mismatched beams. Designers and experimenters who work with A-G systems need simple, accurate ways to predict the performance of matched beams. One usually begins with the coupled Kapchinskij-Vladimirskij (KV) equations [1]. For weak focusing fields, solution by the smooth approximation [2], [3], [4] is reasonably accurate. Methods of solution for strong fields [5], [6], [7] have tended to be indirect or complex. Lee, using a double expansion in focus strength and emittance, demonstrated high accuracy for a special case [8]. Our previous paper [9] presented general results in a simple explicit format, but with only moderate accuracy.
The present paper, which was summarized in Ref. [10], uses a convenient iteration method to obtain explicit results. High accuracy is demonstrated for a wide range of parameters, which are: the A-G field strength; beam emittance; and beam charge or current. These would be small parameters if normalized (cf. App. J), but we use physical units throughout. We assume the lattice to be periodic with the beam matched to the lattice, i.e., having the same periodicity. The lattice is assumed symmetric. (Asymmetric cases are treated with a novel numerical technique by Lund et al. in Ref. [11]. ) We solve explicitly for average radius [Eq. (30)], peak radius [Eq. (36)], and the phase advances [Eqs. (40) and (43)]. We give these results for three levels of accuracy and complexity. All these formulas apply to arbitrary symmetric lattices. They require at most the 3rd and 5th harmonics of the lattice profile, becoming simpler for smooth profiles (App. G). Envelope functions are given in Appendix E and illustrated in Fig. 2.
Summary: We begin with the coupled KV equations and show in Sec. III that the matching assumption decouples them so that only one equation has to be solved. (Although the required symmetry Eq. (3f) has been noted before [8], [11], we are not aware that it has been used in obtaining a solution.) Section IV expands the a(z) envelope about its mean and then splits the resulting equation into its average part A and periodic part ρ. The differential equation for ρ(z) is solved (Sec. V) by iteration. These results are combined to obtain a matching equation for the average radius A. This equation is written to various orders of accuracy; our special definition of "order" is given in Sec. V. The results (Sec. VI) are compared to the results of numerical solution of the full KV equations. The first-order case is usually called the smooth approximation. Second-and third-order terms increasingly improve the accuracy. Over a wide range of parameters, our third-order versions provide a great improvement over those previously published [2]- [6].
Section VII combines results from Secs. V and VI to give the maximum and minimum radii. The phase advances σ and σ 0 are given in Sec. VIII. There is a useful approximation for σ 0 in Sec. IX. Appendices F, G, and H discuss the use of Fourier analysis of the A-G focusing function to facilitate the solution. Some exact formulas for the FODO case are given in Appendix I.

II. THE KV EQUATIONS AND SYMMETRIC LATTICE MODEL
The KV equations for envelopes a(z) and b(z) are [1] with dimensionless perveance Q, emittance ∈, and lattice focus-strength function K(z). The latter is periodic over the lattice cell length, defined as 2L, so that K(z + 2L) = K(z). In practice, the lattice cell has one or two symmetries. The double symmetry occurs if the focus drift spaces are equal in length. In this paper we assume that the function K(z) has odd symmetry at the center of the drift spaces and even symmetry at midelectrode points. These symmetry points necessarily alternate at intervals L/2. It is convenient in what follows to choose one of the even points for the z origin. (The odd choice is appropriate for the case of unequal drift lengths; that case is not discussed in the present paper.) Then the double symmetry is expressed as (3b) It follows that K(z) is anti-symmetric about 3L/2 and that K(z) obeys the shift relations K(z + 2L) = +K(z).
(3d) It also follows that Equation (3d) confirms the lattice periodicity which was assumed above. Equations (3a) and (3b) indicate that K(z) can be expanded as a cosine series containing only odd-numbered harmonics.
In this paper we solve Eqs. (1) and (2) for the beam envelopes in the special case where the beam is matched, i.e., a(z) and b(z) have the same periodicity as the lattice. The initial conditions with the symmetry of Eq. (3e) are For given Q and ∈ and given lattice parameters, the beam will only be matched for specific values of a 0 and b 0 , which we denote by a m and b m . Then the initial conditions for a matched beam are (1) and (2) and use Eq. (3c), we find that a L (z) and b L (z) also satisfy the KV equations with transposed initial conditions: so that the denominator of the last term in Eq. (1) can be written a(z) + b(z) = a(z) + a(z + L). Then b(z) does not appear and Eq. (1) is decoupled. In Sec. III and App. A we find an expression for the function [a(z) + a(z + L)] −1 in terms of integrals involving the given lattice function K(z). After that, Eq. (2) is not needed for our matchedbeam case. The main work of this paper will be to find a m and b m -see Secs. IV through VII. First we introduce some definitions.

Operators, Functions, Parameters, Derived Quantities
To aid the solution of Eqs. (1) and (2), we define in Table I {f } ≡ f −˙f¸.
For even ψ(z) ψ = 0: and A ≡˙a(z)¸, In Table I, the operator . . . performs an average over a cell length 2L. The operator {. . . } removes the average part of a periodic function: e.g., 2 cos 2 x = {1 + cos 2x} = cos 2x. The operator operates on even periodic functions that have zero average. It gives the repeated indefinite integral and removes the average part, if any, of the result. This removal can be implemented by constructing a suitable lower limit for the outer integral. To illustrate, Eq. (10) can be written which subtracts the value at L/2, so that g(L/2) = 0. Alternatively, one can start both integrals at zero and then apply the operator {. . . }, as in Eq. (7). For example, cos z = {sin z} = {1 − cos z} = − cos z. Note that, from Eq. (9), h(0) = 1.

III. DECOUPLING THE KV EQUATIONS FOR SYMMETRIC MATCHED-BEAM INITIAL CONDITIONS
For a matched beam with the symmetries of Eq. (3), we showed in Sec. II that b(z) = a(z + L), implying that Eqs. (1) and (2) are decoupled and that a = b ≡ A. We write with ρ(z) > −1 for all z. Then the Q term in Eq. (1) can be written where F is obtained from an expansion and iterations. All terms consist of integrations of the focus-strength profile h(z) [Appendix A]: Additional terms are shown in the appendix. With the lattice symmetries of Eq. (3), all the series terms for F (z), including δ(z), have only even-numbered cosine harmonics-in contrast to h(z), which has only odd ones.
Equations (1) and (2) are decoupled to all orders so that Eq. (2) is superfluous from here on. Equations (1) and (2) are replaced by

IV. EXPANDING AND DECOMPOSING INTO AVERAGE AND PERIODIC PARTS
Substituting a = A 1 + ρ(z) in the first three terms of Eq. (1), expanding 1/a 3 , dividing by A, and using (22) and (14), Eq. (1*) is equivalent to To solve for the ripple ρ(z) and for the mean radius A (which appears in the definition of α and q), we decompose Eq. (23) into a pair of equations. Averaging Eq. (23), (24) Subtracting Eq. (24) from (23), with the {· · · } operator defined by Eq. (5). There are now two equations, each containing ρ(z) and A-the latter represented by α. These represent the KV equations (1) and (2) for our matched beam. We will obtain A and ρ(z) to good accuracy in Secs. IV-VI and App. E.

V. ITERATIVE SOLUTION: ARBITRARY SYMMETRIC FOCUSING PROFILE
On the right-hand side of Eq. (25), the kh(z) term dominates the terms involving the unknown function ρ(z). Therefore, we use only kh(z) in the initial integrations which give ρ (0) . Then we insert ρ (0) into (25) and integrate again to obtain ρ (1) . (A miniscule term, qk 2 δ(z), in ρ (0) is omitted.) This process is repeated to get ρ (2) : In (26c), stands for two applications of the operator defined by Eq. (7).
To complete the approximate solution of the KV equations, ρ(z) from Eq. (26) is put into the matching equation (24). In Eq. (26) we included no items e.g., 2αk 2 {g} 2 that would give terms in (24) higher than third power in the parameters k 2 , α, and q. That is, we go no higher than third order as defined below.
Some terms vanish by orthogonality, since h(z), g(z), g 3 , etc., possess only odd harmonics while δ(z), g 2 , g 4 have only even ones for symmetric quadrupoles. A term not shown, qk 2 δ(z), involves multiple integrations of an already small function and would contribute < 0.04% to the maximum radius a max even at σ 0 = 120 • and affect A by less than 2 parts in 10,000. This is much smaller than the two parts per thousand criterion used for significant terms in Appendix E.
Order: We reckon the order of a term in the matching equation by counting the number of factors k 2 , α, and q. These would become small parameters in a nondimensional formalism (Appendix J). We prefer to retain physical units for quantities such as the axial coordinate z. Note that Lee [8] defines order using k and ∈ so that our third order corresponds approximately to his fifth order. Appendix B evaluates the combination of Eqs. (24) and (26). The result, in simplified form, is shown in the next section. Input quantities: Q, ∈, and quad voltage VQ(∝K). Other parameters: see Table II. VQ, fixed at 20 kV, gives phase advance σ0 of 83.37 • ; ∈ and Q are varied so that depressed tune σ ranges between 0 • and 76.5 • ; exact σ0 and σ are obtained numerically. where In Eq. (28) (27) is in the standard form of the smooth approximation, Eq. (33), and can be solved to find the third-order A: If the input quantity is the mean radius A inp , then Eq. (27) gives the allowable Q to third order, Second Order: Eq. (B7) has two second-order terms. One yields the correction to K eff seen in Eq. (28). The other term is αk 2 g 2 , or, using definition (18), αΦ/3. We define and get Eq. (32) can be solved for A II or Q II in the same way as for the third order, giving useful approximations when K(z) and ∈ produce σ 0 and σ less than about 80 • .
First Order: In Eq. (B7), the three terms of lowest order in α, q, k 2 produce what is called the first-order matching equation in this paper (Ref. [9] used another terminology). This is the classic smooth approximation. These terms give k 2 [ h] 2 = α/3 + q, or, using the definitions (14), (16), and (17), The equations derived in this paper apply to arbitrary focusing profiles satisfying Eq. (3), such as the smooth profiles in Appendix G, but for illustration (Figs. 1, 2, and 3), we use the FODO model (Appendix H) with parameters from Table II. First-, second-, and third-order results for A, from (33), (32) and (30), are plotted in Fig. 1a for FODO. The smooth approximation is seen to be relatively inaccurate except near the point where its error curve crosses the 0% line.

VII. EXPLICIT THIRD-ORDER RESULT FOR a max
Knowing the matched mean radius A, one can complete the solution for the beam envelope a(z) = A 1 + ρ(z) using ρ(z) from Eq. (26); b(z) can be found by changing the sign of the terms that contain odd powers of k.
Some terms of Eq. (26) can be written in exact form [Appendix I] for models such as FODO, but Fourier expansion is more useful in general: (34) Values (usually of order unity) of h 1 and c n for both smooth and FODO profiles are given in Appendices G and H. With the definition we have  (37) and (the smooth approximation) The time dependence of Eq. (36), from App. E, is plotted in Fig. 2

VIII. PHASE ADVANCES: ARBITRARY SYMMETRIC FOCUSING FUNCTIONS
From the well-known phase-amplitude result [13], the phase advance per quadrupole cell of length 2L is We approximate a(z) by A III [1 + ρ(z)] with A III from Eq. (30) and ρ(z) to third order from Eq. (26). Subscripts  are omitted for brevity. Expanding a −2 and taking the average gives (The 2ρ term has zero average by definition.) Appendix D shows that to third-order accuracy for arbitrary symmetric A-G lattice functions. Errors with respect to exact values from simulations are illustrated in Fig. 3a for the particular case of a FODO lattice. Useful accuracy is retained after dropping two terms and using lower-order A II from Eq. (32): Figure 3a shows large errors for the first-order result (smooth approximation): The undepressed tune σ 0 for arbitrary symmetric focusing profiles is found by setting Q = 0 in Eq. (27), then eliminating ∈ from Eq. (40). Details are in Appendix D.
The third-order result is  for arbitrary symmetric focusing functions. To secondorder, The smooth approximation is Figure 3 compares the accuracy of Eqs. (43), (44), and (45) for the special case of FODO focusing (with occupancy η = 0.5), but gives an idea of the relative accuracy for arbitrary cases. At σ 0 = 100 • , the smooth approximation is off by 12%; this is improved to 1% for third order.
We will find that only lowest order results are needed and drop the R 2 term, etc. We subtract the average of Eq. (A1), using definitions from Table I, and obtain which we solve by iteration. In the initial step, the only known quantity is kh(z), so that Substituting −kg for ρ in the previous equation, where we used the definition of δ. Using the shift Eq. (3c), Then Eqs. (A2), (A6) and (A7) give resulting in Eq. (1*) in Sec. III. The expansion and iteration could be extended to produce more terms [all derived from h(z)] but these would give even smaller corrections to our results.

APPENDIX B: DETAILS OF DERIVATION OF THE MATCHING EQUATION
Rearranging Eq. (24) slightly gives (B1) Inserting ρ from Eq. (26), the left-hand side of this equation is where we have dropped the subscript on ρ. The orthogonal k 2 δ term is absent. We simplify by changing the order of integrations, using the h(z) symmetries [Eq. (3)]. For example, − h g = − g g = + [ g] 2 . Applying this technique throughout gives For the right-hand side of (B1), The very small k 4 δ 2 term was omitted from ρ 2 .
Again changing the order of integrations, the right-hand side of (B1) becomes Four of the terms of (B6) combine with terms of (B2), so that The three k 4 terms on the rhs of Eq. (B7) are: Adding (C3) through (C5) gives 3 2 g 2 2 (1 + (29/243)c 3 + · · · ), where the small c 3 correction can be neglected since it corrects a term which is already third order. Using all these results along with definition (14), the right-hand side of (B7) (without the q term) becomes as in the main text. In the last term, uses the lowest-order value for A because this term is already of the highest order that we retain.

APPENDIX D: DEPRESSED AND UNDEPRESSED TUNES, GENERAL CASE
Here we evaluate the expansion terms in Eq. (39). From Eqs. (B3-B5) in App. B, From App. C, Eqs. (C2)-(C5), When the last three are added, the c 3 coefficient is only 25/729. Dropping this and using the definitions gives to third-order accuracy Undepressed tune: Setting Q = 0, (27) is The factor ∈ 2 /AI 4 in the last term can be replaced by Comparing with the definition of Φ and Fourier expanding as before, the last term, to lowest order, is 3Φ 2 for Q = 0. Thus, altogether, Making a similar replacement in Eq. (40) for the case Q = 0, Using Eq. (D10) to eliminate ∈/A III 2 , or, finally, to third order.
for arbitrary symmetric A-G profiles. Sometimes it is convenient to work with the squares of σ 0 and σ, which are for third order and Using Fourier expansion, written as in Appendix C, and recalling ρ m ≡ h 1 kL 2 /π 2 , the terms of Eq. (26) are The small final two terms from Eq. (26) have been omitted here for simplicity.

The Significant Terms
We drop small quantities in the above equations. The criterion is that they contribute less than two parts per thousand to the final result for a bad-case scenario: large focusing voltage (giving phase advance of 112 • ) and large β. This leaves Adding all these gives ρ(z) to third-order accuracy. A few small terms were omitted as mentioned before. Setting z = 0 gives ρ max and a max = A (1 + ρ max ) as presented in Section VII. Setting z = L changes the sign of all terms, except the even term k 2 δ(z), and yields a min . α and β are only needed to lowest order, i.e., α I and β I ; we omit the subscripts here. Results from Eqs. (E6)-(E10) are shown in Fig. 2 in the main text.

Fourier Coefficients
Recall from Section II that the focusing force K(z) in the KV equations is written as The Fourier coefficients are It is often convenient to define where c 1 = 1 by definition and where |c 3 | usually turns out to be of order unity-see Tables III and IV

Solution of Envelope Equation
In the solution for ρ(z), Eq. (26), the largest term is The mean square of the integral of Eq. (F2) gives the effective force which is used in the matching equation and for calculating undepressed phase advance. The correction term Φ (used in evaluating phase advances, average radius or transportable current, etc.) is since the harmonics contribute practically nothing. Dividing this into Eq. (9), we find which could be useful in certain calculations.
APPENDIX G: SOLUTION FOR SOME SMOOTH LATTICE PROFILES All our results apply to an arbitrary symmetric focusing profile once its Fourier coefficients are known. Some simple but representative smooth models are discussed here. In each case, we use only the lowest required harmonics.
This is the model used by Lee [8].
with all the other coefficients zero except c 1 , which is unity by definition. This choice gives a flat field at the midpoint of the quadrupoles, without the discontinuities The third-harmonic correction for Φ can be neglected in most cases.
which gives a fairly realistic profile (Fig. 4) and corresponds to FODO with η ∼ 53%. The third-and fifthharmonic corrections are well under 1% for this case. Table III summarizes the above results.  In the case of the popular FODO lattice model (Fig. 5), the Fourier coefficients h n are readily calculated from Eq. (F4):  Fig. 5.) Because of the symmetries expressed in Eq. (3), the calculation of averages is simplified, requiring only integration over one-fourth of a cell. One finds for a specified One can do the integrals in the first two terms of ρ(z), Eq. (26). For the integrations, we divide the cell into five zones: For even-numbered zones, the first integral is (I3) with n the zone number. P is defined as follows: for any integer m, For odd-numbered zones the integral is − h = P n + 1 2 As required by the definition of , the average has been subtracted. The maximum value of ρ (0) (where z = 0, implying n = 0) is The next term in Eq. (26) includes the integral g. For even-numbered zones and for odd-numbered zones The maximum value of g is These results, for any value of η, may be compared with those from Appendices G and H to determine the number of Fourier terms needed for a given accuracy in each case.

APPENDIX J: NON-DIMENSIONAL FORMALISM
As mentioned in Section V, we have chosen in this paper to use physical units for quantities such as the axial distance z. The result is that the expansion parameters k 2 , α, and q are not necessarily small numbers. To show that they are nevertheless appropriate parameters, we replace z by The KV equations become where the primes now indicate differentiation with respect to θ.
In Table I, averages are now defined by Similarly, the symbols and refer to integration with respect to θ, not z. Other replacements in Table I are All these items are now dimensionless. In Eq. (15), the quantity β is supplanted by α. Equations (35) and (C7) no longer apply; β must be replaced by α in Eqs. (29), (36), (37), (40), and (C6) as well as four places in App. D and six places in App. E. The equations throughout this paper are otherwise unchanged, and one can verify that the expansion parameters k 2 , α, and q are now dimensionless small numbers.

APPENDIX K: ALTERNATIVE MATCHING AND INITIAL CONDITIONS
Another method of beam matching is convenient for some simulation studies. The beam is launched at the midpoint of a quadrupole gap. Using the same lattice models as before, the matching conditions would become This choice would be essential for the case of unequal drift spaces (only a single symmetry) but we do not consider this case here.

APPENDIX L: WORK FOR THE FUTURE
It would be useful to investigate the modifications that would be required to extend the methods of this paper to non-symmetric A-G profiles, which are often used in practice.