Stability of longitudinal bunch length feedback for heavy-ion synchrotrons

In heavy-ion synchrotrons such as the SIS18 at Helmholtzzentrum für Schwerionenforschung, Helmholtz Centre for Heavy Ion Research (GSI), coherent oscillations of the particle bunches are damped by rf feedback systems to increase the stability and to improve the beam quality. In the longitudinal direction, important modes are the coherent longitudinal dipole and quadrupole oscillation. In this paper we present a new and rigorous approach to analyze the longitudinal feedback to damp these modes. The results are applied to the rf feedback loop at GSI that damps the quadrupole mode. The stability analysis is compared with simulations and is in good agreement with results of a beam experiment. Finally, we summarize practical implications for the operation of the feedback system regarding performance and stability.


I. INTRODUCTION
To obtain high quality beams that are provided reliably for experiments, many efforts are made to refine and optimize accelerator components to prevent beam oscillations or instabilities. Feedback systems are an indispensable tool to stabilize the beam and to damp oscillations that occur steadily during operation. One large new accelerator complex that is currently built is the Facility for Antiproton and Ion Research, called FAIR. This center will expand the facilities of the GSI Helmholtzzentrum für Schwerionenforschung GmbH. The existing heavy-ion synchrotron SIS18 at GSI will be used as a preaccelerator for a new synchrotron ring with a circumference of about 1084 m [1]. For FAIR, new digital feedback systems based on field programmable gate arrays (FPGA) and digital signal processors (DSP) are planned to stabilize the beam in the longitudinal direction. The main goal of the feedback is to prevent oscillations of the shape of these bunches, as these oscillations lead to a degradation of the beam quality [2]. Beam oscillations are usually classified in different modes. This standard theory has been developed in [3][4][5].
This paper deals with longitudinal oscillations of the bunch length. This type of oscillation is usually classified as the quadrupole mode with mode numbers m ¼ 2 and n ¼ 0. Early investigations of feedback for this mode are found in Refs. [6,7]. Since then, considerable progress of digital hardware such as FPGAs has been made and most of the newly developed feedback systems rely on digital hardware. Recent works for the damping of the quadrupole mode are Refs. [8][9][10] for electron machines, Ref. [11] for a proton synchrotron, Ref. [12] for a heavy-ion collider, and Ref. [13] at GSI for a heavy-ion synchrotron.
The new contribution of this paper is a rigorous approach for the modeling and analysis of the feedback loop, taking into account also nonlinearities due to the accelerating rf voltage. The structure of the paper is as follows. Section II describes the single-particle longitudinal dynamics in heavy-ion synchrotrons and basic notations. The bunch shape oscillations are defined and the control problem is formulated. In Sec. III, a feedback model for the bunch length oscillations is derived. The modeling approach is based on moments and is applicable for quite general nonlinear accelerating voltages. Because the moments are not readily available for measurements, it is shown how they can be obtained from the measured beam current signal. The presented results are used to derive the closed-loop dynamics of a specific setup for the heavy-ion synchrotron SIS18 at GSI. The stability and performance of the feedback is analyzed in Sec. IV. The analytical calculations are compared with nonlinear macro particle tracking simulations (cf. [14][15][16]) and experimental measurements. A tuning rule for the feedback is presented in Sec. V, the practical implications of the results are discussed in Sec. VI, and a conclusion is given in Sec. VII.

II. NOTATION AND BEAM DYNAMICS
In the following, a single-harmonic rf voltage is assumed with the (angular) frequency ! RF . The reference revolution frequency is denoted by ! R ¼ 2=T R and the harmonic number is h ¼ ! RF =! R . Every particle in the bunch is characterized by its time delay at the cavity and its energy deviation ÁW, both with respect to the reference frame R. For the sake of clarity, we will focus on the stationary case, i.e., the reference voltage U R is assumed to be zero and the reference period T R is constant. The longitudinal equations of motion can be written as [16] where Á k and ÁW k are the time delay and the energy deviation of particle or ion k with respect to R, Q > 0 is the charge per ion, T R is the reference period for one turn, U R is the reference voltage, R and R are relativistic Lorentz factors, W R is the reference energy, and tr is the transition gamma.
The rf voltage at the cavity is chosen as a sinusoidal function whereÛ 0 denotes the target rf amplitude, ! RF is the rf frequency, R is the reference arrival time at the cavity, and the amplitude and phase modulations u 1 and u 2 can be used as inputs to control the beam. In this paper, the stationary case with U R ¼ 0 V, R ¼ 0 s, and R < tr is considered. Figure 1 shows the trajectories in the phase plane ðÁ; ÁWÞ for this case. The origin of this plane is the reference R. It has to be noted that the calculations may be extended to the more general acceleration case below or above the transition gamma tr and also to Q < 0. For convenience, the coordinates Á k and ÁW k are replaced by This leads, with the assumptions of the stationary case, to the single-particle dynamics For small oscillation amplitudes Á k , the dynamics may be linearized and this leads for R < tr to a harmonic oscillation with the synchrotron frequency At GSI, this frequency is typically in the range of a few hundred Hz, whereas the revolution frequency ! R is larger by 2-3 orders of magnitude. For larger amplitudes, the oscillation is nonlinear and for the stationary case, the synchrotron frequency is given by where KðkÞ denotes the complete elliptic integral of the first kind, andx k is the maximum oscillation amplitude on axis x of particle k. Qualitatively speaking, the frequency decreases monotonically and tends to zero forx k ! . It is necessary to distinguish between the motion of single particles with individual synchrotron frequencies ! syn ðx k Þ and the coherent motion of the whole bunch with corresponding coherent frequencies. In a real beam, each bunch may consist of 10 9 -10 12 , or even more, ions. Figure 1 shows a bunch with a matched shape, i.e. its shape is compatible with the trajectories. The individual particles in the bunch will perform synchrotron oscillations with frequencies ! syn ðx k Þ, but the bunch shape itself will be stationary.
A different situation is shown in Fig. 2 by means of a macroparticle tracking simulation. In this figure, the particle density is visualized with a color scale, where red corresponds to higher densities and blue to lower densities. An ensemble of about 2:5 Â 10 5 macroparticles is simulated in the phase plane. At time t 1 , the bunch is assumed to be mismatched. In addition to the phase plane, Fig. 2 also shows the corresponding beam current signal I b ðÁÞ, which is the projection of the phase space density on the axis Á. Because of repetitive circulation of the beam in the ring and the slow variation of the bunch shape, the beam current signal may be regarded as being periodic and can be written as a Fourier series with dc I 0 , higher harmonic amplitudes I k , and phases ' k . The mismatch in Fig. 2 leads to an oscillation of the bunch length, also called quadrupole or breathing mode. At time t 2 ¼ T syn;0 =4, the bunch shape has rotated by about 90 in the phase plane and the signal I b has become higher and more narrow. Obviously, the bunch shape will regain its original shape approximately after half a synchrotron period, thus the coherent oscillation frequency ! bl of the bunch length is roughly 2! syn;0 ; cf. [2]. For increasing bunch sizes, an increasing amount of particles in the bunch will oscillate with a synchrotron frequency ! syn < ! syn;0 and it can be expected that therefore ! bl < 2! syn;0 holds. A further effect of the nonlinearity of the rf voltage is the filamentation of the bunch. This filamentation is already visible at time t 2 and will finally lead to a new stationary distribution as shown in Fig. 2 for a large time t 3 . Filamentation leads to an emittance increase and a deterioration of the beam quality. One goal of the feedback will be to prevent an excessive emittance increase by introducing additional damping.

III. THEORY AND CALCULATIONS: MODELS OF BUNCH SHAPE OSCILLATIONS
Typically, bunches in heavy-ion synchrotrons are not short compared to the rf period T RF . Therefore, the nonlinearity of the rf voltage should be taken into account during the feedback analysis or design. To the authors' knowledge, existing feedback analyses in this field are based on simplified models that linearize the rf voltage and are thus valid only for very short bunches. Therefore, Sec. III A introduces an approach based on moments. Because beam moments are not accessible for a direct measurement, Sec. III B discusses their relation to the beam current. Section IVA finally deploys the feedback model for a specific setup at the heavy-ion synchrotron SIS18 at GSI.

A. Models of bunch shape oscillations
Moments have been used for the simulation of beam dynamics in accelerators, for example, by the authors of [17,18]. It will be shown in the following that a similar technique can be used to derive feedback models. Assume a particle ensemble with N particles. The center of gravity of this ensemble in phase space can be estimated by whereas raw moments % n;m and central moments n;m of total order n þ m may be defined as x n k w m k ; n;m 1 N X N k¼1 ½x k À " x n ½w k À " w m : In particular, " x ¼ 1;0 and " w ¼ 0;1 . The uniqueness theorem as given in [19] states that under mild assumptions, a probability density function and its moment sequence are uniquely determined by each other. The requirement that the bunch shape and thus the particle density in the phase plane should be stationary therefore implies that the moments should be stationary as well. If the moment dynamics are known, the control problem is a stabilization problem.
To incorporate the nonlinear particle dynamics (3), the polynomial approximation is made, where p < 1 is the maximum order in x k . By using the trigonometric theorem sinðx À u 2 Þ ¼ sinx cosu 2 À cosx sinu 2 and Taylor series for sinx and cosx, Eq. (3b) can be rewritten in this form and this leads to x with respect to time leads to where Eq. (3a) was used. Accordingly, using Eq. (6) leads to Raw moments % n;m may be substituted by central moments k;l by [19] which finally leads to i.e. the dynamics depend on central moments up to the order p. The dynamics of the central moments are derived in the same way. This results in nonlinear functions The moment dynamics (7)-(9) for a maximum order n þ m n max represent a closed set of differential equations for p ¼ 1 only. In case of p > 1, the resulting system dynamics can be written in vector notation as This is a system of nonlinear first-order ordinary differential equations with the vector x containing the moments up to order n max and the vector x Ã containing the additional moments of higher order. In this paper, the center of gravity and the length of the bunch are of interest and the state vectors are chosen as To obtain a closed system, a truncation method is necessary that establishes a relation Introducing this into (10) finally leads to the moment dynamics that describe the shape of the bunch. We propose the assumption that the bunch density in the phase plane ðx; wÞ is approximately a Gaussian density distribution with ellipsoidal contour lines. This can be verified by the measurements presented in Sec. IV. This normalized distribution can be written in its most general form as ðx; wÞ with vectors and matrices The bunch is shifted in the phase plane by ðx 0 ; w 0 Þ and has the orientation È. The standard deviations of its half-axes are denoted by x and w . The density function is defined on a finite domain that is exactly one rf period. The higher order moments of this density function may be calculated as x n ½w À " w m ðx; wÞdxdw: Under the assumption that the density is small outside of the defined domain, i.e. ( 1 for both x > and x < À, this calculation can be performed analytically. This simplification will be assumed implicitly in the following. An obvious result is that the bunch center is Furthermore, all moments with n þ m > 2 may be expressed as functions of the moments of second order For example, 4;0 ¼ 3 2 2;0 holds. This establishes the necessary truncation relation (11).
Because of the assumption of a Gaussian density, the complex dynamical process of filamentation will not be included and the obtained approximate model (12) will not incorporate Landau damping. In the stability analysis, this will lead to conservative results with smaller stability domains. Therefore, Landau damping will be modeled approximately by an additional damping term.

B. Observation of oscillations
The moments of the bunch are not readily accessible for measurements. However, as shown in [20], the center of gravity x and the variance 2;0 can be recovered from the beam current signal Here, Q bunch is the total charge of the bunch. As illustrated in Fig. 2, this corresponds to the projection of the phase plane density on the axis x ¼ ! RF Á. For the Gaussian density, this leads to with the dc component I 0 ¼ Q bunch =T RF and the variance 2;0 . The calculation of the spectrum (5) leads to the relations To obtain a feedback system, the dynamics (12) have to be completed with a so-called output vector. This vector contains all quantities of the state vector x that are measured and that can be used for feedback. In case of the beam feedback, the measurement consists of the phase ' 1 and amplitude I 1 of the main harmonic of I b : Altogether, Eqs. (12) and (15) describe the nonlinear dynamics of the bunch shape.

C. Linearized dynamics
The calculation for a specific approximation order p can be performed using a computer algebra system such as MATHEMATICA. For zero inputs, the nonlinear system (12) and (15) has the equilibrium This means that the shape of a stationary or matched bunch is a Gaussian density with ellipsoidal contour lines centered in the origin of the phase plane ðx; wÞ and with variances 2;0 e and 0;2 e . With Eq. (14) this implies that x 0 , w 0 , and È are zero and The parameter x has a strong physical significance, because the bunch length is usually defined in terms of the standard deviation, e.g. as 2 x or 4 x . Therefore, x will occasionally be used instead of the variance 2;0 e in the following.
In general, large p > 1 are desirable to obtain dynamics (12) with low approximation errors. The comparison for different p showed that for p þ 1 > 14, the relative change of the oscillation frequencies of the bunch center and bunch length with respect to p þ 1 ¼ 14 is less then 2% for x < =2. These frequencies will be introduced in the following. Therefore, p þ 1 ¼ 14 is chosen in the following and this yields the condition This shows that for very small bunches with 2;0 e ( 1, the variances are equal and the matched shape of the contour lines are circles. This is in agreement with the fact that, provided that the axes are scaled accordingly, the trajectories in the phase plane are circles in the linear regime near the origin; cf. Fig. 1. For larger bunch sizes, the variance 0;2 e is smaller than 2;0 e . This corresponds to the flattening of the shape of the trajectories towards the separatrix. A linearization of the nonlinear dynamics (12) around the equilibrium (16) leads to the system of linear ordinary differential equations of first order with small deviations Áx ¼ x À x e and Áy ¼ y À y e from the matched shape. The linearized system falls apart into two single-input single-output systems. First, the dynamics of the bunch center " x d dt with the input u 2 and This is also the case for the nonlinear equations that will not be listed here explicitly due to space restrictions. The coefficients a i and b 1 are polynomials in x : ðÀ1Þ n n!2 n 2n x ¼ 0;2 e 2;0 e ; (19a) ðÀ1Þ n ðn þ 1Þ n!2 n 2n x ; (19b) ðÀ1Þ n n!2 n 2ðnþ1Þ It is possible to derive a transfer function directly from the representation (18); cf. [21]: notes the unit matrix with appropriate dimension. For the bunch length deviation Á 2;0 , the transfer function is obtained. One pole and zero at s ¼ 0 cancel; this is due to an invariant in the system, namely, the bunch area. This follows from a comparison of the third and fifth row of (17a), which shows that the relation holds, i.e. a relative increase in the variance of one axis implies the same relative decrease in the other axis. The poles of (20) are given by AEj! bl : For 2;0 e ( 1, the bunch shape oscillates with twice the synchrotron frequency ! syn;0 , since a 1 ¼ 1, a 2 ¼ 1, a 3 ¼ 2 holds. For larger bunch sizes, the frequency decreases. These results are in agreement with the qualitative considerations of Sec. II. In a similar way, the poles of the transfer function of the bunch center are obtained as s bc ¼ AEj! syn;0 ffiffiffiffiffi a 1 p AEj! bc ; and due to (19), the coherent frequencies of the feedback model (17) are thus given by ðÀ1Þ n n!2 n 2n x v u u t ðÀ1Þ n ðn þ 2Þ n!2 nþ1 2n x v u u t : (21b) Figure 3 shows these oscillation frequencies of the bunch center and length as functions of the bunch size. In addition, they are compared to the frequency that is obtained from Eq. (4) for the amplitudex k ¼ 2 x , i.e. this is the synchrotron frequency at twice the standard deviation of the bunch. It is similar to the coherent frequencies (21), as can be seen from Fig. 3, and will therefore be called the effective synchrotron frequency. Figure 3 also shows the coherent frequencies obtained from tracking simulations. For very large bunches ( x > 1), the deviation between theory and simulation increases. This is due to the simplifying assumption that was made in Sec. III, i.e., a small particle density for jxj > .
As mentioned above, the transfer function (20) does not include inherent damping (Landau damping), because filamentation has been neglected due to truncation (11) with the assumption of a Gaussian density. A more general transfer function with damping is given by with the damping coefficient d 2 ½0; 1. For d ¼ 0, transfer function (20) is obtained. The poles are given by i.e. the oscillation frequency is ! bl , independent of d.
Transfer function (23) can also be written as and this shows that the relative voltage modulation u 1 leads to a relative bunch length modulation. Figure 4 shows simulation results comparing transfer function (23) with a nonlinear tracking simulation. The black curve shows the tracking result with about 4:8 Â 10 5 macroparticles. The bunch is initially matched at a voltage amplitudeÛ 0 of 5 kV with a bunch size of x ¼ 1:135 rad; this corresponds to a bunch length of 4 x ¼ 19:6 m. At t ¼ 0, the voltage is raised to 10 kV to induce bunch length oscillations. Without filamentation, i.e., assuming a constant bunch area, the new equilibrium would bẽ Using these values, the initial value x is obtained and the resulting response of the linear system Á 2;0 ðtÞ leads to the bunch length which is shown as the blue curve in Fig. 4. Because of filamentation, the result of the tracking simulation has a 4% larger final bunch length x;track ¼ 0:991 rad. This difference in the final value is of minor importance for the feedback analysis, because the feedback loop rejects the dc component. The comparison with transfer function (23) shows a good agreement for a damping factor of d ¼ 0:15. The damping factor has been chosen such that similar envelopes are obtained for both curves. Note that both the damping and the frequency of the nonlinear simulation increase with time and are thus time dependent. This is a clear demonstration of the fact that filamentation is a nonlinear process. All necessary simulation parameters are given in Table I; the bunch sizes are summarized in Table II.

A. RF setup of SIS18
A simplified scheme of control loops in the heavy-ion synchrotron SIS18 at GSI is shown in Fig. 5. A complete description of the rf systems is given in [22]. An amplitude modulator drives the accelerating cavity that produces the voltage U RF with amplitudeÛ 0 and amplitude and phase modulations u 1 and u 2 ; cf. (2). The voltage determines the longitudinal dynamics of the beam. A beam position monitor is used to measure the beam current I b .
In the control loop for the stabilization of the beam phase, the phase difference " x À u 2 between the bunch  and the voltage U RF is determined. The controller consists of a finite impulse response (FIR) filter that is implemented on a FPGA. The time delay introduced by this digital hardware is modeled as the time delay T d . A direct digital synthesizer is used to produce the analog phase modulation signal u 2 and can be modeled as an integrator with negative sign as described in [23]. The control loop for the bunch length has a similar structure. The only difference is that the amplitude I 1 of the first harmonic is used and instead of the direct digital synthesizer, an integral controller is included in addition to the FIR filter. The amplitude is normalized by its equilibrium I 1 e ; cf. (16b).
In the following analysis we will focus on the control loop for the bunch length. Using (15) and (16b) yields ÁI 1 =I 1 e ¼ ÀÁ 2;0 e =2; and using (23) leads to the linearized control loop as shown in Fig. 6. The loop contains a FIR filter y fir ðtÞ ¼ À 1 4 I 1 ðtÞ þ 1 2 I 1 ðt À 1=2f pass Þ À 1 4 I 1 ðt À 1=f pass Þ I 1 e : This FIR filter was first designed and implemented on FPGA for the bunch phase loop [23]. It rejects the dc component of the measurement and has a bandpass characteristic at the frequency f pass , which is typically chosen close to the frequency of the bunch shape oscillation. Because of the dc rejection, it is possible to use the small-signal variable ÁI 1 in Fig. 6. For the implementation of the filter, a DSP system is used consisting of analog narrowband preprocessing, automatic gain control, analog-to-digital converters, FPGA, DSPs, and digital-to-analog converters. More details on this hardware can be found in [22], p. 6. The FIR filter is realized on a DSP system with floating point arithmetic, enabling almost arbitrary coefficients.
The filter transfer function is Thus, the feedback loop is a linear time-delay system [24]. The frequency response of the filter is and this is a bandpass filter with a first maximum at ! ¼ 2f pass . The phase at this frequency is À. In the time domain, the total length of the filter is 1=f pass . This corresponds to a discrete filter length (total number of taps) of N fir ¼ f samp =f pass with the sampling frequency f samp ¼ 375:44 kHz. Because the length N fir is an even integer, the filter frequency f pass is discrete as well and this will become apparent for the simulation results in the next section; cf. Fig. 9. Two parameters are available for the tuning of the feedback. On the one hand, the integrator gain K I determines the open loop gain. On the other hand, the filter frequency f pass , i.e., the filter length, changes the open-loop phase by À!=2f pass and the gain.
To obtain the stability domain depending on the degrees of freedom f pass and K I , a Nyquist-type stability criterion is used. The open-loop frequency response is given by with the normalized gain Its Bode plot is shown in Fig. 7 for K I;norm ¼ À0:1425, f pass ¼ 9 kHz, T d ¼ 10 s, d ¼ 0:15, and ! bl ¼ 28:9 kHz. This corresponds to the settings of experiment C, cf. Table II, and a ratio of 1.95 between filter frequency and bunch oscillation frequency. The amplitude margin of the feedback loop equals 14.9 dB. A necessary condition for a transition between stability and instability is the crossing of the Nyquist plot with the critical point À1: jG ol ðj! cr Þj ¼ 1; ∡G ol ðj! cr Þ ¼ Àp; where ! cr denotes the critical frequency at the crossing and p 2 f1; 3; 5; . . .g is an odd integer. The gain at the crossing will be denoted by K cr . In the following, boundaries K cr > 0 and K cr < 0 will be determined such that stability corresponds to K cr < K I;norm < K cr . Using the normalized variables leads to the condition jK cr j ¼ jj for the amplitude. Considering only positive frequencies ! cr > 0, the phase condition can be written as with > 0 and the abbreviations For the special case d ¼ 0, the phase condition can be solved explicitly for and introduced into (25) to obtain K cr as a function of and different p. The stability limit for a specific is then the smallest absolute gain with respect to p. This analysis was used in [23,25] to obtain stability diagrams. It turns out that for a specific , the feedback is only stable for either positive values of K I;norm up to a certain maximum or for negative values up to a minimum value. The distinction of these two cases can be easily made by analyzing the direction of traverse of the Nyquist plot or by calculating the net angle that is circumscribed by the Nyquist plot with respect to the critical point À1.
For d Þ 0, a numerical evaluation is necessary. In addition, stability for a specific is now possible for both positive and negative gain values. The stability borders are again obtained by finding the smallest absolute gain with respect to different p. The result is shown in Fig. 8. As could be expected, a nonzero damping coefficient d increases the stability domain.
Because of the modeling assumptions and the linearization, a comparison of the previous results with nonlinear tracking simulations is useful. To evaluate the feedback performance, the settling time T settling of the closed loop is calculated. This settling time is defined as the time for which the bunch length 2;0 ðtÞ finally remains inside a 10% tube around its equilibrium value 2;0 e , i.e.
2;0 ðtÞ À 2;0 e 2;0 e <0:1 for all t > T settling holds, where 2;0 e was chosen as the mean of 2;0 ðtÞ. Figure 9 shows the result of a simulation scan based on nonlinear tracking with a simulation length of 400 turns. The parameters of the simulation are again given in Table I. The performance is visualized with a color scale from blue (unstable, no settling) to red (stable, smallest settling time). The beam experiment with largest damping (run C; cf. Sec. IV B) is indicated by a Â and is inside the region of best performance of the tracking simulations. In addition, Fig. 9 shows the experimental results of 14 additional operating points. It can be observed that the nonlinear simulations and most experimental results agree well with the stability limits for d ¼ 0:15.

B. Beam experiment
A beam experiment will be used in this section to verify the newly developed theory and the simulation results. The experimental setup is described in detail in [25]. The parameters of the experiment dated October 24, 2007, are identical with Table I. During the experiment, several runs were performed (cf. Fig. 9). Three runs will be analyzed in more detail in this section and they will be denoted by A, B, and C. Run A represents the open-loop case, where only Landau damping is present. During run B, the bunch length feedback was switched on, and during run C both the bunch position and bunch length feedback were active.
The bunch sizes of the three runs have minor variations that are summarized in Table II. One reason for these variations are fluctuations in the ion source that cannot be avoided completely. The bunches initially have the equilibrium size x . At time t ¼ 0, the rf voltage ampli-tudeÛ 0 is doubled stepwise to induce bunch length oscillations. This leads to the new equilibrium sizes x that are summarized in Table II. Figure 10 shows the beam current signal I b for run C at the time instants t 1 ¼ À0:2 ms, t 2 ¼ 0:112 ms, and t 3 ¼ 0:6 ms. As initial conditions of the simulations, particle ensembles with Gaussian densities as (13) are used. The size x is obtained by fitting the density to the measurement of the beam current I b at t ¼ À0:2 ms. The size w follows then from (19d) and the fact that 0;2 e ¼ 2 w . The comparison between measurement and simulation at t 1 shows that the assumption of Gaussian densities is approximately justified. The simulation results for t > t 1 are obtained by nonlinear tracking simulations using the initial particle ensemble. A good agreement for both the amplitude I 1 and beam current I b can be observed.
The measured amplitude I 1 ðtÞ is shown in Fig. 11 for all three runs. A closed-loop frequency domain graph of I 1 can be found in [13]. In addition, the measurements are compared to the corresponding nonlinear tracking simulations. In the simulations, a time constant of T cav ¼ 20 s is taken into account for the voltage jump ofÛ 0 . This time constant is due to the cavity and its low level radio frequency system. The feedback gain in the simulations equals K I;norm ¼ À0:1425 and the filter frequency is f pass ¼ 9 kHz, which corresponds to ¼ 1:95. This point is marked in Fig.s 8 and 9 by a Â and is clearly in the stable region.

V. TUNING RULE
For future automatic operation for experiments with different ion species and bunch sizes, a tuning rule is necessary for the adjustment of the feedback parameters f pass and K I . Based on tracking simulations, the tuning rule f pass 2f syn;0 % 1:1; is proposed for bunch sizes x ¼ 0:4 . . . 1:1 rad, which corresponds to a AE2 bunch length between 90 and 250 . The central control room may deliver the synchrotron frequency f syn;0 and an estimation of the bunch length x . For the normalized variables of the stability diagram, this implies  Fig. 12 with the Nyquist stability limits for d ¼ 0:045.
As an alternative to the tuning rule, the linear feedback model of Fig. 6 can be used to obtain a good estimate for the optimum settings. To show this, the linear feedback model was discretized with respect to time, which leads to a linear model of higher order that includes the delays of the FIR filter. The eigenvalues of the resulting system matrix for a given filter can be calculated numerically and the distance of the largest eigenvalue to the unit circle is an estimate for the damping of the closed loop. This distance is shown in Fig. 13 for scan 1, where red is the most stable configuration and dark blue is unstable. The result is in very good agreement with the Nyquist stability limits. It can also be observed that the optimum of Fig. 13 is very close to the optimum of  linear model for all scans are given in Table III. For very large bunch sizes (scan 3), the estimate of the linear model becomes less precise, whereas the tuning rule is quite accurate for all scans.

VI. DISCUSSION
In Fig. 11 the first maximum of run A occurs at t % 0:11 ms. Taking into account also the time constant T cav , this implies an oscillation period of about 2 Á 0:11 ms À 20 s ¼ 0:2 ms. According to the developed model and the given bunch size in Table II, the frequency of the bunch length oscillation is ! bl ¼ 1:463 Á ! syn;0 . This implies a period of 0.206 ms, which is close to the observed value.
The feedback algorithm consisting of the FIR filter and the integral controller is implemented in a discrete way on a combined DSP/FPGA system [13]. From (24) and (26), the discrete integral gain is obtained, where f samp ¼ 375:44 kHz is the sampling time of the feedback loop. Figure 9 shows that the stability limits for d ¼ 0:15 agree very well with the nonlinear tracking simulations. This verifies that Landau damping may approximately be described by an exponential damping term. Expectedly, the effect of this damping is a considerably larger region of stability. An important practical implication of (26) is that the open-loop gain depends on the bunch size. It is thus necessary to adapt the gain to this parameter or to provide a sufficient stability margin.

VII. CONCLUSIONS
A new approach for the modeling of bunch length oscillations in synchrotrons has been proposed based on a moment method. With an additional damping term, the obtained models are valid also for bunches with a large length compared to the period of the rf voltage. The models generalize existing models in literature that are typically based on linearized single-particle dynamics. The analysis in this paper shows how the oscillation frequencies of the bunch modes decrease for larger bunch sizes.
The models have been used to analyze the stability of an rf feedback setup at GSI. Experimental results are in agreement with the simulations and the developed theory. The proposed approach is useful for further optimizations of the rf feedback loops. In particular, the models can be used for the design of more complex feedback algorithms or the robustness analysis of the existing feedback with respect to changing parameters during the acceleration of the beam. In this paper, only the stationary case with a single-harmonic rf voltage has been studied, but the approach may be applied for the acceleration case and for the case of more complex rf voltages.