Recent advances in the CERN PS impedance model and instability simulations following the LHC Injectors Upgrade project

Transverse instability growth rates in the CERN Proton Synchrotron (PS) are studied thanks to the recently updated impedance model of the machine. Using this model, macroparticle tracking simulations were performed with a new method well-suited for the slicing of short wakes, which achieves comparable performance to the originally implemented method while reducing the required number of slices by a factor of 5 to 10. Furthermore, dedicated beam-based measurement campaigns were carried out to benchmark the impedance model. Until now, beam dynamics simulations based on this model underestimated instability growth rates at injection energy. Thanks to a recent addition to the impedance model, namely the kicker magnets' connecting cables and their external circuits, the simulated instability growth rates are now comparable to the measured ones even when neglecting the impact of the space charge force. Finally, the space charge force is included in simulations and its impact on the instability growth rate and intra-bunch motion is studied.


Introduction
A horizontal instability regularly occurred in the PS at injection energy on operational LHC-type beams (of longitudinal emittance  l = 2 eV.s) during the gradual beam parameters ramp-up to reach LIU (LHC Injectors Upgrade project) specifications [1].Once the instability was mitigated, by correcting chromaticity, comprehensive instability growth rate measurements were performed with the aim of continuing the benchmark of the PS impedance model after its recent hardware changes, with beam-based measurements [2].It turned out that the current model, when injected into PyHEADTAIL [3] simulations, was not able to reproduce the instability, in particular growth rates were underestimated [4].Nonetheless, it provided valuable insight on the origin of the discrepancy and pointed to a missing impedance source.
-1 -In this contribution, the transverse impedance frequency range responsible for the instability is narrowed down thanks to the measured instability intra-bunch motion.Then, a candidate impedance contribution missing in the model for this frequency range is studied: the connecting cables and external circuits of the kicker magnets.An analytical model based on the transmission line theory is derived to model this impedance contribution.Before performing macroparticle tracking simulations to assess the impact of this new contribution, a novel method to compute the wake kick in macroparticle tracking codes is proposed to speed up simulations.Finally, simulations using the updated model are compared with beam-based measurements.Then, the space charge force is introduced in simulations and its impact discussed.

Measurements of the flat bottom horizontal instability
The instability characterization and the mitigation strategies to tackle it were presented in ref. [4].Following the instability mitigation efforts, numerous instability growth rate measurements were performed, without transverse feedback and using a single bunch.The horizontal chromaticity was set to positive (in particular,  ′ x = 0.6 when the instability was observed for the first time), zero and negative values during the measurements.Simulations failed at reproducing the beam-based measurements as they underestimated the instability growth rates by a factor of 5 to 10 as shown in ref. [4].Subsequently, an investigation started aiming at understanding the origin of this discrepancy.

Horizontal instability characterization
During the measurements campaign, instability growth rates were measured as well as the instability intra-bunch motion.In the PS, the bunch centroid and intra-bunch motion can be observed by means of a wide-band pick-up.This device allows acquiring the longitudinal and transverse motion within a specified time interval spanning between a single bunch length to a revolution period.Observing in parallel the longitudinal and transverse signals allows us to identify the number of nodes in the transverse pattern, as well as a potential asymmetry between the head and the tail.The longitudinal and transverse motions of a bunch at the beginning of the exponential increase of the beam centroid were measured with the wide-band pick-up over 50 acquisitions every 3 turns.Additionally, the instability power spectrum is computed by means of a Fourier Transform using the NEFFINT [6,7] algorithm.Both the intra-bunch motion and power spectrum are displayed in figure 1.Each acquisition spectrum is computed and their average is plotted in white dashed lines in the plots.The instability intra-bunch motion is characterized by a mode-0-like envelope (no node) featuring an asymmetry between head and tail of the bunch.Its resulting power spectrum exhibits a main peak around 3 MHz.
One hypothesis is that the peak observed in the power spectrum corresponds to resonant modes in the impedance spectrum.As the power spectrum is the frequency domain equivalent of the transverse oscillations, it allows quantifying the frequencies of interest in the impedance spectrum.In the following, we will investigate this possibility.

Identification of the missing impedance source
According to theory [8], the effects of transverse impedance on beam dynamics, such as tune shift or instability growth rates, result from the sum along the betatronic lines of the product of the impedance and the unstable mode spectrum.Thus, by overlapping the unstable mode power spectrum from  The dominant contribution in the impedance spectrum in the power spectrum frequency range is the vacuum chamber impedance.All other contributions, displayed in the 'Others' category, remain small for frequencies lower than 10 MHz.Noticing that the frequency range of the unstable mode power spectrum (0 MHz to 10 MHz) matches the one of the PS kicker magnet BFA09S connecting cables impedance [9], the kicker magnets connection cables and their external circuits are a candidate to explain the impedance missing in the PS impedance model.

Impedance of the PS kicker magnets connection cables and their external circuits
The electromagnetic fields induced by a bunch of charged particles traveling through a kicker magnet can couple with the kicker core and its external circuit through the busbar.The external circuit of a kicker magnet includes the cables, filters, some ferrite pieces and the termination load.The coupling with the external circuit leads to an additional contribution to the beam coupling impedance of the kicker magnet.This contribution contains a longitudinal and horizontal dipolar component.As explained below, these two components' impedance originates from the self-inductance dependency with the horizontal position of the bunch in the kicker magnet, thus leading to zero quadrupolar and vertical impedances.Due to the attenuation properties of the cables, the impedance caused by the external circuit is limited to the MHz range.The external circuit impedance can be calculated analytically by approximating the magnet as an ideal transformer and using the resulting circuit model of the magnet and connecting cables as illustrated in figure 3.
Schematic of the circuit model of a kicker magnet as an ideal transformer [10].
In the following, the impact of the termination load and inductance of a kicker on its cable impedance is studied on the BFA09S and BFA09P kicker magnets, based on the formalism developed by C. Zannini [10].Their resulting impedances are discussed and the impact of their termination is assessed.

Transmission line theory
The impedance of the cables is calculated from the kicker inductance , self-inductance , cable impedance  TL , and  g the external circuit impedance.The latter is expressed using the transmission line theory [14].Electromagnetic waves traveling through a transmission line can be described by wave equations, which depend on the characteristic impedance and propagation constant of the line.
The characteristic impedance ( C ) depends on the physical properties of the transmission line, such as its geometry and the electromagnetic properties of the materials it is made of.Besides, it also accounts for a reflected wave, thus making it different from the line impedance.
The propagation is dictated by the complex propagation constant , the real part of which characterizes the wave attenuation along the transmission line, while its imaginary part governs the phase shift - is strictly imaginary for a lossless line.We will write it as: where  is the attenuation constant in Nepers/m and  is the phase constant in radians/m, both quantities being frequency-dependent.
-4 - The attenuation constant represents the rate at which the signal power decreases as it travels along the transmission line.In a coaxial cable, the attenuation constant follows a logarithmic behaviour with frequency.The phase constant represents the phase shift of the signal as it propagates through the cable.The phase constant can be expressed using the phase velocity inside the transmission line: with  p ≡  √  0  0 / √  the phase velocity,  the speed of light in vacuum,  0 the vacuum permittivity,  0 the vacuum permeability,  r = / 0 the relative permittivity of the cable,  its permeability, and  the angular frequency.In eq.(3.2) we have assumed the magnetic permeability of the cable is the same as the one in vacuum ( =  0 ).
The impedance arising from the connecting cables and their termination can be calculated analytically by resorting to the Telegrapher's equations [14].A formulation for the transmission line impedance based on the telegrapher's equations can be written as: where  L is the load or termination impedance of the line, and  is the position along the line.By replacing  with  the length of the line, the total impedance of the line can be calculated (i.e. g in figure 3).The formulation used is general and allows losses along the line.To go to the lossless case the terms tanh( ) can be replaced by  tan(): The lossless line impedance corresponds to the one derived in ref. [10].Based on the circuit model of the kicker shown in figure 3, the generalized impedance arising from the coupling of the external circuit with the magnet [12] can be expressed as: where  is the horizontal position at which the impedance is evaluated and  0 the horizontal position of the bunch passing through the kicker magnet with respect to its centre line.The impedance depends on the magnet inductance  and self-inductance .Estimations of the magnet inductance and self-inductance can be found in ref. [11] for H and C-shape magnets: ) where  and  are respectively the horizontal and vertical half apertures of the magnet winding.
The longitudinal impedance can be calculated from eq. (3.5) by plugging the inductance and selfinductance expressions derived in eq.(3.6) and eq.(3.7) respectively, assuming the bunch passes through the centre of the magnet: g  + g (C-shape magnet). (3.8) -5 -It can be observed that only a C-shape magnet demonstrates a nonzero longitudinal impedance.Using the Panofsky-Wenzel theorem [15], the horizontal impedance can be calculated from eq. (3.8) as follows: The resulting horizontal impedance is independent of the magnet's geometry and yields the same formula for H-shape and C-shape magnets.The formalism is applied to the PS BFA09 kicker magnets in the following section.

Calculation of BFA09 Kicker Magnet Connecting Cables Impedance
BFA are lumped C-shape magnets designed with six conducting cables for the pedestal module (suffix P) and two for the staircase module (suffix S).They were used in the discontinued Continuous Extraction scheme.The kickers BFA09P and BFA09S modules are identical but the dimensions and shapes of their conducting cables as well as their termination load differ.The BFA09 module is characterized by a semi-horizontal aperture  of 79 mm and a semi-vertical aperture  of 26.25 mm.After the second long shutdown of the CERN complex (LS2, 2019-2021), the BFA21S and BFA21P were removed from the ring while the BFA09S and BFA09P were left in place.The BFA09S cable termination is open while the BFA09P is matched to the characteristic impedance of the cables.The reason for the different terminations is due to the BFA09P still being in operation whereas the BFA09S is not powered anymore.It is considered to remove the BFA09S from the PS during the next Long Shutdown (LS3, 2026-2029).The circuit model of the BFA09S can be seen in figure 4, the difference with the current BFA09S is the absence of a termination load.Both kickers are connected to the same kind of RG220 ( C = 50 Ω) cables whose technical specifications were obtained from the manufacturer.Fitting the attenuation values at different frequencies from table 1 yields  =  log (  ) where A = 2.1 × 10 −7 Np/m.The phase constant is obtained using the relative permittivity  r of the cable, which is listed in table 1.
The ratio of the mutual inductance at  = 0 over the inductance is equal to 1/2 based on eqs.(3.6) and (3.7).
The resulting longitudinal and dipolar impedances for BFA09S and BFA09P cables are shown in figures 5 and 6.The main differences between their impedances originate from their respective -6 - In addition, the attenuation in the connecting cables is considered low compared to its phase, leading to  ≪ .This leads to very large coth( ) when   = 0 mod .This occurs at multiple of 542 kHz which is close to the PS revolution frequency at injection ( 452 kHz).This particular spacing reveals itself to be the most critical one regarding beam stability.Depending on the machine tune, the bunch spectrum betatronic lines are expected to overlap with some of the resonance peaks and may cause beam instabilities.The most important overlap between the betatronic lines and the resonance peaks is obtained for the horizontal tune  x = 6.5.In practice, the horizontal tune at injection is bounded between 6.1 and 6.24, thus minimizing the potential impact of the BFA09S cables.

Contribution of PS Kicker Magnets Connecting Cables and External Circuits to the Impedance Model
The PS contains seven additional kickers, each of them with a specific external circuit and connecting cables.Some of the kickers' external circuits include ferrite and RC filters, strategically placed between the magnet and its transmission cables to effectively dampen undesirable signals.In any case, as they are all characterized by a matched termination, their impedance is comparable to the one of the BFA09P and can be computed in the same way.The overall PS kickers cables impedance can be computed and compared with the rest of the impedance model -this is shown in figure 7. The impedance contribution is characterized by a broadband behaviour resulting from the sum of all the remaining kicker cables with matched termination but different external circuits.In addition, its maximum amplitude frequency coincides with the 3 MHz peak of the power spectrum.Furthermore, the impedance is responsible for close to two thirds of the real impedance and one third of the imaginary impedance in the 100 kHz -5 MHz frequency range, surpassing the vacuum chamber -8 -

Macroparticle tracking simulations of a rapidly oscillating wake
At injection, a PS bunch length typically varies between 140 and 200 ns whereas the vacuum chamber wake function varies considerably over less than 0.5 ns.One of the shortcomings of the current wake kick implementation in the macroparticle code PyHEADTAIL [3], when dealing with such different time scales, is the very large number of longitudinal slices within a bunch, which leads to unmanageable simulation times.In the following we address this limitation, before assessing the impact of the kickers' cables additional impedance contribution.

Wake kick calculation in tracking codes
The effect of wakefields on beam dynamics is expressed by the Lorentz force of the resonating electromagnetic fields produced by a bunch of charged particles after their passage in an accelerator element.The integral of the Lorentz force over the element length yields the longitudinal and transverse momentum kicks.For example, the momentum kick caused by a dipolar force is calculated from the convolution of the dipolar wake function with the bunch longitudinal distribution as follows (here in horizontal): where ⟨()⟩ is the average horizontal position of the particles along the bunch,  the longitudinal position of the particle experiencing the force,  the elementary charge,  p the mass of the proton,  the relativistic beta,  the Lorentz factor, () the dipolar (horizontal) wake function and () the beam longitudinal distribution.Macroparticle tracking codes rely on a discrete representation of a bunch instead of a continuous distribution.On top of that, the bunch is typically sliced longitudinally to compute the wake convolution, as we will see below.

Stepwise wake kick computation
A macroparticle is defined by an ensemble of particles sharing their macroscopic representation in terms of overall charge, energy and position.Once a beam is expressed as a body of macroparticles, it is longitudinally sliced along an equidistant grid and each macroparticle is placed in the slice corresponding to its longitudinal position.As a result, for dipolar wakes, a unique kick can be applied to each slice.Rewriting eq.(4.1) with macroparticles and slices, yields: where  slices is the number of slices considered,  n expresses the number of macroparticles in the slice ,  n the horizontal centre of mass position of the slice  and  n−i the wake function for  corresponding to the distance between the centre of the slices  and  (the latter being the one for which the kick is evaluated).Note that here, the wake from behind  is neglected.
-9 -A slice located in the tail of the bunch experiences a wake kick summing the effect of all the slices from the head of the bunch to its location, excluding the slice where the wake kick is computed.The wake kick is dependent on the slice transverse center of mass, and number of macroparticles.Nevertheless, using this approach requires expressing the wake function by a single value in each slice.In PyHEADTAIL, currently, this value is defined as the wake function evaluated at the centre of a slice.Two limitations arise from this solution.First, the slicing must be chosen sufficiently fine to evaluate the wake function behaviour correctly, the key point being that the granularity of the approach does not go below the slice level.Then, an adequate number of macroparticles per slice needs to be simulated to prevent statistical noise effects from undersampling, as warned in ref. [16].A general rule of thumb to ensure statistical-noise-effect-free simulations is to choose around 700 macroparticles per slice.
When simulating the PS with macroparticle tracking codes, the most critical case is at injection energy, where the indirect space charge impedance dominates and translates into an extremely narrowly peaked wake function located around  = 0.After extensive studies aiming at finding the required number of slices to ensure adequate sampling of the indirect space charge wake, it was shown that 8000 slices are needed.To follow the recommendation of 700 macroparticles per slice, a simulation calls for 5.6 × 10 6 macroparticles.As a result, a single simulation lasts more than two weeks on the CERN cluster service.Such a time span becomes unmanageable when numerous scans of various beam parameters need to be carried out.For this reason, efforts have been made to rethink the way the wake function is evaluated in each slice.

Novel integrated wakefield technique
To address the above limitation, the integration over the slices is rethought.The main challenge lies in finding a strategy to account for the full information of the wake function inside a slice instead of evaluating it at the center position of the slice.An elegant method, initially proposed by G. Iadarola, is to make use of the integrated wake function over each slice.Equation 4.1 can be rewritten as a sum of integrals over each slice.Assuming  n  n can be considered constant within a slice, it can be taken out of the integral, leading to: where Δ is the slice width.The wake function integral over the slice width can be regarded as a local kick factor normalized by the slice width.This assumption holds provided that the longitudinal and transverse oscillation patterns of the bunch are longer than the slice width, i.e. that the number of slices is significantly larger than the number of oscillations.Now, thanks to this assumption, the only remaining term in the integral is the wake function.The integral can be calculated once at the beginning of a simulation and replace the previous interpolated value of the wake function at the centre of a slice.Even if a coarse slicing is chosen, the average effect of the wake function oscillations inside a given slice will be reflected in the integrated value.

Comparison between the two methods
To benchmark both methods, they are applied to the PS at injection energy, including only the vacuum chamber wake, which is both a dominant contribution and the one requiring most slices to be fully -10 -resolved.The example configuration is depicted in figure 8 where the bunch distribution is compared to the vacuum chamber wake function, highlighting the need for very thin slices in the stepwise method as the wake varies sharply within less than 0.5 ns.The wake kick calculated for an increasing number of slices with the stepwise and integrated methods is shown in figure 9.The aim is to study when the wake kick amplitude versus  converges with the number of slices.On one hand, the stepwise method reaches convergence only above 5000 slices, with a maximum amplitude close to 3 × 10 −8 .On the other hand, the integrated method is almost converged for 50 slices and fully converged for 500 slices, where the ripples around the head and the peak of the wake kick vanish.Similarly to the stepwise method, the wake kick peak amplitude converges around 3 × 10 −8 .In this example, a tenfold reduction in the required number of slices is achieved.To use the integrated method in simulations, the wake function needs to be integrated first with respect to the chosen slicing, as discussed above.Therefore, additional computation time for that integration is needed at the beginning of a simulation, contrary to the stepwise method.Nevertheless, -11 -after a few turns this integration time (which is done only once) is largely compensated by the time gained from the smaller number of slices and number of macroparticles required, as confirmed by figure 10.

Computation time [s]
Stepwise method Integrated method During the tracking, the wake kick computation time is equivalent for both methods.The integrated method is compatible with single and multi-turns wake.Its validity is benchmarked against the stepwise method for the single and multi-turns wake cases in the following section.

Comparison of the two methods for impedance-induced beam observables
In the previous section, the integrated method was benchmarked against the stepwise one by looking only at the single-turn wake kick.To fully validate the new method, the growth rate of the instability triggered by the PS vacuum chamber wake is studied for several intensities and number of slices, as displayed in figure 11. -12 -We observe that the required number of slices to reach convergence agrees with the results from figure 9.This observation opens up the general possibility of running a preliminary convergence check by studying the wake kick on a single turn, as done above, instead of running a full simulation.This drastically reduces the number of needed simulations and the overall time required by the convergence studies.
Previous simulations were run with a wake function that is considered fully decayed after one revolution in the machine -this is well-suited only for rapidly decaying wakes.We now turn to multi-turn wakes, i.e. considering that the wake function persists after one or several revolutions.For instance, the PS vacuum chamber wake requires 10 turns to decay to a point where its impact on the beam dynamics becomes negligible.
As visible in figure 12, the effect of such wake accumulation is identical for both methods, as long as a converged number of slices is used -both methods yield similar instability growth rates and impedance-induced tune shifts.The agreement between both methods for the tune shift confirms the usefulness and validity of the integrated method in configurations where broad-band, inductive impedances, such as the indirect space charge at low energy, are dominant, as is the case here.Indeed, the integrated method captures the tune shift from the indirect space charge impedance, similarly to the stepwise method with significantly fewer slices.
The integrated method is now available in PyHEADTAIL and also compatible with Xsuite [17], both with CPUs and GPUs for each code.The integrated method has been routinely employed for PS simulations, effectively reducing the computation time from a week to 8 h.Using the integrated method, simulations of the flat bottom instability can be performed and compared with measurements.

Beam-based measurements and simulations of the PS flat bottom horizontal instability 5.1 Comparison in the absence of direct space charge
The observable of interest used to compare measurements and simulations is the instability growth rate.It is obtained by fitting an exponential function to the evolution of the transverse beam position over -13 -the number of turns.To do so, the envelope of the transverse oscillations is retrieved by calculating the Hilbert transform [18] of the signal.Then, a fitting procedure is carried out by applying a modified least-square method with an exponential function [19].An example of turn-by-turn data of a measured instability and the fitted growth rate can be found in figure 13.It can be noted that, due to the presence of space charge and non-linearities in the accelerator, the instability growth rate evolves with the number of turns.The fitted instability growth rate reaches its maximum at the beginning of the instability when its amplitude remains small.In some cases, an instability can be damped through Landau damping before causing any beam intensity losses.In the following, the start of the fitting window is set to a doubling of the signal envelope baseline value (i.e. when the beam is stable, before the instability starts), and its end to a twofold increase in the signal amplitude.
Simulations were performed with PyHEADTAIL for the three chromaticities for which instability growth rates were measured in the PS.All the relevant measurement parameters and PyHEADTAIL settings can be found in appendix A, and the PS impedances and wake functions in appendix B. To start with, to confirm the role of the kicker connecting cables and external circuits impedance, simulations were performed with the PS impedance model but excluding this contribution, then including it back.Results are shown in figure 14.
An agreement is obtained between measurements and simulations only when the kicker connecting cables and external circuits impedance is included.Without them, the instability growth rate (which essentially results from the vacuum chamber and beam instrumentation impedances) is underestimated by one order of magnitude.Overall, the horizontal instability is the joint consequence of these three contributions.
Using the machine and beam parameters related to the measured intra-bunch motion presented in figure 1, the intra-bunch motion was also simulated and is shown in figure 15.
-14 - First, a peak around 3.4 MHz is visible in the simulated intra-bunch power spectrum, as was the case in the measured one.At the same time, the simulated intra-bunch motion exhibits no head-tail asymmetry, which were observed in measurements.One hypothesis is that the asymmetry could be caused by the effect of direct space charge [20], which is currently not included in the simulations.

Comparison including direct space charge
The space charge force is characterized by its dependency on the beam energy, intensity, transverse and longitudinal sizes.However, the PS transverse emittance increases with the beam intensity due to space-charge-driven emittance blow-up in the PSB.Therefore, the transverse emittance should be known over an adequate range of beam intensity, as seen in figure 16, to realistically portray its effect in simulations.The emittance was measured in the horizontal plane and the beam is expected to exhibit equal horizontal and vertical emittances.It can be noted that the longitudinal emittance remains constant with beam intensity.
-15 - Then, the direct space charge effects and octupolar imperfections introduce some amplitude detuning according to the variation of the beam size along the accelerator.Therefore, the ring needs to be split into an adequate number of segments to correctly depict the amplitude detuning.The PS is composed of 100 main units (combined function magnets), each composed of a focusing and defocusing half-unit.Thus, the beta functions and dispersion sampling needs to be finer than the length of a half-unit.For this reason, the ring is split into 240 segments.It was found in ref. [19] that, for this number of segments, the beta functions and dispersion are realistically reproduced with respect to MAD-X.
The space charge force is applied after each segment, using a 3D Particle-In-Cell (PIC) solver [21] and taking into account the transverse beam size calculated from the beta functions and dispersion.The PIC solver mesh parameters are summarized in table 2. Apertures were added around the PIC grid to avoid unphysical behaviour when particles go out of the grid and are simulated without the space charge forces.As a result, particles located in the tails of the bunch distribution reaching the apertures are considered lost.As the particles located in the tails are closely tied to the beam stability [22], it is of paramount importance to simulate properly their behaviour.
Each of the PIC solver parameters was scanned and its convergence was ensured.The transverse grid extent was set at ±8 x,y to allow for some space-charge-driven emittance blow-up before the start of the instability, without causing any beam intensity losses due to the PIC grid apertures.In addition to the direct space charge force, the octupolar non-linearities are also accounted for in the simulations.Their impact is quantified by the anharmonicity coefficients [23]  The octupolar non-linearities act on the beam after each segment through a transverse momentum change depending on the particles' transverse actions and phase advance at this location.Considering the PS optics for LHC beam production [24], the detuning coefficients are set to  xx = 25.6 m −1 , Contrary to simulations without space charge and as pointed out in figure 13, the instability growth rate has to be obtained from the earliest stage of the instability, rather than from the full signal.Indeed, as the amplitude of the bunch transverse oscillations increases, so do transverse emittance and amplitude detuning, hence stabilizing the growth throughout the simulations.Conversely, subsequent instabilities may arise from the perturbed bunch distribution.The measured turn-by-turn signal baseline cannot be reproduced in simulations as it is mainly dictated by the acquisition equipment resolution when the bunch transverse oscillations remain small.Instead, a different criterion is used in simulations.The exponential fit is performed from the start of the exponential increase of the transverse oscillations until the transverse emittance has increased by 10%.It is assumed that, in this regime, the first observed instability in the simulation can be studied while its impact did not yet perturb significantly the bunch distribution.
The simulations accounting for space charge are compared against previous simulations and measurements in figure 17.The agreement between measurements and simulations with space charge varies depending on the chromaticity, worsening as it decreases.Nonetheless, a fair agreement is obtained for intensities under 1 × 10 12 p.Above this value, the addition of space charge in simulations seems to lead to an overestimation of the instability growth rate.In any case, the simulated growth rates remain comparable to the measured ones within 30%.A potential explanation for this discrepancy could be related to the fitting procedure and the different fitting windows used in both cases, leading to different fitted growth rates, as highlighted in figure 13.Fitting window boundaries applicable to both measurements and simulations are required to realistically compare -17 -instability growth rates.Another possible criterion for the end of the fitting window is to stop when the second derivative of the growth rate turns from positive to negative, indicating an inflection point and that its maximum value has been reached.This fitting procedure is inspired by the FITX library [25] On top of the instability growth rate, the intra-bunch motion can also be retrieved from simulations including space charge.The simulated intra-bunch motion is recorded over 128 turns when the transverse emittance reaches an increase of 10%.
The resulting measured and simulated bunch profiles and horizontal intra-bunch motions for the three chromaticities studied previously, are presented in figures 18 and 19.In each case, the tail of the distribution is located on the right-hand side of the plots and the head on the left-hand side.-18 -

JINST 19 P04018
The measured and simulated intra-bunch motions agree only for  ′ x = 0.6, where the best agreement is also found with the instability growth rate.For the remaining chromaticities, measurements demonstrate transverse oscillations with a fine structure characterized by fast oscillations towards the tail of the bunch.While the head/tail asymmetry is reproduced by simulations, the various nodes along the bunch remain missing.
The outcome of including the space charge force in simulations leads to an overestimation of the instability growth rate and a poor agreement with the measured transverse intra-bunch motion.On one hand, the instability growth rate discrepancy between measurements and simulations could be attributable to the different fitting windows used.On the other hand, the different intra-bunch motions hint either at a missing or incorrectly modeled mechanism in the simulations or a missing impedance in the model.
In PyHEADTAIL, the impact of the optics is included by supplying the beta functions and dispersion at each segment's location.By doing so, the impact on tracking of the generalized gradient expansion of the accelerator elements is neglected, thus preventing the presence of resonances or effects from non-linear optics on beam dynamics (aside from amplitude detuning, which is included).Furthermore, non-linear chromaticity was not experimentally measured nor accounted for in the simulations.Yet, second-order chromaticity has been proven to impact instabilities' intra-bunch motion [26] and growth rate [19,26].Consequently, the missing mechanisms may influence the beam dynamics to the point of impacting the validity of the space-charge force modeling.Hence, further studies could include simulations with BimBim [27] and Xsuite to validate the results obtained with PyHEADTAIL, followed by the introduction of non-linear chromaticity in either PyHEADTAIL or Xsuite.

Conclusion
The measured instability power spectrum proved to be a useful tool in narrowing down the frequency range responsible for an unexpected instability observed at low energy in the PS.An analytical model of the PS kicker magnets connecting cables and external circuits was developed and found to generate a significant source of impedance.Macroparticle tracking simulations were performed using a novel wake kick calculation method, which allowed to significantly reduce the simulation time.The impedance contribution from the kickers has a significant impact on beam dynamics and its addition to the PS impedance model allows simulations to reproduce the growth rate measurements with quite good accuracy without space charge.When space charge is introduced in simulations, the instability growth rates are overestimated and the intra-bunch motion is only comparable to the one measured with positive chromaticity.

Figure 1 .
Figure 1.Measured instability transverse oscillations (left) and transverse power spectrum (right) over 50 acquisitions.Horizontal chromaticity ( ′ x = 0.6) is set to the value measured at the occasion of the first observation of the instability in the machine.The average transverse power spectrum over all the acquisitions is represented by a dashed white line.

Figure 2 .
Figure 2. Overlap of the power spectrum measured for  x = 6.225 and  ′ x = 0.6 (black line, from figure 1) with the impedance spectrum.The betatronic lines, spaced by  0 , are represented by dashed black lines.

Figure 7 .
Figure 7. Overlap of the measured power spectrum (black line, from figure 1) with the updated impedance spectrum for  x = 6.225 and  ′ x = 0.6.The betatronic lines, spaced by  0 , are represented by dashed black lines.

Figure 8 .
Figure 8. Longitudinal distribution of a PS bunch at injection (top) compared to the vacuum chamber wake (bottom).

Figure 9 .
Figure 9. Wake kick comparison between the stepwise (left) and the integrated methods (right).

Figure 11 .
Figure 11.Instability growth rate versus beam intensity for a various number of slices, for  ′ x = 1, with the stepwise (left) and the integrated (right) methods.

12 Figure 12 .
Figure12.Instability growth rate and tune shift versus beam intensity with a single turn ( AT = 1) or twelve turns ( AT = 12) of accumulated wake, for both wake kick computation methods.

Figure 13 .
Figure 13.Example turn-by-turn data of a measured instability (Intensity = 5.9 × 10 11 p,  ′ x = 0.6) with its resulting growth rate fit (left) and the evolution of the fitted growth rate with the number of turns (right).

Figure 16 .
Figure 16.Emittance measurements against beam intensity at injection energy (flat bottom).

Table 1 .
[9]09S and BFA09P cable specifications[9].Type  [mm]  [mm]  [µH]  cable [m] For instance, BFA09P cables are matched and demonstrate low longitudinal and transverse impedances around the MHz range, before vanishing quickly for higher frequencies.The peak impedance of BFA09P is also twenty times smaller than BFA09S whose cables are left open.In addition to a larger peak impedance, leaving the cable termination open also introduces a significantly slower decay of the longitudinal and transverse impedances with frequency.The sharp resonances in the impedance spectra can be explained from the asymptotic expression of eq.(3.3) when L goes to infinity (open termination).The equation then becomes:

Table 2 .
Parameters for the PyHEADTAIL simulations with space charge.