Nested Head-Tail Vlasov Solver

Nested Head-Tail (NHT) is a Mathematica-based Vlasov solver for transverse oscillations in multi-bunch beams. It takes into account azimuthal, radial, coupled-bunch and beam-beam degrees of freedom, single- and inter-bunch dipole wakes, an arbitrary damper, beam-beam effects and Landau damping.


INTRODUCTION
Collective instabilities impose a serious limitation on beam intensity in circular accelerators, storage rings and colliders. That is why quantitative description of the instabilities is required for understanding of both existing and future machines. Especially dangerous are transverse instabilities, normally resulting in a beam loss. In principle, there are two approaches for their modelling. The older approach is based on the Vlasov equation, a dynamic equation for a phase space density of collisionless plasma; the equation was suggested in 1938. In fact, this equation was first written not for plasma, but for gravitating cosmic objects by J. Jeans in his publications of years 1913-1915, see a historical article of M. Henon [1]. So far, solutions of the Vlasov equation were limited by rather simple cases (see e. g. Ref [2]), insufficient to tell about complicated reality of multibunch beams with feedbacks, beam-beam effects and octupoles. That is why the second approach to beam stability problems, macroparticle simulations, attracted more and more attention, driven by continuing burst of computational powers. For colliders, such codes as HEADTAIL [3] and BeamBeam3D [4] are known and used for more than 10 years. Attractiveness of the macroparticle tracking programs is related to their similarity to real beams; they appear to be as close to reality as possible, allowing rather straightforward introduction of all the factors influencing beam stability. However, an attempt of these direct simulations of reality has its drawback: convergence typically requires a big number of macroparticles per bunch, at the order ~10 6 or so. For thousands bunches per beam in the collider, it makes a required study of multi-dimensional space of parameters so far impossible by means of the macroparticle trackingeven with the help of parallel computations by powerful supercomputers.
This limitation of the tracking methods brings us back to the Vlasov equation with a motivation to develop more sophisticated methods of its solution, where all important factors of reality would be properly taken into account. Nested Head-Tail (NHT) suggests an attempt in this direction [5]. Its name points to its primary idea: solutions of the Vlasov equation are sought as expansions over conventional head-tail functions, defined at a set of nested rings in the longitudinal phase space. Since the impedance, feedback, and coherent beam-beam are described by linear response functions, their analysis is reduced to a standard eigensystem problem that is solved instantaneously, as soon as the related matrices are defined. After this is done, growth rates and tune shifts of all potentially dangerous modes are known. However, with any feedback, that pure linear system is still unstable: its stabilization requires some Landau damping, an intrinsic self-stabilization mechanism caused by nonlinearity of single-particle motion. In general, these anharmonicities lead to very complicated equations (see e. g. Eq.(6.179) in Ref. [2]). However, for many practical cases, the anharmonicities may be treated as perturbations of the linear system. When Landau damping moves the coherent tune shifts not by much, it can be found as a perturbation. That is how Landau damping is treated in this paper, allowing finding thresholds of the instabilities, with both octupoles and beam-beam nonlinearities taken into account. Fortunately, for many practical cases this perturbative approach to Landau damping is justified. For pure transverse nonlinearities it leads to well-known dispersion relations, see Ref. [6] and references therein. Otherwise, more general form of the dispersion relation has to be applied, see a section "Landau Damping". Numerous examples of NHT applications for LHC are shown in the last two sections.

BASIS FUNCTIONS
Bunch longitudinal distribution inside a linear potential well can be represented by a sequence of concentric rings, or air-bags, as it was suggested by Ref. [7]. It appears to be optimal to keep all the rings equally populated, so the phase space density is reflected by variable distances between them. To do that, the phase space has to be divided onto certain number of equally populated areas; then a weighted-average radius has to be found for each area. The number of rings has to reflect the wake properties. Unless the wake is extremely microwave, the default NHT choice of 5 r n  shown in Fig. 1 should work reasonably well. On these concentric rings, a conventional set of head-tail harmonics l  is defined as a basis for phase space density perturbations (see e.g. Ref. [2], Eq.(6.175)): Here  is the synchrotron phase, l is azimuthal head-tail

SINGLE BUNCH EIGENSYSTEM
Any dipole perturbation of the bunch distribution can be expanded over this nested head-tail basis. Following Ref.
[2], assuming time dependence 0 exp( ) Here, X is the perturbation vector with components along the basis functions l  , matrix Ŝ reflects harmonic oscillations inside the RF potential well, matrix Ẑ reflects a single-turn transverse impedance, g is a damper gain, and F is a flat-wake matrix: Here () Z  is the transverse impedance, flat-wake matrix F describes any dynamic response whose variation over the bunch length can be neglected [8]. Note that flat wakes do not satisfy causality condition of conventional wakes, since they are not related to an instantaneous action of the same-beam particles. Examples of flat wakes are those of dampers whose bandwidth is much smaller than an inverse bunch length, inter-bunch and multi-turn wakes, beam-beam responses for beta-function exceeding the bunch length.
Taken with the damper gain g, the term ĝ F describes the response of a flat damper, which sees only centroids and whose kicks are bunch-flat. The term "flat" as it is applied here to the damper assumes space or time flatness, so it is not to be confused with the idea of a high bandwidth damper, whose response is flat in a frequency domain. To avoid this confusion, the term "flat" is applied below for the time or space domain only. Fourier-image of a flat wake is delta-function, so the expression for the flat-wake matrix F in Eq. (3) follows from the impedance matrix Ẑ after a substitution ( ) ( ) , defining the gain g as a damping rate in units of the synchrotron frequency.
In case the damper is so broadband that its response is not flat over the bunch length, its actual non-flat wake function has to be used, making its matrix different fromF .

MULTI-BUNCH EIGENSYSTEM
Inter-bunch wake functions () Ws are normally flat thanks to suppression of high-frequency cavity modes. If so, inter-bunch terms in the equation of motion can be described by means of the flat-wake matrixF , used in Eq. (3) for the flat damper. Inter-bunch wake fields can be conventionally summarized for equidistant bunches whose oscillation amplitudes differ only by a phase factor and M is the number of bunches. Inter-bunch interaction contributes its own term to the originally single-bunch Eq.  (4) These equations reduce the multi-bunch problem to a set of single-bunch ones since every inter-bunch mode  is treated separately by Eq. (4), where the sought-for eigenvector X has the same structure as for the singlebunch problem of Eq. (2). The multi-bunch problem of Eq. (4) assumes that the damping rate g is identical for all the inter-bunch modes, being the same as for the single-bunch case; in other words, it assumes that the damper is a bunch-by-bunch one. In case it is not so, the damping rate becomes a function of the inter-bunch mode number  . This function can be found in a way similar to the wake interbunch coefficients W  . In the time domain, the damper response can be described by means of its wake function () Gs, associated with the gain frequency profile G  by means of the Fourier transform: assuming this wake function to be even (remember that flat wakes are not causal). With that wake, the inter-bunch mode coefficients G  analogous to the wake coefficients W  can be found: In Eq. (4), the mode-independent gain g has to be changed on a mode-dependent one, gg   , and the latter can be expressed as where 0 g stays for the damping rate of the inter-bunch mode 0   .
In case the bunches are not all equidistant, similar summation of the inter-bunch wakes can be done using the "train theorem" [8]. This procedure requires two steps.
First, an eigensystem { , } p R of the total inter-bunch wake (sum of the conventional and damper wakes) has to found: Second, for the examined inter-bunch mode, its eigenvalue p  has to be used in Eq. (4) by means of a According to the "damper theorem" [8], for a sufficiently high damper gain, the inter-bunch (and beam-beam) collective interaction can be neglected, thus reducing the multi-bunch problems (4), (10), (13), (14) to the singlebunch case (2), provided the inter-bunch and beam-beam wakes are flat.

BEAM-BEAM RESPONSE
For colliders, there is one more source which complicates beam dynamics: beam-beam effects. There are two different aspects in the beam-beam interaction: coherent ("strong-strong") and incoherent ("weakstrong") ones. The coherent beam-beam effect is associated with a collective response of one beam to another. This response works as a specific wake function, coupling the two beams. In case the bunch length is much smaller than the beta-function in the collision point, this interaction can be treated as flat in the sense of this paper. Flat kicks of oncoming bunches are constant over the bunch length, being determined by centroid offsets and being independent of all other details of intra-bunch oscillations. The incoherent aspect of beam-beam interactions causes additional anharmonicities of singleparticle motion, thus contributing to Landau damping. In this section, only the coherent part is considered, while the incoherent one is left for a chapter below where all optics nonlinearities are treated together in a framework of the dispersion relation.
Inclusion of the coherent beam-beam interaction doubles dimension of the problem. Each beam can be described by the same Eq. (4), where an additional beambeam coupling term has to be added in the right-hand side.
A simplest case of beam-beam response is one of a single flat collision in a single interaction region (IR). In this situation, Eq. (4) is modified as Here 1,2 X are perturbation vectors for the two beams, and  is a beam-beam tune shift. The first beam-beam term in the right-hand side, for non-round beams see Ref. [9]. A factor of s Q enters in Eqs. (11), (12) due to a convention of this paper to measure all tune shifts in units of the synchrotron tune. For more than one collision per IR, the beam-beam term has to be modified, taking into account that the interbunch phases of the oncoming beam are not identical; they vary according to the considered coupled-bunch mode  . Note that for equidistant bunches, beam-beam collisions do not break the inter-bunch mode structure exp(2 / ) i k M   since that is just a consequence of the translational invariance. Summation over 21 K  mirror-symmetric kicks results in a following modification of Eq. (10): where  is the total linear part of the incoherent beambeam tune shift. One more complication of the beam-beam coupling appears when there is more than one IR. To avoid unnecessarily cumbersome expressions, let's assume that there are only two IRs. In case the IRs are completely identical, and the betatron phase advances between them are equal for the two beams, the coupling terms of the IRs simply add together in Eq. (13). If one of these identical IRs is tilted by 90° relatively another, both the total incoherent and the total coherent beam-beam terms are cancelled, provided the phase advance between them of the beam one is equal to that of the beam two. If the phase advances are not equal, the incoherent terms vanish, but coherent do not. The last situation was realized at the LHC, with its orthogonal beam crossing planes for the two main IRs and significant difference between the inter-IR phase advances of the beam one and beam two, 1 2     [10,11]. Thus, for the LHC case, the dynamic equations (13) have to be modified to These equations yield an eigensystem for transverse dipole multi-bunch oscillations of one or two coupled beams for the given wakes/impedances, gain frequency profile and the collision scheme. Eigenvalues q give the total coherent tune shifts resulted from a combined action of the single-and multi-bunch wakes, damper and beambeam response.
For the LHC, the Nested Head-Tail program is usually run with 5 radial rings, 21 azimuthal head-tail basis functions, Eq.(1), and 15 representative inter-bunch modes, yielding 5 21 15 1575    collective modes for a single beam and twice more for two coupled ones. After computation of the single-bunch impedance matrix Ẑ , which needs to be done only once for the given impedance model, the solution of the eigensystem problem takes about 1 second with Wolfram Mathematica 7 or higher installed on an average laptop.

LANDAU DAMPING
The coherent tune shifts q are found above under a condition of pure harmonic oscillations for all the three degrees of freedom. In general, anharmonicities significantly change mode structure, driving it far from the harmonic case. If so, the considered harmonic solution appears to be useless, and the problem has to be solved from scratch. However, if the anharmonicities are not that large, they can be treated as perturbations to the eigenvalues q; the latter can be used as zero approximation. At first approximation, a tune shift driven by the perturbation can be conventionally expressed as a diagonal matrix element of the perturbation term with the unperturbed basis. If the absolute value of this tune perturbation is small enough, the perturbation approach is justified.
For arbitrary non-degenerate matrices, this problem was solved in Ref. [12]. This solution can be expressed as follows. Let Â be a non-degenerate matrix, and , VU be eigenvectors of this matrix and its Hermit conjugation: †*; With all different eigenvalues  , the two sets of vectors  This algebraic result suggests treating of incoherent nonlinearities as perturbations, giving rise to the question about its justification. Nonlinear terms yield an incoherent frequency spread causing Landau damping. When non-linearity is high enough, the beam is stable. The instability threshold corresponds to a case when the growth rate computed for a purely linear situation is equal to the Landau damping rate. Thus, at the threshold, the nonlinear perturbation of the coherent tune is equal to its imaginary part. That is why justification of the perturbation method for the threshold computation requires the imaginary part of the coherent tune shift q to be small compared to its real part. Assuming this is true, the method of perturbation is justified for threshold computation.
where  stays for a perturbed mode frequency of the nonlinear system, matrix M combines all the collective response matrices generating a coherent field ˆ MX , and the coherent vector X is a beam-average of the beamlet perturbations x .
From here, the beamlet perturbation x can be expressed in terms of the collective perturbation X. After the beamaveraging, ... , this yields where k X and k Y are the unperturbed eigenvectors, †*( the normalization condition † 1 kk  YX was assumed. Since all the matrices of Eq. (20) are diagonal, the dispersion equation (20) can be further simplified. For every unperturbed eigenvalue q, the corresponding mode tune with Landau damping taking into account can be found from the following dispersion equation Here, ( , , ) not by a naively expected ...Fd  . For an ensemble of non-harmonic oscillators, that rule was derived by H.G. Hereward [13], and explained by him as resulting from an incoherent frequency shift when the coherent field acts on a given particle. In general, Eq. (22) tells that the mode tune shift  is determined not only by its harmonic approximation q, the incoherent spectrum  , and the phase space density F, but also it depends on the eigenvectors X and Y. ;.
For this, and apparently only for this case, the dispersion relation, and thus a stability condition, does not depend on the eigenvectors. Then, the stability condition can be expressed in terms of the stability diagram independent of the eigensystem [14,15]. Namely, Eq. determined by the averaged nonlinearities  for the phase space density F. The same stability diagram is valid for all the modes: as soon as the eigenvalue q  is located above the diagram, the mode is unstable; otherwise, it is suppressed by Landau damping. For Gaussian transverse distribution, and with negligible spread of the synchrotron frequencies, the 2D dispersion integral was found by R. Gluckstern [6]: (26)       Ref. [16].
For the LHC impedance model, the highest coherent tune shifts are comparable to the synchrotron tune if the damper is turned off [17], making the weak head-tail approach to Landau damping (24), (25) marginally applicable in this area of the parameters. However, when the head-tail phase and the gain are sufficiently large, the most critical coherent tune shifts q are reduced several times [5], allowing to rely on the weak-head tail stability diagram (25) without any visible loss of accuracy.

BENCHMARKING
Above, the main ideas and formulas of the NHT code are described. Some NHT results for LHC were compared with the Air-Bag Averaged (ABA) Vlasov solver [18] and BeamBeam3D tracking code [4,19], showing a good agreement for all the examined cases.
In Ref. [19], instability growth rate computed from the BeamBeam3D tracking simulations for LHC parameters is presented as a function of the chromaticity and damper gain. For those simulations, a single 3D-Gaussian bunch per beam and single IP were assumed. The intensity and collision parameters were taken close to the end of the beta-squeeze case. Namely, 10 rms beam radius of the beam-beam separation was assumed, and the computed beam-beam long-range kick was additionally enhanced by a factor of 10, thus simulating 10 identical long-range collisions instead of one. The IP optics was taken as perfectly round, all the octupoles were zeroed, the potential well was supposed to be ideally parabolic, and the doubled nominal impedance of the LHC was implemented.
To make a comparison, NHT computations were fulfilled for the same conditions. Without octupoles and longitudinal nonlinearity, the only source of the Landau damping is a long-range beam-beam octupole term yielding the incoherent tune shift per collision, which can be presented as  In Ref. [19], the instability threshold is shown as a threshold chromaticity Q for certain values of the damper gain. In Fig. 2, these results of BeamBeam3D are presented together with corresponding NHT ones.
To appreciate the agreement between the two sets of results, one has to take into account that at high chromaticity, 10 Q , the stability condition is barely sensitive to Q (see the following section), which significantly amplifies initially small computational errors when the chromaticity is that high. Note that the NHT data at Fig. 2 reflect eigenvalue computation with beambeam on with the following analysis of the stability diagram. Figure 3 shows BeamBeam3D single bunch, no damper, no octupoles growth rates [19] compared with corresponding NHT eigenvalues Im q for the most unstable mode, demonstrating agreement within few percent or better.

LHC: EIGENSISTEM
In the rest of this paper, various capabilities of the Nested Head-Tail program are demonstrated for the LHC at 4 TeV. All the computations are done for the horizontal degree of freedom. LHC horizontal wake and impedance functions are presented in Fig. 4 and 5, as they were provided to the author by N. Mounet [18,20]. According to this model, vertical wake and impedance are similar to horizontal, being slightly smaller. Hereafter, this impedance model is referred to as a nominal one. Various beam-based measurements at the LHC are showing, though, that actual impedances are ~2-3 times higher than the nominal ones (see Ref. [17] and following sections). So far, a reason for this discrepancy is unknown. Strictly speaking, it has to be supposed that not only a scale of the impedance is higher, but its frequency dependence may be considerably different from the model as well. However, due to a lack of knowledge about the real impedance of the LHC, the NHT computations are normally performed just with the doubled nominal impedance and wake functions.  A significant difference between positive and negative chromaticity is not surprising, but non-monotonic dependence over the gain at the negative chromaticity area was not expected. Note a stable area at moderate gain and slightly negative chromaticity.   . While nodamper growth rates (blue) are about doubled by the beam-beam interaction, they are changed much less for the standard gain 1.4 g  . This property of the coherent spectrum is a consequence of the "damper theorem", see Ref. [8].
Up to this point, all the plots show eigenvalues for pure linear particle motion. The next section takes into account octupole transverse nonlinearity causing Landau damping and responsible for certain instability thresholds. Note that maximal growth rates * q shown in several plots above are not generally proportional to the octupole strength, required for stabilization, since this strength depends not only on the growth rate, but on the tune shift Re( ) q  as well.

LHC: INSTABILITY THRESHOLDS
For Landau damping, two families of octupoles are installed in the LHC, normally fed with the same currents o I . According to Ref. [20,21]  Note that Figs. (14)(15)(16) do not take into account beambeam nonlinearity. At first glance, Eq. (27) could be used for the octupole component of the long-range tune shift to be taken into account together with the LHC Landau octupoles. To do this, one has to double this tune shift, since there are two interaction regions, and divide it back by a factor of two for the first or the last bunch of the batch (pacman bunches), so that for the pacman LHC bunches Eq. (27) gives the result. However, this result would not be normally correct, since in the reality the beams are initially separated in the orthogonal plane as well, which is not taken into account by Eq. (27). Analysis of the stability diagram of Gaussian beams provided by Ref. [22] demonstrates that with this additional beam-beam separation, the beam-beam nonlinearity is approximately equivalent to 100A of the Landau octupoles at the end of the beta-squeeze. For positive polarity of the octupoles, these 100A of beambeam nonlinearity go with the same sign as the octupoles, thus increasing Landau damping. For negative polarity, the two nonlinear contributions are of the opposite sign, leading to a collapse of the stability diagram at certain beam-beam separation. When this was realized [23], the initially negative octupole polarity has been inverted.
Single-beam measurements of the instability thresholds made at the high gain and high chromaticity plateau never exceeded 200A of the Landau octupoles [17]. Compared with Fig. 15, it leads to the conclusion that the effective single-bunch impedance of the LHC should be 2.5-3 times higher than the nominal one of Fig. 4. Figure 16 shows the effective octupole current required for stabilization of two LHC beams seeing each other at the end of the beta-squeeze. The effective octupole current is the sum of the Landau octupole current and a contribution of the beam-beam nonlinearity expressed in terms of the equivalent octupole current. According to Ref. [22], the oncoming beam contributes 100A to the effective octupole current for the pacman bunches and twice more for the central ones at the end of the squeeze. It follows then, that about 100A of the Landau octupoles should be sufficient for the stabilization, assuming the machine is operated at the high gain and high chromaticity plateau of Fig. 16, as it normally was. Contrary to this conclusion, at the end of the squeeze a transverse instability was permanently observed, notwithstanding the octupole current was kept at its maximum of 550A [24]. An initial idea that this instability is driven by the coherent (strong-strong) beam-beam effect or some hidden two-beam impedance was refuted by these NHT computations [25], confirmed later by a dedicated LHC beam experiment [26]. To explain this instability, a hypothesis of three-beam instability, or beam-beam-beam effect was suggested, where the third beam is an electron cloud accumulated in a high-beta area of the main interaction regions [27].

SUMMARY
This paper describes Nested Head-Tail Vlasov solver effectively used for high energy beams in LHC, where radial modes, couple-bunch modes, feedbacks, beambeam effects and nonlinearities responsible for Landau damping are accurately handled. Main advantage of that solver against macroparticle tracking codes is many orders of magnitude shorter CPU time, which allows a fast and efficient analysis of that complicated system in a multidimensional space of parameters.