Time evolution of coupled-bunch modes from beta function variation in storage rings

We present an analytical and numerical study of the equations of motion for bunches coupled by transverse wakefields. We base our study on a recent lattice design for the damping rings in the baseline configuration of the International Linear Collider. Using the macroparticle model, and assuming resistive wall wakefield coupling, we present numerical results on the time evolution of the multibunch modes. Decay modes display growth after initial decay, and mode amplitudes exhibit high-frequency oscillations. These phenomena are not expected if the beta function is assumed to have a constant, averaged value. We show analytically that they can come from coupling between modes caused by variation of the beta function in a real lattice. The effect is shown to be comparable to the effect of a nonuniform fill pattern and significantly larger than that of the higher-order mode wakefield localized in the rf cavities. Turning to the case of constant beta function, we develop a more complete treatment of the equations of motion. We derive general formulas for the bunch trajectories, and show that such formulas can only be valid in the limit of small wakefield coupling.


I. INTRODUCTION
Wakefield coupling between bunches is an important cause of instability in accelerators [1][2][3]. This is a widely studied area and many numerical studies have been carried out [4 -6]. The electromagnetic fields induced by relativistic charged particles traveling through an accelerator beam pipe affect the motion of the following particles. In a synchrotron storage ring under certain conditions, the wakefields acting from one bunch to another can lead to beam instability, manifested as motion with large synchrotron or betatron amplitude, decoherence, and beam loss. Feedback systems are commonly installed in modern, high intensity storage rings to maintain bunches on the design trajectory by providing precisely timed kicks [7,8]. In the case of the damping rings of a linear collider, strong wake forces can occur during the injection process because of the large offsets of these bunches from the design trajectory [4,6], and there are concerns that these forces could induce jitter in the damped bunches ready for extraction.
A good understanding of the behavior of wakefield coupled bunches is therefore important for storage rings, and models have been developed to describe the observed dynamics. Various analytical results based on the macroparticle model are available for the equations of motion of wakefield coupled bunches. For example, there exist analytic formulas for the growth rates that describe the increase in the amplitudes of the normal modes resulting from wakefield coupling [5]. Another example is analytic formulas to describe the time evolution of each bunch [4]. These results are important for the design of feedback systems needed to suppress any instabilities. The results may also be useful for predicting the transient effects that can arise from injection and extraction of bunches.
The concept of normal modes is convenient for characterizing the coherent motion of bunches in a storage ring. The mode amplitudes can be obtained by performing a discrete Fourier transform (DFT) of the transverse or longitudinal bunch positions at any instant. The resulting snapshot, repeated over time, provides a useful description of the evolution of the beam in a real accelerator [3], and is an important means of characterizing the effectiveness of the feedback control system.
The derivation of the analytic formulas referred to above, and the use of the normal modes, are based on the assumption that any effects from variations of the beta function around the ring average out, so that the beta function can be assumed to be constant. The agreement with experimental results and the successful design of feedback systems show that this has indeed been a reasonable approximation in machines built and operated to date. However, the demand for very high levels of stability in machines with intense beams, such as linear collider damping rings [9], make it worthwhile to study this issue in greater detail, both to investigate the effects of variations in the beta function, and to examine numerical accuracy of the analytic results.
The baseline configuration for the International Linear Collider (ILC) [9] specifies damping rings with a circumference of over 6 km, storing between 3000 and 6000 bunches with 400 mA average current [10]. Initial estimates suggest that growth rates from resistive wall wakefields alone could drive coupled-bunch instabilities with growth times of the order of a few tens of turns; only just within the capabilities of modern digital feedback systems. Long-range wakefields from other sources, such as higher-order modes in the rf cavities, could increase the growth rates further. Furthermore, large variations in beam current (of up to 10% of the nominal stored current) will occur in the positron damping ring during the process of extraction and reinjection. Freshly injected bunches can be expected to have significant longitudinal and transverse jitter, which may couple to damped bunches awaiting extraction. A thorough understanding of the dynamics resulting from wakefield coupling between bunches is necessary if bunches meeting the stability specifications are to be extracted from the damping rings. To address these issues, we base our present study on the macroparticle model, which has been widely used and shown to produce reasonable agreement with experiments. We carry out detailed tracking simulations to include the effects of variations in the beta function through a realistic lattice, based on the present ILC baseline design, and reexamine the assumptions of previous analytic results (e.g. Ref. [3]). To understand all the features observed in the tracking simulations, some extension to the existing theory is needed.
The structure of the paper is as follows. In Sec. II, we describe the tracking simulations, discussing the model on which it is based and significant features in the results when applied to the ILC damping rings. In Sec. III, we consider in more detail the mode coupling and the highfrequency oscillation of the mode amplitudes that are observed in the simulation results, and show that these effects may be the result of variations in the beta function around the ring. In Sec. IV, we revisit the analytical solution to the equation of motion for wakefield coupled bunches in a storage ring, considering in some detail two aspects that have not previously been fully discussed: the significance of the ''initial history'' of the system; and the existence of a family of solutions to the equations of motion, beyond the single solution usually identified with the instability growth rate. Finally, in Sec. V, we derive a solution to the equations of motion that is more general than the one usually written for wakefield coupled bunches. Partly for convenience and partly to clarify our general approach, some standard material is summarized in two appendices. Appendix A briefly reviews the model leading to the familiar formula for the growth rates of coupled-bunch modes. Appendix B rederives the kick method used in the tracking simulation by direct integration of the equations and clarify the assumptions made. Comparisons between the beta function variation and some other effects that can impact the growth rates are given in two further appendices. Appendix C reexpresses the formalism in terms of matrices for linear systems, and compares the effect of varying beta function with the effect of an uneven fill pattern. Appendix D calculates the effect of the transverse higher-order mode (HOM) from the rf cavities and shows that, with the kind of cavity expected for the ILC damping rings, the effect is small. Finally, Appendix E repeats the simulation for a simple lattice to show that the effects of beta function variation is not limited to large rings.

II. TRACKING SIMULATION
In earlier papers on coupled-bunch instabilities, the beta function has sometimes been assumed to vary only slightly from its average value [11]. Even when it is known to vary strongly, it has been hoped that the rapid oscillations of the beta function can be averaged out when assessing their effect on the beam [4]. In Ref. [6], for instance, although a tracking simulation has been carried out on a lattice with large variations in the beta function, the results have been compared with analytical results obtained from a model that assumes a constant beta function. The fact that there appears to be a good agreement between the simulation and analytical results justifies the averaging of the beta function in the analytical model.
Various methods are available for numerical integration of the equations of motion, given by and explained further in Appendix A. In this paper, we assume a form of the resistive wall wakefield given by Eq. (A3). More refined treatments are available -see [12] and references therein. The qualitative features of the results presented in this paper are likely to remain unchanged because the wakefield model only determines the constant numerical values of the coefficients of x in Eq. (1), and does not affect the form of the differential equation. We hope to repeat the calculation with improved wakefield models in the future. If the wake force on the right of Eq. (1) is set to zero and the beta function is averaged around the ring, we have the equation of motion for a simple harmonic oscillator, which has the simple solution Constants A and B are obtained from initial conditions at t 0, where the values x m 0 and _ x m 0 are specified. Reference [4] shows that, for small wake forces, it is sufficient to compute an instantaneous change in momentum, or ''kick,'' from the wake forces once every turn around the damping ring. The kick is computed by taking the product of the wake force at a fixed point on the ring and the time of revolution, T 0 . The kick causes a change in _ x m t but not x m t. The new value of _ x m t is then used to define the starting point for the computation of the trajectory x m t over the next turn. This process is then repeated over many turns.
In Ref. [6], action-angle variables are used. The kick is applied to M points around the ring, where M is equal to the number of bunches. Each point is called a time slice. At each time slice, the bunch is given a kick. The Cartesian coordinates x m t and _ x m t are converted into action-angle variables. The wake force is assumed not to act between slices. The action remains constant up to the next slice, while the angle increases by an amount equal to the phase advance. At the next slice, the Cartesian coordinates of the bunch are calculated from the action-angle variables, and the next kick is applied. Here the kick is the product of the wake force and the time taken for the bunch to travel the distance between the slices. The process is repeated for each subsequent slice. This method has the advantage that, if the Twiss parameters are known over the lattice, the effect of the varying beta function can easily be included through the phase advance, and through the conversion between Cartesian coordinates and action-angle variables.
The methods of delay differential equations (DDEs) [13][14][15][16] offer some possibilities for a more careful approach. DDEs are equations which involve the history of the variables, which is the case for the equation of motion we consider here, describing the motion of bunches in a storage ring in the presence of long-range wakefields. The general approach is to integrate a DDE numerically or analytically in the way that would normally be done for ordinary differential equations. For example, Eq. (1) would first be integrated over one time step, . In order to do so, the value of x m t must be known over the time interval. In the case of numerical integration, however, this is normally known over discrete points on the interval. Except for the simple Euler method, this is not sufficient, either for a Runge-Kutta integration or an analytical integration. It is therefore necessary to interpolate x m t over the time interval: the quality of this interpolation is crucial to the accuracy of the integration [15]. Many DDE algorithms have been developed to address this issue.
Let us now apply this approach to Eq. (1). In our case, the distance between time slices is small compared with the wavelength of the betatron oscillation: the horizontal tune x is about 52.4, whereas the number of time slices we use is M 3649. This gives M= x 69:6 time slices per oscillation. If we imagine dividing one wavelength of a sine curve into this number of intervals, we see that a linear interpolation should be quite sufficient. If we apply this to the wake force term, we can replace it completely by a single linear term pt q fitted over each time interval. As it turns out, this form is amenable to a direct analytic integration. The surprising outcome is that the solution resulting from this first-order interpolation of the wake force is exactly equivalent to the kick method using a zeroth-order interpolation of the wake force described above. The full derivation is given in Appendix B. This approach puts the heuristic kick method on firm mathematical grounds, where the underlying assumption is now properly clarified: the assumption is simply that a linear interpolation on the history has been used. This allows us to keep track of the order of accuracy involved. For in-stance, we can imagine improving the accuracy by using a more accurate interpolation. Furthermore, provided that the linear interpolation on the history is sufficiently accurate, the requirement that the wake force be small is no longer necessary.
We now apply the kick method to the ILC damping ring, for the horizontal motion. Neglecting effects resulting from energy variations due to changes in path length, the horizontal and vertical dynamics are essentially the same. The parameters used are as follows: We first carry out the simulation for the case when the beta function is constant, and compare the growth rates with the analytical result. This provides a reference for the case that includes variations in the beta function. The simulation in this section follows closely the method in [6] : (i) The bunches are regularly spaced, and travel very close to the speed of light. (ii) The initial x m t and _ x m t at t 0 are given random values. (iii) The history, or wakefield sum on the right of Eq. (1), is truncated at N 100 bunches. (iv) The initial history, i.e., the values of x m t and _ x m t over the time interval ÿN; 0, are set to zero.
(v) The distance between time slices is taken to be the distance between bunches. (vi) At each time step, the action-angle variables, J m and m , are obtained by transforming from x m t and _ x m t for each bunch. (vii) The mode is obtained by a discrete Fourier transform over the normalized action, 2J m p e i m , with respect to m. The original idea behind using the normalized action in [6] is to include the velocity information in the mode calculation. This turns out to have interesting consequences which will be discussed in Sec. IV. Here we just note that the mode, defined in Appendix A as a form of DFT over x m t, is not exactly the same as a DFT over 2J m p e i m . In [6], the growth rate is obtained for each mode by plotting the mode amplitude against time, and fitting the graph with an exponential function. There, it is observed that some of the modes are not simple exponentials. We shall later demonstrate a few situations where such nonexponential behavior may appear; but for the moment we simply assume that the modes evolve exponentially, as suggested by the analytical procedure leading to the growth rate formula [5]. Then, if a mode has frequency ! i= , the evolution of the mode has the simple exponential form e ÿi t e t= e ÿi! t . The derivative of the mode amplitude e t= at t 0 is 1= , which is just the growth rate. We therefore obtain the growth rate by computing the initial gradient of the mode amplitude -over the first turn. Figure 1 shows the simulated growth rates plotted against mode number, for the case of constant beta function. The growth rates calculated from the analytical formula in Appendix A are plotted on the same graph. The agreement is good, which seems to justify our assumptions so far. In Sec. IV, we shall show that the ''oscillations'' in the simulation result for the growth rate as a function of mode number are due to truncation of the wakefield sum in Eq. (1). Note that radiation damping is not considered in this paper. The way to include this is explained in Appendix A.
We now apply the simulation to the varying beta function, with N 100. Since the output positions from MAD [17] are not regular, we need to interpolate to obtain the values at equally spaced intervals, corresponding to the distance between adjacent bunches. The interpolation method we use is to fit a cubic equation to the beta function between every two adjacent time slices, with the constraint that the gradient of the fitted cubic matches the gradient of the beta function (obtained from the alpha function output) at each time slice.
The beta function clearly has large oscillations. In the absence of any wakefield, the position of a bunch is given by In our case, the wake force is weak compared to the focusing force. If we consider the wake force to be a perturbation to the betatron motion, the mode amplitude should follow closely the trends of the beta function, and oscillate strongly. Figure 2 shows the mode amplitude of mode number 100 plotted against time. As expected, it oscillates strongly. This immediately presents a problem when we want to determine the growth rate, since the initial gradient is no longer equal to the growth rate because of these oscillations. If we suppose that the long term, averaged behavior is still exponential, we may try to extract a growth rate by smoothing out the oscillations before taking the initial gradient. We do so by taking the average of the amplitude from turn 0 to 0.1, and the average over turn 1 to 1.1. These averaged amplitudes are used to determine the gradient over one turn, which is taken as the growth rate. The result is shown in Fig. 3. There is a significant increase in scatter compared with Fig. 1, but the results still follow the trends of the analytic result fairly closely. This appears to support the idea that the effect of strong oscillations in the beta function can be averaged out.
However, there are a number of modes that show radically different behavior between the case in which we use an averaged beta function, and the case in which we use the real variation of the beta function around the lattice. For example, if we look at mode 3600, its growth rate in the case of an averaged beta function (Fig. 3) is negative, so it should decay exponentially. However, when using the real beta function, the growth rate appears to be positive (Fig. 3). Does this mean that mode 100 would grow? To answer this question, we simulate the long-term behavior of the mode amplitude. Figure 4 shows the results. The amplitudes are displayed every 10 turns, for 0.1 turn. The amplitude in the first 0.1 turn has been expanded in Fig. 2.  Mode 100 decays initially, but starts to grow after 40 turns. This is surprising, because a mode that evolves with simple exponential behavior should either always decay, or always grow. An inspection of the other ''decaying modes''those with negative growth rates as shown in Fig. 3shows that all of them grow after about 100 turns. For comparison, we have included the result for nonuni-form charge distribution (uneven fill) in the case of constant beta function, as well as for the transverse HOM in the rf cavities. These are discussed in more detail in Appendices C and D.
Let us turn our attention to another feature in this simulation. Figure 5 shows three samples of mode amplitudes, taken at turn numbers 0, 60, and 90. Superimposed upon the variation of the mode amplitude is a small oscillation with a much higher frequency compared to the oscillation that results from variation in the beta function.    This oscillation becomes increasingly obvious as the turn number increases. To search for the cause of this oscillation, we compute the frequency spectrum of each mode at turn number 90, where the oscillation is most prominent. Specifically, we perform a fast-Fourier-transform of the normalized action from turn 90 to 90.1. This contains 365 time slices or data points. The result is displayed as a color scale plot in Fig. 6. The maximum height of the spectrum for each mode has been normalized to 1. The betatron oscillation shows up as a horizontal line near the top. The high-frequency oscillation falls exactly on the diagonal through the origin. In the following section, we shall seek an explanation for the two unexpected features observed in the simulation results: the long-term growth of modes that are expected to decay, and the high-frequency oscillation in the mode amplitudes.

III. DYNAMICAL EFFECTS OF VARIATIONS IN THE BETA FUNCTION
In this section, we consider the growth of ''decay'' modes and the high-frequency oscillations that are observed in the tracking simulations of the ILC damping ring in the previous section. These features are not expected from the standard analytical theory, and in order to explain them, we need to make some development of the theory. We shall focus our attention on the effects of the strong oscillation in the beta function. The original equation of motion (1) can be written in the following form: for bunches m 0 to M ÿ 1, where the constant factor and wakefield coefficients have been absorbed into the constants b n . We multiply this by expÿi2m=M, and sum over m. After some manipulations, we get the following set of M equations: b n e i2n=Mx t ÿ n: (5) This is the equation for mode , which is defined by These equations are completely decoupled from one another. This derivation is exact as long as the beta function, and hence the betatron frequency ! , is constant.
To include the effects of variations in the beta function, ! is treated as a function of time. First, define So that the equation of motion can be written: b n x mn t ÿ n: (8) Notice that the index m enters into the argument of Kt because each bunch sees a different beta function. Let K be the average of Kt over time, and define the deviation of Kt from K: kt Kt ÿ K: (9) Then move the oscillating part of Kt to the right-hand side in Eq. (16): Now transform to modes. Multiplying by expÿi2m=M and summing over m, we obtain Compared with Eq. (5), a new term has appeared in the equation for the modex t. This expresses the oscillating part of Kt as a driving force. Another form for this term can be obtained by defining the following: This is essentially a discrete Fourier transform. Applying the convolution theorem, we get The new term now expresses the coupling between different modes. The driving force term that we derived above, ÿ P Mÿ1 m0 e ÿi2m=M kt mx m t, arises as a result of variations in the parameter K. For convenience, we shall call it the parametric driving force, after the concept of parametric resonance [18]. If the parametric driving force is responsible for the high-frequency oscillations observed in the tracking simulations, then this driving force should oscillate at the same frequency. We can show that this is indeed the case, by the following argument. Consider a hypothetical situation in which all bunches are at a constant unit displacement, i.e. x m t 1 for all m. The parametric driving force is then ÿ P Mÿ1 m0 e ÿi2m=M kt m, which is just ÿk t. This is essentially the spectrum of kt. From the property of Fourier transforms, we know that a shift in time corresponds to a shift in phase in the frequency domain. We can therefore expressk t in terms ofk 0. For our purpose, it is sufficient to sample the time t at the time slices defined around the damping ring. So let t p, where p is an integer. Substituting into the parametric force term, kp m: Since km repeats itself around the damping ring, it is periodic in m with period M. Therefore, Since t p, we can also write ÿk t ÿe i2t=Mk 0: Thus the parametric force has an oscillation frequency of f M . This is the same frequency as that of the highfrequency oscillation observed in the simulations: if we were to plot f against the mode number , we would see exactly the same diagonal line that we see in Fig. 6.
We thus have a plausible explanation for the highfrequency oscillation. As the bunches travel round the ring, they are perturbed by the oscillating part of Kt. Each mode has a characteristic frequency. This is also given by f M , because is the number of wavelengths present in that mode, and M is the period of revolution around the ring. Each mode would therefore resonate with the corresponding frequency component of Kt, which is k 0. For convenience, we shall call this oscillation the parametric oscillation.
We now turn to the mode coupling term in Eq. (13) and consider the possibility that this may lead to exponential growth of decaying modes. Consider two modes, a growing modex 00 and a decaying modex . As a result of the mode coupling term, the growing mode would make a contribution 1 Mk ÿ 00x 00 to the equation for the decaying modex . In the long term, this could cause the decaying mode to grow. To demonstrate this mode mixing, we carry out a simulation for the case where only one single mode, , is present. The purpose is to see if the other modes would in fact appear after some time. In order to do so, we must construct the bunch displacements for a single mode, and use this as the initial condition at t 0. This is given by [5] x m t e 2m=M e ÿi t : The time derivative is then given by _ x m t ÿi e 2m=M e ÿi t : Setting t 0 and taking the real parts, we obtain the required initial conditions: We carry out the simulation for mode 500. Figure 7 shows mode mixing that takes place when the actual ILC damping ring beta function is used. Notice that all the other modes start to appear after a fraction of a turn. We have also carried out the simulation result when Kt (or the beta function) is constant. As expected, there is only ever one mode as time evolves-the spectrum of mode amplitudes remain exactly the same as the top panel in Fig. 7  From the results in this section, we conclude that, when the beta function varies, the Fourier mode is no longer an eigenmode. An eigenmode may vary in amplitude as time passes, but the amplitude would scale by the same amount at every point in space. As Fig. 7 shows, a Fourier mode can change into something quite different if the beta function is not constant. We note here that other disturbances, such as nonuniform charge distribution and localized wakefields, can also have the same effect. We consider these cases in Appendix C and D, respectively, and compare their effects with the varying beta function. The above effects are unexpected and have not previously been observed. In order to confirm that they arise from a simple variation in the beta function and are not due to the complexity of the ILC damping ring, we repeat the simulation for a simple lattice in Appendix E.

IV. THE INITIAL HISTORY
We now turn to the possibility of deriving an analytic solution for the equation of motion, at least when the beta function is constant. To do so, we must address the issue of initial conditions, which has been overlooked up to now. The initial conditions include both the values of x m t and _ x m t at t 0, and the initial history. The initial history refers to the values of x m t and _ x m t before t 0. In a tracking simulation where the history is truncated at N bunches, the initial history has to be specified over the time interval ÿN; 0. In fact, it has been completely ignored so far-the initial history has been set to zero in all of the simulation in the previous sections, as well as in earlier publications [4,6]. The solution of the equation of motion obviously depends on values of x m t and _ x m t over the whole time interval. Since there is no obvious reason why all values should be neglected except those at t 0, it is perhaps surprising that we have good agreement between simulation and analytical results, particularly since we have observed that modes that are expected to decay according to the analytical theory can, after a period of time, start to grow exponentially. We suggest that good agreement is usually observed because, for some reason, x m t is very insensitive to the initial history. This hypothesis will be checked numerically in Sec. V, where we also derive an analytical formula for x m t for the case where the beta function is constant. In the present section, we prepare the necessary groundwork. We do so by characterizing the behavior of the solution for the single mode equation of motion, as a function of the initial conditions at t 0. The initial history will continue to be set to zero in this section.
When the beta function is constant, the mode coupling term is zero, and Eq. (13) simplifies to x t Kx t X n1 b n e i2n=Mx t ÿ n: (20) This equation contains a single mode, fully decoupled from all other modes. This could offer a fast simulation method if we just want to simulate one mode only. To carry out the numerical integration of this equation, we are immediately faced with the problem of initial conditions. It turns out that, if we use random values for the initial conditions, as we have done for the tracking simulation in Sec. II, we may not find the expected behavior. For instance, if we integrate the equation for mode number 100, which has a negative growth rate according to Fig. 3, we would expect a decay mode. However, if we use random initial conditions, the decay mode may in some cases start to grow again, as we shall show in this section. This is not because of mode coupling, which is absent for constant beta function. We do know that the normalized action gives the ''correct'' decay mode behavior. This suggests that we try initial conditions that match with the normalized action. At this point, we must address a technical detail that we have postponed up to now. In the algorithm for the tracking simulation in this paper and in [6], the mode has been obtained by performing a DFT on The mode, however, is defined by Eq. (6), which really corresponds to the inverse DFT. This is the form that has been used in our analysis, as well as in Ref. [5]. A DFT on 2J m p e i m would give results with the same amplitude as an inverse DFT on 2J m p e ÿi m (with a negative imaginary exponent). This has not made any difference to the growth rates simulation because only the amplitudes have been computed. However, to use Eqs. (5) and (6) to analyze the effect of using the normalized action, we should look at 2J m p e ÿi m . We drop the 2J m p initial amplitude, which can be inserted into the solution later if desired, because the equation is linear. When the beta function is constant, m ! t. It turns out that we must choose the initial conditions to be the same as those of the normalized action before we can get a mode that decays all the time. To do so, we can take the initialx to be e ÿi m e ÿi! t . _x would then be ÿi! e ÿi! t . This places a constraint on the initial conditions. Initial conditions that satisfy this constraint are said to be properly ''matched.'' To understand the significance of matching the initial conditions, we shall integrate the equation both for the case when the initial conditions are properly matched, and when they are not.
A set of properly matched initial conditions can be specified as follows: For t 0;x 0 1 and _x 0 ÿi! : (23) With these conditions, we can start the integration. We shall take N 100 as we have done in Sec. II. Figure 8 shows the result for mode number 500 over 600 turns. The inset shows the result over the first 0.1 turn. Over this short interval, as expected, the motion is very nearly simple harmonic. As we have hoped, it gives a nicely decaying mode. Next, we should check if this equation gives the correct growth rate for all of the modes. We compute the initial gradient of the amplitude of everyỹ over the first turn, and the result is the same as Fig. 3. So the single mode equation (20) works for simulation of a single mode. Now we shall show that, if the initial conditions are even slightly mismatched, then a mode that initially decays can eventually start to grow. We introduce some error into the initial condition used above and see what happens when we integrate the equation. For simplicity, we shall add the error to the velocity only. The initial conditions are now x 0 1 and _x 0 ÿi1 "! : The initial history over ÿN; 0 is again set to zero. The integration results for mode 500 when " 0:1 are shown in the Fig. 9. Notice that the decay mode grows again after some time. We have repeated the calculation for " ÿ0:1, ÿ0:01, ÿ0:001, 0:001, 0:01, and 0:1. According to the analytic result, the growth rate for mode 500 is negative, so we should expect to see a decaying mode. In all cases except " 0, the decay mode grows again after some time. This result suggests that modes are only observed to decay without later growth if the initial conditions are properly matched. The question then naturally arises as to whether this is indeed the case, and if so, why?
To begin to answer this question, consider the elementary solution x e ÿi t : (25) Substituting into the single mode equation (20), we get the characteristic equation Equations of this form have multiple roots [13][14][15][16]. To get an idea of the distribution of the roots, we can move all terms to one side and make a contour plot of the absolute value of the resulting function. Figure 10 is a contour plot of jf j 2 when N 100, where FIG. 9. Single mode simulation result for the amplitude of mode 500 when an error is introduced into the normalized action initial conditions. The amplitudes are displayed every 10 turns, for 0.1 turn.
In the contour plot, we have transformed the time variable according to so that the time interval between adjacent slices is now 1. This makes the scale easier to read. From such a snapshot, we can learn a few things about the roots: (i) The root familiar from the analytic growth rate formula is the one that is close to ! . This is labeled (a). (ii) There is a root in an apparently symmetric position, labeled (b). (iii) There exists a series of roots with much larger and negative imaginary parts, and that appear to run along the direction of the real axis. These are labeled (n), where n is an integer.
We have plotted the contour over a much larger area, and the pattern appears to be the same. There could well be other roots, but we shall start with these.
Based on this knowledge of the roots, a plausible form of the general solution is x Ae ÿi a t Be ÿi b t X 1 nÿ1 C n e ÿi n t : (29) In principle, the coefficients A, B, and C n are obtained by fitting this equation to the initial conditions. An infinite number of C n are required because the initial history can be any function over the time interval ÿN t < 0. Thus, the infinite sum on the right-hand side of Eq. (3) is similar to a Fourier sum. We shall not touch on the issue of completeness of the elementary solutions, an area of active research in the field of delay differential equations [19]. Since n has large, negative imaginary parts, they will be damped out rapidly. A quick calculation for mode 500 shows that, for the roots visible in Fig. 10, the time constant is about 0.01 turn. However, we do not know if they can be completely ignored because when we repeated the calculation for N 1000 and 10 000, we find that n move closer to the real axis, so that the damping time increases. In the next section, we shall use numerical means to determine the errors that result from ignoring these roots. For now, we shall assume that the dominant terms on the right-hand side of Eq. (30) are just the first two: x Ae ÿi a t Be ÿi b t : Roots (a) and (b) are expected if the wakefield is small. They are close to ! , which are the roots of the simple harmonic equation obtained when the wakefield sum on the right-hand side of Eq. (1) is set to zero. In particular, root (a) is the one that is used in the derivation of the analytical formula for growth rates given in Appendix A. For mode 500, we know that the first term is decaying. When we magnify the contour plot above around root (b), we can see that it has a positive imaginary part. So the second term is growing. Drawing on the similarity with the simple harmonic equation, we can see that both terms should be important. Yet, only the first term has been considered so far in this paper, as well as in Refs. [4 -6]. We shall show that it is the second term that caused the mode amplitude to grow when mismatched initial conditions were used in the integration of Eq. (20). In order to show this, we shall derive the analytic solution with mismatched initial conditions. First, note that root (a) is very close to ! , and root (b) is very close to ÿ! . For small t, the modes evolve as given by Eq. (30): x t Ae ÿi! t Be i! t : We consider mismatched initial conditionsỹ 0 1 and _ỹ 0 ÿi1 "! . Solving for A and B at t 0, we get For large t, this evolves into It is possible to obtain accurate values for roots (a) and (b) from the contour plot using a minimum search algorithm. However, the known analytic formula for root (a)-in Appendix A-is sufficiently accurate for small wakefields as we shall show, and a corresponding formula will be derived for root (b) later in this section.
Using this analytic solution, it is easy to plot the modes for the ''matching error'' used above. Figure 11 compares the simulation result for " 0:1, with the analytical result from Eq. (33). The agreement is excellent, with an error of about 0.3% at 490 turns. This means that we now understand why we cannot choose random initial conditions if we want to see the decay mode for mode number 500: (i) The matched initial conditions give the expected results because in that case, the solution is overwhelmingly dominated by root (a), the same one that is used in the analytic formula for growth rates. (ii) There is another root, (b) that is equally important, if we want to determine the correct trajectory for each bunch.
We now have an explanation as to why using matched initial conditions gives the expected damping or growth, and why more general initial conditions (including, for example, an initial history) leads to growth in modes that are expected to be damped. However, this does not mean that the earlier results on the growth rates are wrong. It just means that, for each mode, we have extracted only the contribution from root (a) from the tracking simulation. This has taken place through the use of the normalized action in Sec. II. The actual results of the tracking simulation must necessarily contain contributions from root (b) also, which we have not considered up to now.
We shall complete this section by deriving the analytical formula for root (b). The steps are nearly the same as those that have led to root (a). To make contact eventually with the growth rate formula in Appendix A, we shall make the resistive wakefield more explicit. From the wakefield formula also in Appendix A, the wake coefficients in the equation of motion are proportional to 1 z p . Substituting the elementary solution e ÿit into the single mode equation of motion, Eq. (5), we obtain the characteristic equation Factorizing, For small wakefields, ! , so This gives the formula for a [5]: However, there is also another root when ÿ! . Starting again from when ÿ! , Then b is given by The symmetrical nature of Eqs. (38) and (41) means that graphically, the imaginary part of b is a simple reflection of the imaginary part of a , which is the growth rate. These are plotted in Fig. 12(a) for N 100. The real part corresponds to the frequency, and is plotted in Fig. 12 An interesting point to note is that the imaginary parts of a and b are both negative for the following mode numbers:

< and M ÿ < < M:
These modes will always decay, even if there is an ''error,'' i.e., even if the initial conditions differ from the normalized action, and will therefore not produce the growth seen in mode 500 above.   We are now ready to explain the small oscillations in the growth rates in Fig. 3 obtained from the tracking simulation. Figure 13 shows the growth rates obtained from the above analytical formula for a in Eq. (38), truncated at N 100 (labeled ''Analytic N 100''). This clearly shows the same oscillation as the simulated growth rates, also shown on the same figure labeled ''Simulation N 100.'' On the same graph, the analytic results for N 500 are plotted. The amplitudes gradually decrease as N increases. In the limit that N goes to infinity, it should converge to the analytic result given by the Fourier transformed expression in Appendix A, labeled ''Analytic N infinity'' in the figure. We have also repeated the calculations of growth rates by doing a minimum search for the root (a) on the contour plot in Fig. 10. These results, when plotted on the same graph in Fig. 13, are visually indistinguishable from those already obtained from Eq. (38). The spiraling curves displayed in Fig. 14

V. ANALYTIC SOLUTION
In this section, we derive an analytic solution for the equation of motion for the case of constant beta function. An analytic solution has been derived in Ref. [4]. However, this only makes use of the solution closest to ! . It does not include the root b , nor is there an explanation of why the rest of the roots n should be neglected. The use of only root a means that the analytic solution will match the actual solution only when the initial conditions are properly matched, as we have demonstrated in the previous section, and not for general initial conditions. Here, we shall derive a solution that is valid for any initial condition, and measure its accuracy by comparing it with the case when the history is neglected, i.e. set to zero. The importance of initial history in delay differential equations has been widely studied (see, e.g., Ref. [16]). In the cases of wakefield coupled bunches that we have considered so far, it appears that correct growth rates may be obtained even without the initial history being considered. However, the initial history may become important if the system parameters are changed significantly, and it is important to know for what parameters this may happen.
We approach the problem as follows. We know from the example in Fig. 11 and fitting the coefficients A and B to the initial conditions at t 0 only. If this solution is sufficiently accurate, we can apply it to all modes and derive an analytic formula for the trajectory of every bunch. For such a formula to be useful, we must first obtain an estimate of the actual errors involved. We know that there is an error of 0.3% in Fig. 11. There are two possible sources of this error: one is the kick method used to perform the numerical integration, and the other is the neglect of the initial history. The reason why we expect the latter error is as follows. In the simulation, the initial history has been set to zero. However, if the true trajectory is in fact given by Eq. (42), the initial history over the time interval ÿN; 0, must also be given by Eq. (42). The reason is that, if the initial history is not given by Eq. (42), then it must be given by Eq. (29) which is reproduced here: x Ae ÿi a t Be ÿi b t X 1 nÿ1 C n e ÿi n t ; (43) where the coefficients C n are nonzero. However, if C n are nonzero over the time interval ÿN; 0, they must continue to be nonzero for t > 0, since these coefficients are constant. Then Eq. (42) would not give the actual trajectory. Therefore if Eq. (42) is the actual trajectory for t > 0, then the initial history must also be given by Eq. (42). This reasoning gives us a way to check the numerical error due to the kick method. All we need to do is to calculate the initial history using Eq. (42) for some values of the coefficients A and B, carry out the simulation, and compare the results with Eq. (42). As an example, we set A 1:05 and B ÿ0:05, which would give the same initial conditions at t 0 as in Fig. 11. We compute the initial history at t ÿN; ÿN ÿ 1; . . . ; ÿ2; ÿ; 0. This and the subsequent simulation are illustrated in Fig. 15 for mode 500. The simulation result up to turn number 500 looks the same as in Fig. 9. If the simulation method is perfectly accurate, it should give exactly the same results as Eq. (42); any difference must represent the error due to the method.
The relative error from the kick method for mode 500 is shown in Fig. 16. It appears that the simulation method gives a systematic error that increases with time. In order to determine the error due to the neglect of the initial history, we must separate out this systematic error. We shall take the following approach. Ideally, if we can compute the trajectories from all possible initial histories and compare them, we can determine how strongly they depend on initial history. This is of course not possible because the initial history is infinite dimensional, i.e., an initial history can be any function on the time interval ÿN; 0. The next best option then may be to take a small random sample. One way is to construct the initial history as a Fourier sum of random coefficients, and carry out the simulation. To obtain the error due to the initial history, we should not compare the result with Eq. (42), because this would include the systematic error due to the kick method. Instead, we can compare it with the simulation result obtained for an initial history given by Eq. (42), which we used in the previous paragraph. In this way, the systematic error can be removed. The error due to the initial history, for one set of random Fourier coefficients, is shown in Fig. 17 for mode 500. This error is much smaller than the systematic numerical error in Fig. 16. It is interesting to note that this history-induced error approaches a constant at large turn numbers. This suggests that, in this case at least, the error induced by neglecting the initial history does not grow, which is reassuring. We shall use this history-induced error to provide a measure for the accuracy for the analytic formula.
We repeated the calculation for a number of random Fourier coefficients, as well as for zero initial history, and obtained similar results to those shown in Fig. 17. Using the zero initial history, we compute the history-induced error for every mode. The result is plotted in Fig. 18. The error peaks at mode numbers close to the betatron tune. Notice that the peak error is just over 1%, which is not as small as we would like it to be. Fortunately, the peak is narrow, and for most of the other modes, the errors are of order 0.001%. Together with the property that the error remains constant at large turn numbers, this result could provide enough confidence for us to construct an analytic solution.
Before writing down the analytic solution, we should look at the parameters that affect this history-induced error. First, combining Eq. (5) and the formula for the resistive wake function in Appendix A, we write the equation of motion in the following form: All the issues with the infinite number of roots and the initial history arise from the presence of the right-hand term. Therefore, the effect of the initial history depends  is large, the effect of the initial history will be large. To quantify the effect of b 0 , we insert a factor, w, on the right-hand side of the equation: x The parameter w scales the strength of the wakefield; we shall call w the ''wakefield strength.'' We compute the history-induced error for mode 500 for different wakefield strengths. The result is shown in Fig. 19. As expected, the error increases with wakefield strength. The log-log plot follows a fairly straight line, at least for this particular mode, which is convenient for error estimation. Recall our estimate above of the 4% error for the analytical solution of a bunch trajectory. This figure shows that, for the ILC damping ring, it may just be possible to construct an analytic solution with acceptable accuracy based only on a and b . The accuracy would decrease for larger wakefields.
We are now ready to write down the analytic solution. Based on the above results and discussion, we first make the assumption that the wakefield must be sufficiently small. The effect of the initial history may then be neglected, and only two elementary solutions are significant for each mode. For a very small wakefield, the dominant solutions are obtained from the following roots of the characteristic equation. Instead of using Eqs. (38) and (40) where we have to truncate at some N, we use the Fourier transformed version that can be obtained from Appendix A: which converge much more rapidly. Equation (47) can be obtained from Eq. (46) by comparing the sign changes between Eqs. (38) and (41). These give us the accurate limit when N goes to infinity, as required by the original Eq. (1). Next, we compute the mode amplitudes of the initial conditions x m 0 and _ x m 0: The solution for each mode is then The unknown constants are obtained with the help of the initial conditions. At t 0, we have and Solving, we get and Finally, we can invert the discrete Fourier transform using and obtain the analytic solution for every bunch: 1.E-11 1.E-08 1.E-05 1.E-02 where the imaginary part is to be neglected since x m t is real.
We have arrived at this solution heuristically, and a test on its validity should be carried out. As before, we do a comparison between tracking simulation and analytic solution. For convenience, we combine the errors in both displacement and velocity into one measure, by using the distance between points in a phase space. Let displacement and velocity be the two coordinates on a Cartesian plane. To give equal emphasis to both quantities, we scale the velocity by a factor of 1=! . The reason is because, for small wakefield, if the displacement is given approximately by e i! t , then the velocity is given by the same expression multiplied by i! . At a particular time t, let the point x; _ x =! be given by the actual values of the displacement and velocity, and let x e ; _ x e =! be given by the estimated values. We then define the distance between these two points as the error. We also define the relative error as the error divided by the distance from the origin to x; _ x =! . The relative errors for bunch number 100, shown in Fig. 20, are computed in a similar way as those for the single mode case in the previous section. The roots used in the analytic solution are obtained using Eqs. (38) and (41) instead of Eqs. (46) and (47), since we are comparing with simulation for N 100. The numerical error is obtained by comparing the analytic solution with simulation results using a history created from the analytic solution, starting with a set of random x m and _ x m at t 0. The history-induced error is obtained by comparing the simulation results for the case where the history is created from the analytic solution, with the case where the history is set to zero. Thus, the accuracy of the analytic solution up to 100 turns for the bunch trajectory is about 0.07%. If we were to plot the analytic and simulation results for x m and _ x m =! on the same graph at about 100 turns, they would look very similar to Fig. 11.

VI. CONCLUSION
We have demonstrated that the Fourier modes are in fact not the eigenmodes in a real lattice as a result of beta function variation. The concepts of modes and growth rates remain approximately valid in a real lattice in the short term. For mode 100 in the ILC damping ring (Fig. 4), short term would mean less than 30 turns. This is of the order of the time required for injection and extraction [20], as well as the damping time of the feedback system. Over the long term, the strong oscillations in the beta function give rise to coupling between modes. The mixing of growth modes into decay modes can cause the decay modes to grow. This observation casts some doubt on the accuracy of the growth rate formula in Appendix A. In Appendix E, for instance, we can see on Fig. 24(a) that, even for a very simple lattice, the actual growth rate is larger than that predicted by the growth rate formula. However, if we are mainly interested in the maximum growth rate, the formula may still provide a good estimate. The periodic variation of the beta function, with period equal to the time of revolution around the damping ring, also gives rise to additional oscillation. The frequency of this oscillation is equal to the number of wavelengths of a mode passing through a given point per unit time.
We have also shown that the growth of decay modes is not restricted to large, complex rings like those of the ILC, but could also take place in small, simple rings consisting of no more than repeating FODO cells and a small number of bunches. A FODO lattice consists of an alternating sequence of focusing (F) and a defocusing (D) quadrupoles, separated by drift spaces (O) of equal length. In general, we may imagine that the eigenmode can be changed if any assumption in the macroparticle model is changed. Practical examples include nonuniform fill patterns and localized wakefield. In the case of the ILC damping ring, we have calculated the corresponding effects due to nonuniform fill patterns in Appendix C and shown that the effects of the beta function variation is comparable. In Appendix D, we have also computed the effects of localized wakefield due to the transverse HOMs in the kind of rf cavities expected for the damping rings, and shown that these effects are quite small in comparison to the beta function variation.
Although the Fourier mode is only an eigenmode in the case of the constant beta function, it remains useful because this case is amenable to analytic solution. It is therefore important that the solution be properly characterized. As previously published solutions [4] are not applicable to general initial conditions, we have derived more complete solutions to address this gap. We have also shown that the validity of such solutions depends directly on the strength of the wakefield. The reason is that the wakefield gives rise to a dependence on initial history of the bunch trajectories, which is not likely to be available in practice. The analytic solutions would therefore be useful only for small wakefields. The ILC damping ring parameters we use corresponds to a current of about 580 mA. (This is rather higher than the average because the damping rings tend to have lower bunch number or particles per bunch, whereas we have taken the worst case in our simulation.) Increasing the current beyond this value would quickly increase the error in the solutions, as Fig. 19 shows.
Simulation using the kick method is required when the beta function oscillates strongly, but it is computationally intensive and the numerical error grows with turn number. The analytic solution is fast but is useful at large turn number only if the beta function is constant. Otherwise, it cannot account for the growth of decay modes resulting from mode mixing. Nevertheless, it is important as a theoretical reference for the simulation. For machines requiring very high levels of beam stability, particularly in the presence of strong transients (as in the case of the ILC damping rings), a thorough and rigorous understanding of the dynamical effects associated with long-range wakefields will be essential. Together, simulation and analytic solution can provide a more complete picture of the bunch trajectories, and will be useful tools for determining the specifications of the feedback systems and estimating the impact of transient effects.

ACKNOWLEDGMENTS
We are grateful to the Science and Technology Facilities Council for funding through the Cockcroft Institute.

APPENDIX A: STANDARD FORMALISM FOR COUPLED-BUNCH INSTABILITIES
We summarize here the standard formalism for coupledbunch instabilities relevant to the work in this paper. In the macroparticle model for bunches in a storage ring, a number of assumptions are often made to simplify the equations of motion: (i) Each bunch of charged particles (in our case, electrons or positrons) is treated as a single, rigid particle. Only the motion of the centroid is considered. The distortion of each bunch is neglected. (ii) Bunches travel at the same speed around the ring, and distances between adjacent bunches are equal. In our case, we also assume that this speed is very close to the speed of light. (iii) The displacements of bunches from the design trajectory is small, so that a linear approximation of the wake force is valid. This wake force is proportional to the displacement of the exciting bunch.
The resulting equation of motion for each bunch takes the following form [5]: The notations are as follows: x m t Transverse displacement of the mth bunch ! Betatron frequency W 1 z Wake function Time taken by light to travel from one bunch to the next N Number of particles in each bunch r 0 Classical electron radius c Speed of light T 0 Time taken for each bunch to travel one round around the ring Energy of each particle in units of its rest mass The argument, z, of the wake function is the distance between the exciting bunch that produces the wakefield, and the spectator bunch that experiences it. ! is assumed to be a constant here, but variation of the beta function can have significant effects, as we see in Secs. II and III. When the betatron frequency is assumed to be constant, we use the averaged value given by the following expression: where is the transverse (vertical or horizontal) tune. An important contribution to the long-range wakefield in many storage rings comes from the resistive wall of the beam pipe. For a beam pipe with conductivity and circular cross section of radius b, the wake function is given by [5] Equation (A1) may be understood qualitatively as follows. The left-hand side has the same form as the equation for a simple harmonic oscillator. This represents the effect of the focusing forces arising from the focusing elements. The right-hand side is the sum of the wake forces arising from all of the bunches in front of bunch m. This sum is to be taken to infinity. Note that if we start at bunch m and go around the ring, we would come back to bunch m itself, so that each bunch can see its own wakefield. In the absence of any wakefield, the right-hand side is zero and bunch m undergoes simple harmonic motion at the betatron frequency, ! .
It should be mentioned that we have left out damping effects, such as those coming from synchrotron radiation, decoherence, and action of a feedback system. We are mainly interested in the qualitative behavior produced by the interaction between the focusing term ! 2 x m t, and the wakefield coupling on the right-hand side of Eq. (A1). A damping term, 2& _ x m t, can be inserted on the left-hand side when required. The effect will be to subtract & from the growth rate formula given below. The modification to the analytic results is straightforward and has been given in [4]. Next, it is often more convenient to work with modes instead of the bunch displacements directly. A mode is essentially a discrete Fourier transform (DFT) of the bunch displacements. Suppose there is a total of M bunches. The mode number is defined by This has been described as a normal mode [4], a multibunch mode [5] or a Fourier mode [6]. The advantage of using modes is that the equations of motion, which are fully coupled in Eq. (A1), become completely decoupled, as in Eq. (5). However, decoupling only really occurs in the case of constant beta function. We show in Sec. III that if the beta function varies around the ring, the modes remain coupled.
For constant beta function, it is possible to derive an elementary solution for each mode in the following form: x t e ÿi t : (A5) By substituting this into the transformed equation of motion, and solving the characteristic equation under certain approximations (see Sec. IV), the following formula is obtained [5]: where Z ? 1 ! is the Fourier transform of the wake function in Eq. (A3). The growth rate is the defined as the imaginary part of : Im : (A7)

APPENDIX B: REDERIVING THE KICK METHOD
The equation of motion given by Eq. (A1) is a delay differential equation. If we truncate the sum on the righthand side at the Nth term, it can be written in the following form: Suppose we want to integrate to obtain x m t from t 0 to . In order to do so, we need to know the history of each bunch represented on the right-hand side. Specifically, we need to know x m1 t from t ÿ to 0; x m2 t from t ÿ2 to ÿ ; . . .
x mN t from t ÿN to ÿ N ÿ 1: The history of x mn t is ''stored'' at the time slices, i.e., its values are known at t 0; ÿ; ÿ2; . . . ; ÿN. If the time interval is short enough, we can make a linear interpolation at each time interval to obtain the value of y mn t in between the slices. Thus, for t ÿn to ÿn ÿ 1, we may write the interpolated function as x mn t a n t n c n : Since the values of x mn t are known at t ÿn and ÿn ÿ 1, we have x mn ÿn c n ; x mn ÿn ÿ 1 a n c n : The constants a n and b n can then be calculated. The nth term on the right-hand side of Eq. (B2) is approximated by a straight line: x mn t ÿ n a n t c n : The equation of motion can be written in the form b n a n t c n : Defining p X N n1 b n a n and q X N n1 b n c n ; (B6) the equation of motion becomes The solution can be obtained analytically and is given by Let the initial conditions be Since q is in fact the value of wake force on the righthand side of Eq. (B1) when t 0, q is exactly the same as the kick that is used in the kick method described in Sec. II. Therefore, we have demonstrated that the kick method is equivalent to integrating the equation of motion with a linear interpolation on the history of the bunches.
In general, ! 2 in Eq. (B1), which represents the focusing strength of the magnets, can also be 0 or negative [21]. These correspond, respectively, to drift space or defocusing region. For these cases, the above steps can be repeated in a similar way and would lead to the same first-order result.

APPENDIX C: MATRIX METHOD
We express problem of the varying beta function in matrix formalism, and compare this with the case of nonuniform charge distribution. The objective is to highlight the difference between these two cases, as well as the difference from standard linear system analysis that arises because of the time delay involved in the coupled-bunch problem.
In the varying beta case, the beta function is fixed with respect to the lattice as the bunches move through the lattice. In the nonuniform distribution case [22], the distribution moves with the bunches. As a result, although the Fourier modes no longer form the eigenmodes for the nonuniform distribution case, it is nevertheless possible the find a new set of eigenmodes by diagonalizing the matrix that describes the dynamics with the wakefields. In the varying beta case, although this matrix can also be diagonalized at any instant in time, the matrix elements change with time because of the relative movement between the bunches and the beta function. This means that a different eigenmode may be required for each time step, which is not very useful.
In either case, because the equations of motion resemble closely those of a system of linear differential equations [23], the direct approach is to express it as a matrix eigenvalue problem. Solving for the eigenvalues and eigenvectors using standard methods then gives the complete solution to the problem. This is essentially the approach that has been taken for both the uniform distribution case [5] and the nonuniform distribution case [22]. In the uniform distribution case, the eigenvectors are the Fourier modes, and the eigenvalue for each mode can be obtained by making the approximation that ! , as is done in Eq. (37). However, if the system depends on the history of the dynamical variables, it is possible to have more than one eigenvalue for each Fourier mode, as shown in Eq. (29). In the matrix formalism, as we now show, this is because the matrix itself is also a function of the eigenvalues.
We start with the case of uniform distribution and constant beta function. The equation of motion is given by Eq.
Notice that the matrix elements in each row are exactly the same as those in the previous row, but shifted one place to the right. This type of matrix is called a circulant matrix, and it is known that the eigenvectors of such a matrix are the Fourier coefficients, e i2n=M . The matrix equation above forms the linear system that we have mentioned.
In standard systems, we can solve the eigenvalue problem and obtain the complete solution. In this case, however, we notice one complication -the matrix C is a function of , which also appears in the eigenvalue ÿ 2 ! 2 . If we trace this back to the original equation of motion, we see that this arises from the term x mn t ÿ n, where the value of the dynamical variable at an earlier time is required. Since we do not know , we do not know the matrix elements in C, and we cannot solve the problem directly. As explained earlier, the approach taken in Ref. [5] at this point is to approximate the in C by the known quantity ! , so that we have The matrix C is now specified, and we can proceed to solve the eigenvalue problem. The eigenvectors then give Fourier modes that behave exponentially with time.
In this paper, we have been motivated by the observation on exponential behavior to solve the problem exactly without making the above approximation. Although we do not know the value of the elements in C, we do know that this matrix is circulant. We can write down the solution to the problem, in terms of the unknown , using the results for circulant matrices from Ref. [24]. This gives the eigenvectors x 0 1; e i2=M ; e i22=M ; . . . ; e i2Mÿ1=M T (C8) for 0; 1; . . . ; M ÿ 1, and the eigenvalues On substituting Eq. (C6), this can be shown to be equivalent to the characteristic equation (26), which has an infinite number of solutions. Thus, unlike the standard linear system where each eigenvector corresponds to one eigenvalue, we have a time delay system in which each eigenvector has an infinite number of eigenvalues [16].
We define a few matrices here that will be useful in subsequent discussion. Let where U is the Hermitian conjugate of U. Next, we consider the case of nonuniform distribution pattern, but still constant beta function. The equation of motion may then be written as where f m is the number of particles in bunch m, relative to the number of particles in each bunch in the original uniformly distributed bunches. Following through the above reasoning, the matrix equation is then given by with the off-diagonal coupling terms given by ÿk ÿ 0 =M, as has been discussed in Sec. III. Here, we have used D n to denote the diagonal elements in D.
Comparing Eqs. (C16) and (C17), it is clear that the coupling terms in the case of varying beta function depend on time, and the diagonalization could in general be different from one time step to the next. In principle, we can compare the off-diagonal elements in the two matrices to see which of the two effects-varying beta or nonuniform distribution-is stronger. While Eq. (C16) is fairly straightforward to compute, Eq. (C17) would be difficult given the complexity of the lattice. Instead, we can simulate the time evolution of each and every mode in the same way as we have done for Fig. 4. We do this for both the varying beta and the nonuniform distribution. For the nonuniform distribution, we simulate approximately one of the cases that are considered for the ILC damping ring [20]. We consider minibunch trains of 22 bunches each, separated by gaps that correspond to the distance of 20 bunches. Thus we have 22 consecutive bunches, followed by 20 empty spaces, and this is repeated all around the ring. The last part of the ring, into which this repeating unit of minitrain plus gap cannot fit, is left vacant. We selected all of the modes that would be expected to decay exponentially in the uniform distribution, constant beta case, as predicted by the growth rates in Fig. 3. These are all plotted on the same graph, and the result is shown in Fig. 21. The graph shows that many of the modes that are expected to decay will eventually grow. Figure 21(a) shows the decay modes for the case of the varying beta function, and Fig. 21(b) those for the case of the nonuniform distribution. The numbers in the legend are the mode numbers obtained from the modes in Fig. 3 with negative growth rates. For illustration, the graph for mode 100, in particular, has been plotted in Fig. 4. Figure 4 also contains the graph for the case of constant beta and uniform distribution, where the amplitude decays monotonically -also observed for all other decay modes in the simulation. It is clear from Fig. 21 that the effect of varying beta function is comparable with the effect of nonuniform charge distribution, if not stronger.

APPENDIX D: HIGHER-ORDER MODES
We estimate here the effect of a main source of localized wakefield, that of the transverse HOM in the rf cavity of the ring, by carrying out the tracking simulation to determine the effect on the Fourier modes. In order to do so, we first calculate the wake function of a typical rf cavity recommended for the ILC damping ring [10], and obtain a formula for the transverse kick on each bunch by the HOM.
The cavity parameters are obtained from the KEK-B superconducting rf cavity [25], listed in Table I. We assume that there are two sets of nine cavities each, arranged in diametrically opposite positions of the OCS6 damping ring, and that the cavities in each set are fairly close, of the order of 1 m apart from each other. Each set of cavities is approximated as a thin lens, and gives a single kick to any bunch passing through the set. The wake function is obtained by Fourier transforming the impedance, which is given by [10]  where N cav 9. Note that this is in SI units. The wake function is then given by [5] Wz ÿ i 2 for z < 0, and zero for z > 0. In order to relate this to the earlier formulas, it is first converted into cgs units using Then, following the procedure in [5], we obtain an expression for the kick (change in velocity) on the mth bunch when it goes through the cavity: The factor t=T 0 arises because W 1 is the wake function for the whole circumference, whereas the kick is only for one time step. At this point, it is useful to compare the relative magnitudes of the two wake functions. The appropriate quantities to compare are W HOM and the resistive wall wake function, W 1 . These are plotted together in Fig. 22. The HOM wakefield is significantly weaker than the resistive wall wakefield. Note that the HOM wakefield is often expressed as W HOM divided by the length of the cavity, so that the unit is V=pC=mm=m. We have not expressed the wakefield in this form because we need to compare the resistive wall wake-  field for the whole circumference with the HOM wakefield for one set of cavities, which have different lengths.
From the above comparison, we would expect that the effect of the HOM wakefield is small. We have carried out the simulation to compare the effect of the HOM wakefield on the Fourier modes and found that this is indeed the case. The HOM kick is applied once to each bunch when it reaches the middle of each set of cavities. Each kick is a superposition of the wakefield of all the bunches that have gone through the cavity earlier. The kick on the mth bunch is given by _ x m ÿ Nr 0 c fW HOM ÿcx m1 t ÿ W HOM ÿ2cx m2 t ÿ 2 g: The simulation is carried out with constant beta function. Over the first 100 turns, the plotted results of the mode amplitudes with and without the HOM wakefield are visually indistinguishable. The result for mode 100 is shown in Fig. 4.

APPENDIX E: SIMPLE LATTICE
We demonstrate here that the phenomenon of growth of decay modes can occur even for a very simple lattice. We use a lattice consisting only of ten FODO cells. Filled with four evenly spaced bunches. The beta function is plotted in Fig. 23.
The machines parameters are summarized as follows: In order to increase the wakefield coupling between the bunches so that growth and decay of the modes may be observed more quickly, we have used a smaller radius and a lower conductivity for the beam pipe. We assume that the beam pipe is made from the type of stainless steel with the lowest conductivity listed in [26]. The simulation follows the same procedure as the kick method described in Sec. II, and uses the resistive wall wake function in Eq. (A3). We have checked for convergence with respect to the number of time slices in the ring, and the number of terms to keep in the wakefield sum in Eq. (1). We plotted out the displacement x 0 t at the 1000th turn and checked visually that it has converged for 100 time slices and 400 terms in the wakefield sum. In the case of four bunches and a constant beta function, with ! given by Eq. (A2), there are four modes. According to Eq. (A7), modes 1 and 2 are growth modes, and modes 0 and 3 are decay modes. The results for modes 2 and 3 are shown in Fig. 24. On the same graphs are plotted the analytic corresponding results when the beta function is a constant. Notice that, for the actual beta function, the mode 3 amplitude follows the analytic result initially. It decays for the first 100 turns, and then grows after that. Mode 0, not shown, behaves in a similar way to mode 3. Modes 1 and 2, the growth modes, both tend to have slightly higher growth rates than the analytic results over the 1000 turns. This behavior suggests that the variation in beta function in any lattice could lead to decay modes growing, and that this feature is not restricted to large, complex rings like the ILC damping ring. This is a result of the breaking of translational symmetry along the lattice; since bunches no longer experience the same focusing strength at different times, the Fourier modes are no longer eigenmodes. As we have see from Eqs. (8)- (13), the equations of motion are not decoupled by transforming to Fourier modes.
Finally, we should mention the interesting case of using two bunches instead of four. In the case of two bunches and 10 FODO cells, our simulation shows that the decay mode amplitude does not grow. The reason is as follows. Since 10 is an even number, one bunch is exactly 5 FODO cells away from the other bunch, so that the beta function experienced by the two bunches are identical. In Eq. (11), this means that kt m are equal for all m even though they can all vary with time. Equation (11) can therefore be rewritten as x t Kx t X n1 b n e i2n=Mx t ÿ n ÿ ktx t; (E1) so that the equations of motion can become fully decoupled. Therefore the Fourier modes are still the normal modes for this special case. This does not necessarily mean that the mode amplitude must be exponential in time, although in this case it is. Our simulation shows the mode amplitudes follow the analytic results fairly closely, and the decay mode does not grow again-at least not within the first 1000 turns.