Multibunch and multiparticle simulation code with an alternative approach to wakefield effects

The simulation of beam dynamics in the presence of collective effects requires a strong computational effort to take into account, in a self-consistent way, the wakefield acting on a given charge and produced by all the others. Generally this is done by means of a convolution integral or sum. Moreover, if the electromagnetic fields consist of resonant modes with high quality factors, responsible, for example, for coupled bunch instabilities, a charge is also affected by itself in previous turns, and a very long record of wakefield must be properly taken into account. In this paper we present a new simulation code for the longitudinal beam dynamics in a circular accelerator, which exploits an alternative approach to the currently used convolution sum, reducing the computing time and avoiding the issues related to the length of wakefield for coupled bunch instabilities. With this approach it is possible to simulate, without the need for large computing power, simultaneously, the single and multibunch beam dynamics including intrabunch motion. Moreover, for a given machine, generally both the coupling impedance and the wake potential of a short Gaussian bunch are known. However, a classical simulation code needs in input the so-called “ Green ” function, that is the wakefield produced by a point charge, making necessary some manipulations to use the wake potential instead of the Green function. The method that we propose does not need the wakefield as input, but a particular fitting of the coupling impedance requiring the use of the resonator impedance model, thus avoiding issues related to the knowledge of the Green function. The same approach can also be applied to the transverse case and to linear accelerators as well.


I. INTRODUCTION
When dealing with the study of high intensity beam dynamics, we must take into account, in a self-consistent way and in addition to the external guiding fields, the effects of the self-induced electromagnetic fields, called wakefields [1]. The theory of collective beam instabilities induced by the wakefields is a broad subject and it has been assessed over many years by the work of several authors, such as Sacherer [2], Chao [3], Laclare [4], Zotter [5], Pellegrini [6], Sands [7] and others [8]. In this paper we will focus only on the longitudinal beam dynamics in circular machines, even if the proposed approach can be applied to the transverse case and to linear accelerators as well.
To simplify the theoretical study, in general it is convenient to distinguish between short range wakefields, which influence the single bunch beam dynamics, and long range wakefields, where high quality resonant modes excited by a train of bunches can last for many turns exciting, under some conditions, coupled bunch instabilities. In both cases the bunch distribution is considered as a sum of coherent oscillation modes perturbed by the wakefields. A linear perturbation theory is then generally used to study the resulting instabilities.
If we want to analyze the behavior of the beam under the influence of wakefields also in the nonlinear regime, a simulation code is a valid approach [9]. The equations of motion which are integrated in this case are quite simple, and the main issue is related to the inclusion of the wakefield effects avoiding the introduction of numerical noise and nonphysical phenomena. There exist several simulation codes specialized in circular accelerators for the single bunch longitudinal [10][11][12] and transverse [13,14] beam dynamics, and for multibunch [15,16], which have been proved to be reliable tools in the comprehension of the wakefield effects.
In order to include the wakefield in a simulation, we need a convolution sum taking into account the electromagnetic fields acting on a charge and produced by all the others. There arises a series of issues, from the necessity of a strong computing power for coupled bunch simulations to the problem of knowing the Green function of an accelerator instead of the more commonly known wake potential.
In this paper we present a novel simulation code for the longitudinal beam dynamics in a circular machine, which uses an alternative approach with respect to the previously mentioned codes, allowing one to simulate, simultaneously, the effects of short and long range wakefields without the necessity to distinguish between the two, and avoiding the issues related to the use of wakefield and convolution sum. The algorithm is quite light and does not need high computing power, even with the inclusion of intrabunch motion.
In the following section we recall the common approach used in simulation codes to take into account the effects of the wakefields, and the main issues arising with such a method. In Sec. III we describe the new multibunch/particle simulation code MUSIC and its main features, and in Sec. IV we benchmark the simulation results against theories and another code which uses the classic wakefield convolution sum, either for single and multibunch cases. Finally, concluding remarks and outlook will end the paper.

II. COMMON APPROACH IN WAKEFIELD SIMULATIONS
We start our analysis by considering the longitudinal equations of motion of a single particle in a circular accelerator and by taking into account the effect of the wakefields. We suppose that the energy exchange between the charge and the surrounding environment is localized in a single point in the machine and that the minimum time step is the revolution period T 0 . However, the same procedure can be extended to cases in which it can be necessary to update the longitudinal dynamics several times per turn, as, for example, in a strong rf focusing regime [17], or in some cases when the longitudinal broadband impedance due to space charge requires it. By supposing β ¼ 1, we can then write where z is the longitudinal position of a single particle with respect to the synchronous one, and positive behind it, L 0 the machine length, η the slippage factor, ε the relative energy variation of the charge with respect to the nominal energy E 0 , V rf the rf voltage, U 0 the energy lost per turn by the synchronous particle, V wf the wakefield effect, RðT 0 Þ and τ s respectively a stochastic variable which takes into account quantum fluctuations and the longitudinal damping time, two terms related to synchrotron radiation effects and important, in particular, in electron machines.
The common approach used in a simulation code models the single bunch evolution as an ensemble of particles each one ruled by the above two coupled equations. Since a bunch consists of 10 8 -10 12 and even more charges, it requires a strong computing effort and a parallel cluster to keep track of all the particles. In general macroparticles, gathering a given number of charges, are used, and this permits one to reduce the number of equations to be solved step by step.
The term V wf , which has to take into account the wakefield acting on a charge and produced by all the others, depends on the longitudinal bunch distribution λðzÞ according to the relation with Q tot the total charge of a bunch, and w ∥ ðzÞ the wakefield of a point charge. The convolution integral is also called wake potential. It represents the energy gained or lost by a charge due to the entire bunch. In writing the equation we have used the causality property that w ∥ ðzÞ ¼ 0 for z < 0, and we have ignored the effects of long range wakefields. In order to also include these, we have to add a sum over previous bunches and turns to the above convolution integral. Let us consider, for the moment, only the single bunch effects. In order to include the convolution integral in a simulation code using N m macroparticles, we have to transform it into a sum, such that Eq. (3) becomes with z i the position of the ith macroparticle. Equation (4) needs some comments. It has to be evaluated at each turn for each macroparticle, and it depends on the position of all the others preceding it. This means that, at each turn, the effect of wakefield requires about ðN m − 1ÞN m =2 operations. If we want to simulate 10 6 macroparticles, which is a reasonable number to have reliable results, at each turn we need to do more than 10 11 operations, which can be possible nowadays only on parallel computing clusters.
In order to reduce the computing time and only for the evaluation of the wakefield effect, the bunch is generally divided into N s slices, or bins, of width Δ and with center z iΔ , each one containing n i ðΔÞ macroparticles. By supposing that slices act as point charges, the wake potential at the center of each slice is then evaluated by using the following relation: Once the wake potential is known in the positions z iΔ , a linear interpolation permits one to evaluate the wake potential acting on any macroparticle of the bunch. Since in general the number of slices is between few hundreds to some thousands, this reduces quite a lot the number of operations.
This kind of approach has been widely used in simulations, and through the years it has proved to give reliable results. However, some cares have to be taken when deciding the size and the number of the slices (and, of course, of macroparticles). A low number of slices reduces the computing time, but it could suppress some microstructures in the bunch leading to instabilities. Slices represent a nonphysical artifice, which can introduce additional numerical noise to that of macroparticles, making necessary, in some cases, a parametric study. Once the main parameters of a simulations have been decided, it is important to investigate any dependence of the results on the number of slices and of macroparticles inside the slices. As an example, let us consider the wakefield of a broadband resonator given by [1] 2Q z=c cosðω n z=cÞ− ω r 2Qω n sinðω n z=cÞ HðzÞ; where ω r is the resonant angular frequency, R the peak impedance, Q the quality factor, c the speed of light, ω n the natural angular frequency and HðzÞ the Heaviside step function. By using ω r ¼ 15 × 10 9 rad=s and Q ¼ 1, reasonable parameters for an accelerator impedance model, and a Gaussian bunch with an rms length of 2 cm, we show on the left side of Fig. 1 the normalized wake potential obtained with the convolution sum of Eq. (5) with 100 slices compared with the analytical solution of the convolution integral of Eq. (3). The two results coincide. However, if we take a bunch length of 2 m (as, for example, a proton bunch of the CERN Proton Synchrotron [18]), as shown on the right side of Fig. 1, 1000 slices=sigma, which in general represents quite a high number, are not enough and the convolution sum gives a wrong result. Only when the number of slices is higher than 10 4 , the analytical solution and the slice reconstruction of the wake potential are close to each other, but still they do not coincide. This means that if a bunch is very long compared to the wakefield variations, a very high number of slices could be necessary, thus making almost worthless the use of slices.
Another issue, which arises by the use of Eq. (5), is related to the knowledge of the Green function, that is the wakefield produced by a point charge. In order to obtain the impedance model for an accelerator, in general the theoretical approach requires one to solve the Maxwell's equations in a given structure taking the beam current as the source of fields. This could become a very complicated task for some accelerator devices. Therefore dedicated computer codes have been developed, which solve the electromagnetic problem in frequency or in time domain. There are several powerful codes for the electromagnetic design of accelerator devices, such as CST PARTICLE STUDIO® [19] and GDFIDL [20]. However these codes give, as result, the coupling impedance, which is the Fourier transform of the wakefield, and the wake potential of a Gaussian distribution, but not the Green function. In order to obtain the wakefield needed in a simulation code, one has to do numerically the inverse discrete Fourier transform of the impedance, which does not necessarily turn out to be a real function as it should be, and some artificial modifications of the results must be devised to obtain a function that could be used in Eq. (5), with the risk of producing nonphysical results.
A last, but just as important, issue is related to the study of coupled bunch instabilities by taking into account also the internal motion of the particles. This could be of interest for quadrupolar instabilities, and when we want to account for the effect of rf nonlinearity and bunch shape. The instabilities are generated by resonant modes, given by Eq. (6) with a high quality factor Q. Depending on its value and on the resonant frequency, these modes can last for several hundreds of nanoseconds up to microseconds, influencing many bunches for many turns. In this case we need to know the wakefield with a small step and until it becomes zero. This slows down the simulation making it necessary to parallelize the code [16].

III. MUSIC code
In order to deal with the last issue of the previous section, which is the coupled bunch instability, let us suppose, for the moment, that each bunch is a rigid macroparticle performing dipolar oscillations, and a single resonant mode, with the wakefield given by Eq. (6), is excited.
A first bunch with total charge q 1 and position z 1 , interacting with the mode, loses an energy that, due to the fundamental theorem of beam loading [21], is equal to A second bunch of charge q 2 , at a distance z 2 behind the first one, loses an energy due to the wakefield generated by the first one and to the beam loading theorem, and so forth for all the following bunches. For example, for the nth bunch at a distance z n from the first one, we have This is exactly the procedure leading to the convolution sum for coupled bunch beam dynamics. With N b bunches we need to evaluate the wakefield, at each turn, The same result can be obtained with a matrix formalism and a lower number of operations. As shown in the Appendix, we can define a "wakefield" matrix [22] of the kind cosðω n z=cÞ − ω r 2Qω n sinðω n z=cÞ − Rω r Qω n sinðω n z=cÞ ω r Q Rω n sinðω n z=cÞ cosðω n z=cÞ þ ω r 2Qω n sinðω n z=cÞ and a "charge" matrix q as to obtain the value of the wakefield of a resonant mode at any position behind a charge that excites the mode itself. Let us suppose that a charge q 1 , in a position z 1 , is exciting the resonant mode initially unloaded. After its interaction, according to the results of the Appendix, the initial conditions for the matrix evolution can be written as where we have put q 1 ¼ q 0 as new initial conditions. In this way, the value of the wakefield in any position behind the charge can be obtained from the 2 × 1 matrix whereqðzÞ ¼ MðzÞ · q 0 .
If a second charge q 2 is following at a distance z 2 the first one, the energy change of Eq. (8) can also be obtained as where with the notation ½qðz 2 − z 1 Þ 1 we mean the first term of theq matrix. At this point we can continue the evolution of the wakefield by considering, as new initial conditions for the matrix evolution, the quantity which takes into account the voltage across the capacitance and the current in the inductance at the passage of the second charge, plus the excitation produced by the second charge itself. A third charge at a distance z 3 from the first one has then an energy variation due to the wakefield that can be written, starting from these new initial conditions, as and, after its passage, the new initial conditions become and so forth. The advantage of this formalism is that it is not necessary any more to use the convolution sum and, therefore, to keep track of the entire wakefield, but we can just transport the matrixq from one charge to another and add the new initial conditions. This reduces the number of operations to evaluate the wakefield for This approach has been already used to study the longitudinal coupled bunch instabilities in DAΦNE [23] with the inclusion of a time domain feedback system [24]. The code can simulate only a dipole instability because it models a bunch as a single rigid macroparticle, without the possibility to study intrabunch motions. It could also simulate a gap in the bunch train or any filling pattern of the machine.
If we now consider the bunches divided in macroparticles, the same procedure can also be applied between two macroparticles of the same bunch, and between the last particle of one bunch and the first particle of the following one (or the same bunch at the following turn). This idea has been implemented in the new code MUSIC (multibunch/ particle simulation code), which extends the above approach to include multiparticles for each bunch. In this way we avoid the issue related to the length of wakefield when dealing with coupled bunch instabilities, and we do not need to use slices any more: a simulation with N m macroparticles requires N m − 1 operations plus N m logðN m Þ for sorting on the z variable. It is important to note that, as long as the longitudinal positions of the particles do not change too much between updates, this operation is not necessary at each turn. In this way the code can simulate simultaneously single bunch and multibunch effects without the necessity of parallel computing.
Another advantage of the matrix transport formalism of the wakefield implemented in MUSIC is its extension to include a generic impedance, not necessarily a resonant mode. In the above procedure, if we have more than a single resonator, the superposition principle allows us to repeat the same method for all the modes, each one evolving separately from the others with its own wakefield matrix. We can also use the same idea to fit a generic impedance with a sum of resonant modes. In this way we do not need, as input, the wakefield, but, once we have the machine coupling impedance, we can fit it with a finite sum of different resonant modes and use them with the transport matrix formalism. In principle there is no limit to the number of modes that can be used for the fit.
The coupling impedance of a resonant mode, obtained by the Fourier transform of Eq. (6), and given by follows the properties of a generic coupling impedance, that is its inverse Fourier transform is real, and it respects the causality principle for which is equivalent to the Hilbert transforms [3] where P:V: means taking the principal value of the integral. The above expressions indicate that, in principle, knowing either the real or the imaginary part of a coupling impedance, one can construct the whole impedance. This observation is important when we want to determine the fit of an impedance as a sum of resonant modes, since it is sufficient to fit only its real or imaginary part. Moreover, we also notice that, since in general the integrals in Eq. (20) may be difficult to evaluate, and the operation should be applied with care [3], if we know the real or the imaginary part of a coupling impedance, and we fit it with a sum of resonant modes, the other unknown component of the impedance must have a behavior very close to the corresponding one determined by the sum of resonators. This method could then be used to obtain, with some approximation, the whole impedance, once we know either its real or imaginary part.

IV. ALGORITHM BENCHMARKING
The transport matrix of the wakefield described in the previous section is basically the idea behind the algorithm implemented in the simulation code MUSIC. There is no need of using slices, nor wakefield with the problem of the total length and minimum step. Moreover, it is not necessary to do the inverse Fourier transform of the machine coupling impedance. We need to fit it with a sum of resonant modes, and to propagate the wakefield matrix, for each mode, between two consecutive macroparticles. It is also possible to take into account, at the same time, short and long range wakefields by using, in the same simulations, broadband and narrow band resonators.
In this section we compare the results of MUSIC, under some particular conditions, with the theory of coupled bunch instability, we test a feedback system to cure the instability, we check the possibility to use the code for single bunch simulations by benchmarking it against a simulation code which uses slices and the convolution integral [18], and we also study the effect of the microwave instability in electron machines compared with the mode coupling theory. We have to remember, in addition, that all the above effects can be simulated simultaneously with MUSIC.

A. Coupled bunch instability and multibunch modal analysis
The analysis of the beam dynamics of N b equally spaced bunches interacting with long range wakefields is performed by computing the coherent synchrotron frequency shift predicted by Sacherer's theory [2]. Starting from the linearized Vlasov equation, and developing any perturbation of a bunch stationary distribution as sum of coherent azimuthal oscillation modes (dipolar, quadrupolar, …), it is possible to obtain an eigenvalue system representing the coherent frequencies of the modes.
By neglecting the coupling between the different azimuthal modes, hypothesis valid, for example, when the beam intensity is weak, and by considering the case in which there is only a single high quality resonator as source of instability, such that the impedance given by Eq. (18) is always zero except around its resonant frequency, which we suppose to be close to an integer number of the accelerator revolution frequency pω 0 AE mω s , with ω s the synchrotron angular frequency and m a fixed azimuthal mode, it is possible to show that, if p is not a multiple of N b , there are two main multibunch oscillation modes [3], one stable and the other unstable, and the corresponding exponential growth rate and damping rate are given by with N p the number of particles per bunch, and G m ðxÞ a form factor depending on the bunch stationary distribution. For a Gaussian bunch of rms length σ z , we have with I m the modified Bessel function of the first kind, while for a parabolic distribution of total length 2z max , we have with J m the Bessel function of the first kind. In Eq. (21) p can be obtained in two ways, corresponding to the two multibunch oscillation modes: identifying the mode number, ranging from 1 to N b − 1. The stable mode, above transition energy η < 0, is given by the lower sign and μ 1 , while the unstable mode is given by the upper sign and μ 2 . The opposite happens below transition.
Since MUSIC is a time domain code, it can simulate the coupled bunch oscillations giving, for every bunch and for every particle inside the bunch, turn after turn, the values of Δz and Δε of Eqs. (1) and (2). In order to compare the theoretical growth and damping rates of stable and unstable multibunch oscillation modes with the results of MUSIC, we need to perform a multibunch modal analysis from the data of the code. For simplicity let us assume that we want to know, at every turn, the amplitude and phase of the μ ¼ 0; …; N b − 1 multibunch oscillation modes given by the dipolar oscillations of the centers of mass of the bunches (m ¼ 1). Since the oscillations of each bunch must be reconstructed with these modes, we have a system of N b equations of the kind with b ¼ 0; …; N b − 1 identifying the bunch, A b and θ b the known oscillation amplitude and phase of the bunches, A μ and θ μ the unknown amplitude and phase of the multibunch modes. The phase factor 2πμb=N b identifies the "pattern" of a given multibunch mode. For example, the mode μ ¼ 0 represents all the bunches oscillating in phase. The above equation system, representing the modal reconstruction of the center of mass oscillations, is valid in the approximation that bunches have the same synchrotron frequency, even for uneven fills. The values of A μ and θ μ of the dipolar oscillations for all the multibunch modes at each turn are calculated by MUSIC. Let us first consider a simulation with all the bunches modeled as simple rigid macroparticles performing dipolar oscillations (m ¼ 1). As a reference machine, we take the CERN Proton Synchrotron (PS) at a fixed energy above transition with harmonic number h ¼ 7 [25]. For the wakefield, we use two resonant modes, one with Q 1 ¼ 500, R 1 ¼ 5 kΩ, ω r1 ¼ 24 × 10 6 rad=s, which corresponds to p ¼ 8, it damps the multibunch mode μ 1 ¼ 6 and excites the mode μ 2 ¼ 1 (p 1 ¼ 2), and the other with Q 2 ¼ 300, R 2 ¼ 3 kΩ, ω r2 ¼ 30 × 10 6 rad=s, which corresponds to p ¼ 10, it damps the multibunch mode μ 1 ¼ 4 and excites the mode μ 2 ¼ 3. The growth rates and the damping rates of the dipolar modes, in the linear approximation, can be obtained with Eq. (21). Regarding the form factor G m ðxÞ, for a point charge, in the limit of x → 0, both Eqs. (22) and (23) give G 1 ¼ 1.
The other parameters used in the simulation, and the results compared with the theory, are given in Fig. 2. In the figure we show the amplitude of the seven multibunch oscillation modes as a function of the number of turns. The amplitude is expressed in terms of rf phase. The modes 1 and 3 are unstable, increasing exponentially with the turns, while the corresponding modes 6 and 4 are damped. The results are in very good agreement with the growth and damping rates predicted by linear theory of Eq. (21). All the other oscillation modes are unaffected by the two resonators.

B. Coupled bunch instability and Landau damping
In order to test also the multiparticle internal distribution, we have simulated the same case represented in Fig. 2 but with a Gaussian distribution. In Fig. 3 we show the amplitude of the multibunch oscillation modes of the center of mass of the bunches as a function of the number of turns. On the left side of the figure we have used a bunch length of 0.07 rad in term of rf phase (rms value). With this bunch length, from Eq. (22), the two form factors, evaluated at the frequency of each resonant mode, are equal to 0.994 and 0.99, so that we expect a behavior very similar to that of a point charge. Actually, even with these correction factors, there is a small difference between simulation and theory, clearly visible in the two growth rates. In order to investigate this difference, we repeated the simulation with a longer bunch, with rms of about 0.21 rad, which gives G 1 ðxÞ ≃ 0.944 and G 1 ðxÞ ≃ 0.915, still very close to 1. As shown in the right side of the figure, in this case all the modes are damped, even if the linear theory predicts instability.
The reason of the stable behavior of the beam lies in the fact that the nonlinearity of the rf voltage, not included in the determination of Eq. (21), produces a spread in the synchrotron frequencies of the bunch particles which stabilizes the multibunch coherent oscillations through the Landau damping mechanism [26]. Indeed, if we use in MUSIC a linear rf voltage instead of a sinusoidal one, with exactly the same parameters used for the simulations of the right side of Fig. 3, we get the results shown in Fig. 4, now in very good agreement with the linear theory.
The effect of the Landau damping due to the nonlinearity of the rf voltage and produced by a resonant mode, in case of N b equally spaced bunches, by neglecting the coupling between different azimuthal modes m, and for a Gaussian distribution, can be evaluated from the dispersion integral [27] 1 ¼ i and h the rf harmonic number. The integral (25) has to be executed in the complex x plane by deforming the contour of the integration, in order to avoid the singularity. The corresponding stability diagram of the real and imaginary part of the impedance for the longitudinal dipolar oscillations ðm ¼ 1Þ, with the resonant mode frequency at ω r ¼ 8ω 0 þ ω s and the parameters of the previous simulations, is found by imposing the imaginary part of y equal to zero, and it is shown in Fig. 5. On the left side we have obtained the threshold curve for the Gaussian bunch with rms length equal to 0.07 rad, and on the right side the same plot for 0.21 rad bunch length. The intersection of the threshold curve with the real axis of the impedance gives the maximum value of the impedance the resonant mode can have, such that the instability growth rate can be contrasted by the Landau damping.
For the shortest bunch length, the threshold at the resonant frequency gives a maximum impedance of about 700-800 Ω, while for the longer one we get about 7 kΩ. Since in the simulations we have used R ¼ 5 kΩ, this explains the Landau damping for the longer bunch.
In order to check the threshold limit of the right side of Fig. 5, we have performed simulations with one single resonant mode at ω r ¼ 8ω 0 þ ω s with three different values of R (6 kΩ, 7 kΩ, and 8 kΩ), which correspond to the cases below, at limit, and above threshold. Since this resonator excites the mode μ ¼ 1, in Fig. 6 we show the amplitude of only that multibunch oscillation mode as a function of the number of turns for the three cases as resulting from the simulations. The mode is clearly Landau damped for R ¼ 6 kΩ while for R ¼ 8 kΩ there is instability. For R ¼ 7 kΩ instead, we are at the limit of the instability. It is not very clear if the growth of the mode amplitude is saturating.

C. Feedback system
MUSIC has been developed from the multibunch simulation code used to study the longitudinal coupled bunch instabilities in DAΦNE, and then it already includes a time domain feedback system (bunch by bunch) for the dipolar oscillations as described in Ref. [24]. However, another way to cure the instability is to use a frequency domain feedback system, such as, for example, the one used to cure the coupled bunch instabilities in the CERN PS [28], which acts, instead of on the oscillations of the bunches, on those of the multibunch modes.
In both approaches to cure dipolar oscillations, with a time domain or frequency domain feedback, the system must give an energy kick to every mode or bunch that has to be damped, with a phase such that the optimum damping effectiveness is achieved when the kick is in quadrature with the phase oscillations [29], and an amplitude proportional to the bunch or mode oscillation amplitude. For both cases the damping rate provided by the feedback is given by [29] α with ω rf the rf angular frequency, and g the feedback gain expressed in V=rad, and related to the energy kick by In the above equation ϕ can represent the phase oscillation (in rad) of either the bunch or the mode. If, for example, we consider a single mode μ, from Eq. (24), we get From the modal analysis described in Sec. IVA, MUSIC evaluates, turn after turn, the amplitudes of oscillation of the modes and of the bunches, so that it can correctly evaluate the kick provided to the bunches, once the gain has been set, acting as either a bunch by bunch or a mode by mode feedback system.
In order to test the validity of the code with the feedback system to cure the coupled bunch instabilities, we have used the same simulation represented on the left side of Fig. 3 but with the frequency domain feedback acting on the modes μ ¼ 1 and μ ¼ 3. The feedback gain for each mode has been chosen such that the damping rate given by Eq. (27) balances exactly the growth rate excited by the resonant modes. In Fig. 7 we show the amplitude of the multibunch oscillation modes as a function of the number of turns given by MUSIC in presence of the feedback system and with a Gaussian distribution with rms length equal to 0.07 rad.
The results are compared with the case without feedback. The modes 1 and 3 are now at the equilibrium as if they were not excited. Actually this is a dynamical equilibrium because it has been obtained by the two opposite effects: the excitation of the resonant modes and the damping of the feedback. All the other modes are unaffected.

D. Single bunch
In the previous sections we have used MUSIC to simulate coupled bunch instabilities generated by high quality resonant modes under different conditions. Here we want to check the code also for the single bunch beam dynamics, in order to test the matrix transport formalism of the wakefield to study the internal motions of a bunch. As a reference, we use a short range wakefield produced by a broadband resonator with a quality factor Q ¼ 1, a peak impedance R ¼ 6 kΩ, and a resonant frequency of 11.4 GHz. We simulate a small electron machine for synchrotron radiation of total length L 0 ¼ 71.95 m, having a longitudinal damping time τ s ¼ 9 ms, nominal energy of 800 MeV, slippage factor η ¼ −0.0148, harmonic number h ¼ 24, and rf peak voltage of 170 kV [10]. The initial, natural, bunch distribution has a Gaussian shape with rms bunch length of about 0.07 rad, and a relative rms energy spread of 6.7 × 10 −4 .
The theory of the longitudinal single bunch beam dynamics for electron storage rings foresees a beam intensity threshold, below which the wakefield distorts the equilibrium longitudinal distribution without affecting its energy distribution. It is known as potential well distortion and the equilibrium bunch shape can be obtained by solving the Haïssinski integral equation [30], which can be written as ΨðzÞ ð30Þ with λ 0 a normalization constant, σ ε0 the natural rms energy spread, and In order to check the results of MUSIC, we have resorted to the theory and to a single bunch simulation code [18] (which we refer here as SBSC) which uses the classical approach with wakefield, slices and the convolution sum given by Eq. (5). In Fig. 8, on the left side, we show a comparison of the bunch shape given by the Haïssinski equation, MUSIC and SBSC when the bunch is in the potential well distortion regime with a bunch population of N p ¼ 9 × 10 9 . All three methods give the same results. On the right side of the same figure we show the comparison of the rms bunch length (top) and energy spread (bottom) as a function of the number of turns. Also these figures show a very good agreement.
We can observe from the figure that also the synchrotron phase shift due to the bunch losses is well predicted by MUSIC. This quantity can also be evaluated analytically in the case of a Gaussian bunch interacting with a resonator. The results of MUSIC have been checked and they coincide with the analytical predictions if we use either broadband or narrow band resonators.

E. Microwave instability threshold
The potential well distortion theory predicts a bunch length increasing with current and an energy spread constant up to a given threshold, called microwave instability threshold, above which also the energy spread increases. This is exactly what is found by MUSIC. In Fig. 9, we show the rms bunch length (left side) and energy spread (right side) as a function of the bunch population. From the figure we can conclude that the microwave instability threshold is at about N p ¼ 1.2 × 10 10 .
In order to compare these results with the theoretical predictions, we have used Ref. [31], in which the linear mode coupling theory is applied by considering an expansion of the potential well distorted longitudinal distribution in terms of Laguerre polynomials. In particular, in Fig. 11 of Ref. [31] it is shown the microwave instability threshold as a function of ω r σ z =c with a Q ¼ 1 broadband impedance We also investigated in more detail the effect of the microwave instability in terms of mode coupling. At different bunch intensities, we have performed a Fourier analysis of the energy spread as a function of the number of turns. As an example, on the left side of Fig. 10, the energy spread spectrum amplitude of the zero intensity bunch is shown. The first four azimuthal modes are clearly evident. When we increase the bunch population, the peak lines shift and become larger. A contour plot of the spectrum amplitude (the log scale) at different bunch populations is shown on the right side of the figure.
The higher values of amplitude are represented in redblack, the lower in blue. There is a clear shift of the azimuthal modes 3 and 4 towards lower frequencies, some radial modes of m ¼ 4 encounter first the modes of m ¼ 3 at about N p ¼ 10 10 , but the explosion of the spectrum amplitude seems to happen at about N p ¼ 1.2 × 10 10 where modes m ¼ 3 encounter the modes m ¼ 2 and we find some red-black zones. These results do not intend to explain the microwave instability, which is a complicated phenomenon in which many effects combine together, and there is an uncertainty in the calculation of the threshold due to numerical noise, but they show that the transport matrix wakefield approach used by MUSIC is a reliable tool also to study the microwave instability without recurring to slicing of the bunch distribution.
As a final check of this regime, in Fig. 11 we compare MUSIC results with the single bunch simulation code SBSC also above the microwave instability threshold, where the bunch longitudinal and energy distributions oscillate without finding any equilibrium. The comparisons between rms bunch length (left side) and energy spread (right side) as a function of the number of turns show a very good agreement.

F. Kickers, pure inductive, and resistive wall impedance
In the previous sections we have tested MUSIC under different conditions relative to multibunch and single bunch beam dynamics, but in all cases we have used resonator impedances, either with high quality factor or broadband. For long range wakefields, this is what we usually find, since high Q resonant modes are the causes of the coupled bunch instabilities, and for short range wakefields, a broadband impedance can be used as a simple model for a machine impedance [32]. However, the spread of reliable electromagnetic solvers in recent years has permitted one to improve the impedance model of an accelerator, and often a numerical function of the impedance versus frequency is known with a good detail, thus giving a more realistic impedance than a simple broadband model.
In order to avoid the convolution sum, the use of slices, and the problem of the wake potential, in MUSIC we need resonators to transport the wakefield according to the matrix given by Eq. (10). In these cases it is necessary to fit the impedance with a sum of resonant modes. As an example, let us consider the longitudinal impedance of the PS kickers [18], which are the main contributors to the total impedance of the PS. The real and imaginary part of the impedance is shown in Fig. 12 compared with the fit of a single resonator. As can be seen, the approximation is bad. However, by using, for example, 18 resonant modes, the fit, shown in Fig. 13 compared to the initial impedance, is very good.
The fit has been done only to the real part of the impedance, but, as we can see in the right side of the figure, also the imaginary part has been, consequently, well fitted.
The fitting procedure of the impedance as a sum of resonant modes needs here some comments. The fit is strongly nonlinear, and, in general, it is not possible to use directly standard least square methods. Our procedure consisted then in, first, identifying the frequency of the peak impedance (or more than one if necessary) and using it as a resonant frequency of a first mode. Then we have chosen Q and R in order to have a function fitting only this peak. We subtracted the obtained resonant mode to the initial impedance, and repeated the procedure until the differences become negligible. It is important to note that, from the mathematical point of view, also negative values of R can be used for the fit and the transport matrix.
By using in MUSIC the values found for the 18 resonant modes in the potential well distortion regime with the beam and machine parameters used for single bunch simulations, the resulting rms bunch length and energy spread are shown in Fig. 14. Even if the approach uses the wakefield matrix for each mode instead of the wakefield with slices, the results are in very good agreement with SBSC. We observe that the effect of the wakefield here is to reduce the bunch length since it acts as a capacitive impedance due the high bunch cutoff frequency used in the simulation. In another test of MUSIC, we have studied the case of a pure inductive impedance, which could be of interest when the space charge gives a high contribution to the machine impedance. The main problem arising in the simulations of this effect is that the corresponding wakefield is proportional to the derivative of a Dirac delta function, so that the numerical manipulation of Eq. (5) with a code using wakefield and slices is impossible. To overcome this problem, we can observe that, if w ∥ ðzÞ ∝ δ 0 ðzÞ, V wf ðzÞ is proportional to the derivative of the longitudinal distribution, and we can profit of this observation by modifying the classical approach evaluating directly Eq. (5) by means of the derivative of λðzÞ at the center of the slices, as it is done in a modified version of SCSB.
As we have underlined in Sec. III, for simulating wakefield effects, MUSIC needs, as input, a fit of the impedance with a sum of resonant modes. In this case, if we use a single resonator with a resonant frequency much higher than the bunch cutoff frequency, the resistive part of the impedance seen by the bunch can become negligible leaving only the interaction of the bunch with the inductive part. As an example, we have represented on the left side of Fig. 15 a bunch spectrum used for the simulations together with a large broadband resonator ðQ ¼ 0.501Þ in order to model the pure inductive effect. From the figure, we see that the real part of the impedance has an almost negligible interaction with the bunch.
With a pure inductive impedance, the potential well distortion theory foresees a lengthening of the bunch distribution without any phase shift of the center of mass. The new distribution, starting from a Gaussian shape, becomes almost parabolic at its center and it remains symmetric. Indeed the MUSIC simulation results are in very good agreement with the potential well distortion theory and SCSB, as shown on the right side of Fig. 15.
In some cases, in particular for proton machines, the finite resistivity of the vacuum chamber gives an important contribution to the total broadband coupling impedance. The resistive wall impedance, when the skin depth is much smaller than the wall thickness, can be written as [3] ZðωÞ with L the length of the pipe, b its radius, Z 0 the vacuum impedance, and σ c the wall conductivity. This impedance model is valid in the frequency range ω ≫ χc=b with χ ¼ 1=ðZ 0 σ c bÞ. An example of the real and imaginary part of the impedance, with a fit obtained with only four resonant modes, is shown in Fig. 16. Of course, by increasing the number of resonators, the quality of the fit can be improved. However, if we are interested to simulate long bunches, which interact with the low frequency part of the impedance, the fit, shown in the zoomed part of the plots, is totally unsatisfactory.
At low frequency, with still the condition that χc=b ≪ ω, the real and imaginary parts of the impedance depend on the square root of the frequency and have a behavior quite different from that of a resonant mode given by Eq. (18). However, by increasing the number of resonators, it is still possible to obtain a good fit up to a given frequency, as done, for example, in Fig. 17, for which we have used 18 resonators. As for the case of pure inductive impedance, the fit can be satisfactory for long bunches, while it could become critical for shorter ones. In the figure two Gaussian bunch spectra, having rms lengths of 2 and 0.2 m are also shown.
In this situation we can say that 0.2 m represents about the acceptability limit of the fit, as also shown in Fig. 18, where we have compared the wake potentials of the resistive wall impedance [33] with those of the fits for Gaussian bunches of 2 m (left), 0.2 m (center), 0.02 m (right). For the shortest bunch case, the two curves clearly start to deviate. In order to obtain a good accuracy of the fit for higher frequencies, it is necessary to increase the number of resonators. This clearly increases the computing time and, even if, in principle, there is no limit to the number of modes that can be used for the fit, in practice it is unrealistic to run simulations with an excessive number of resonators. A compromise between the running time and the degree of precision of the fit has, under these conditions, to be found.

V. CONCLUSIONS AND OUTLOOK
In this paper we have presented a novel algorithm which exploits an alternative approach to the currently used convolution sum for wakefield simulations, and developed a tracking code, MUSIC, which simulates the longitudinal beam dynamics in a circular accelerator in presence of multibunch with intrabunch motions. Issues arising in the use of convolution sum are related to the necessity of slicing the bunch distribution, of requiring the knowledge of the Green function of an accelerator, and of recording very long wakefield for the study of coupled bunch instability.
The new code needs a fit of the coupling impedance of the accelerator with resonant modes, and the algorithm exploits a matrix formalism to transport the wakefield from one charge to the following one, without resorting to the convolution sum, thus gaining in computing speed and avoiding the bunch slices. It allows one to simulate simultaneously single and multibunch beam dynamics.
We have benchmarked the code against the theories in several cases: for the coupled bunch instability, with and without a frequency domain feedback system to cure it, and for the single bunch beam dynamics, below and above the microwave instability threshold. In this last case we have also resorted to a classical simulation code which requires slices and the wakefield convolution sum.
In all cases the code gives reliable results, showing that the wakefield matrix formalism can be used as an alternative to the classical convolution sum to simulate in an efficient way wakefield effects. In order to give some reference values, in case of a single bunch interacting with a broadband resonator with 10 5 macroparticles, the computing time needed by the wakefield transport matrix algorithm is comparable to that of a classical approach with about 200 slices. In such a situation the advantage of MUSIC is that it renders the use of slices unnecessary, thereby avoiding a verification if the latter introduce any unphysical results.
The same approach of the wakefield transport matrix can be exploited for the development of a code to simulate the transverse beam dynamics with collective effects, since also a transverse resonant impedance has a wakefield similar to that of Eq. (6), which can be transported with a matrix formalism. Moreover, either for the longitudinal and transverse beam dynamics, the method can be applied to simulate collective effects in linear accelerators with the same advantages shown in this paper.
Finally, other kinds of wakefield matrix, different from the resonant mode model used here, could be investigated in order to facilitate the fit of the impedance as a sum of functions such that also standard least square methods can be used to fit a generic accelerator coupling impedance.

ACKNOWLEDGMENTS
We acknowledge many helpful and stimulating discussions with N. Biancacci, H. Damerau, S. Gilardoni, E. Métral, B. Salvant, and G. Sterbini. One of the authors, M. M., would also like to thank the BE/ABP group for the warm hospitality during his stay at CERN. This work was done in the framework of the PS-LIU project.

APPENDIX: WAKEFIELD MATRIX
A resonant mode, in the limit case of β ¼ 1, can be schematized with an electric RLC parallel circuit, driven by a point charge current i b ðtÞ ¼ qδðtÞ, with δ the Dirac delta function, as shown in Fig. 19. With this scheme the wakefield produced by q is the voltage VðtÞ across the capacitance divided by the charge itself.
We want to write the time evolution of the wakefield in matrix form. The time behavior of the voltage can be easily found by solving the differential equation with ω r ¼ 1= ffiffiffiffiffiffi ffi LC p and α ¼ 1=ð2RCÞ. The quality factor of the RLC circuit is Q ¼ ω r =ð2αÞ. The solution of this equation gives free oscillations of the kind 2Q t ½A cosðω n tÞ þ B sinðω n tÞ ðA2Þ with ω n ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ω 2 r − α 2 p . As a second variable we use the current in the inductance. It is related to the voltage by so that it obeys a similar equation, I L ðtÞ ¼ e − ω r 2Q t ½A 1 cosðω n tÞ þ B 1 sinðω n tÞ: ðA4Þ Now we determine the four parameters A, B, A 1 , B 1 . As initial conditions we use Vð0Þ ¼ V 0 and I L ð0Þ ¼ I L0 , which give A ¼ V 0 , and A 1 ¼ I L0 . By considering Eq. (A3) at time t ¼ 0 we get and, by noting that the current coming out from the capacitor must be equal to the sum of the currents in the inductance and resistance, we have B ¼ − Rω r Qω n I L0 − ω r 2Qω n V 0 : We can then rewrite the two equations for the voltage and the current in the inductance depending on the initial conditions as with MðzÞ a 2 × 2 matrix of the kind 2Q z=c cosðω n z=cÞ − ω r 2Qω n sinðω n z=cÞ − Rω r Qω n sinðω n z=cÞ ω r Q Rω n sinðω n z=cÞ cosðω n z=cÞ þ ω r 2Qω n sinðω n z=cÞ ! ðA9Þ for which we have considered z ¼ ct since a second charge, coming after the first one, from the adopted convention, must have a positive z. We recognize in the term (1,1) of the matrix MðzÞ a quantity proportional to the wakefield given by Eq. (6). Let us suppose that the circuit is initially unloaded, and, at the time t ¼ 0, a charge q excites the mode. The capacitor is charged with a voltage V 0 ¼ q=C ¼ qRω r =Q, which becomes the new initial condition V 0 , and the initial current in the inductance I L0 is zero. The initial conditions in the matrix form are then V By using these new initial conditions, we can write the value of the wakefield in any position behind the charge as the first term of the 2 × 1 matrix: