Thresholds for loss of Landau damping in longitudinal plane

Landau damping mechanism plays a crucial role in providing single-bunch stability in LHC, High-Luminosity LHC, other existing as well as previous and future (like FCC) circular hadron accelerators. In this paper, the thresholds for the loss of Landau damping (LLD) in the longitudinal plane are derived analytically using the Lebedev matrix equation (1968) and the concept of the emerged van Kampen modes (1983). We have found that for the commonly-used particle distribution functions from a binomial family, the LLD threshold vanishes in the presence of the constant inductive impedance Im$Z/k$ above transition energy. Thus, the effect of the cutoff frequency or the resonant frequency of a broad-band impedance on beam dynamics is studied in detail. The findings are confirmed by direct numerical solutions of the Lebedev equation as well as using the Oide-Yokoya method (1990). Moreover, the characteristics, which are important for beam operation, as the amplitude of residual oscillations and the damping time after a kick (or injection errors) are considered both above and below the threshold. Dependence of the threshold on particle distribution in the longitudinal phase space is also analyzed, including some special cases with a non-zero threshold for Im$Z/k = const$. All main results are confirmed by macro-particle simulations and consistent with available beam measurements in the LHC.


I. INTRODUCTION
Landau damping [1] is one of the most efficient mechanisms of beam stabilization in both transverse and longitudinal planes of particle motion in hadron synchrotrons. Landau damping is considered to be lost when the frequency of the coherent bunch oscillations moves outside the incoherent frequency band.
In the longitudinal plane, Landau damping of the coherent modes is provided by synchrotron frequency spread coming from a non-linearity of the rf field, modified by a selfinduced wake-field. Synchrotron frequency spread can be increased by increasing the rfbucket filling factor or longitudinal emittance and by using a higher-harmonic rf system. All these ("passive") methods are applied in CERN synchrotrons together with some "active" damping systems of beam instabilities (feedback systems).
Undamped coherent beam motion observed at Tevatron [2], RHIC [3], SPS [4] and LHC [5] can be attributed to loss of Landau damping (LLD). In LHC, in the absence of a longitudinal wide-band feedback system and a higher harmonic rf system, single bunch stability is guaranteed by controlled longitudinal emittance blowup during the acceleration cycle. However, it was noticed from the start of the LHC operation that the beams with the nominal parameters are at the limit of stability at different moments of the acceleration cycle due to LLD, with thresholds about four times lower than expected [6]. An agreement of first macro-particle simulations, based on the existing impedance model, with the measurements [7] suggested that the discrepancy cannot be explained by a large missing impedance. The initial choice of the longitudinal beam parameters in the LHC [8] was based on the criterion derived from the Sacherer stability diagrams [9] (see Ref. [10] for more details). While in the longitudinal plane these diagrams are justified only for a coasting beam [11] or a very low-frequency impedance [6], the Sacherer criterion is still widely used for the analytical threshold evaluation, also often assuming a constant impedance ImZ/k, with k being a harmonic of the revolution frequency.
A general way to study beam stability is to use the Vlasov equation linearized for a small perturbation of a stationary particle distribution function. The first self-consistent system of equations suitable for the eigenvalue analysis of longitudinal beam stability was proposed by A. N. Lebedev in 1968 [12]. Converting it to an integral equation, he has also proven that a space-charge impedance cannot cause instability of a bunched beam. Nevertheless, as was shown later, a reactive impedance may lead to LLD.
In the early 70's, F. Sacherer derived another matrix equation, from which then the stability diagrams were obtained using the following assumptions: (i) the rf field is linear; (ii) the actual self-force (from beam interaction with its environment) depends only on the location of the bunch center-of-mass (so-called "synthetic kernel" approach describing the rigid-bunch motion); (iii) the synchrotron frequency spread is neglected in the calculation of coherent frequency shift with intensity. The LLD threshold is reached when the coherent mode crosses the boundary of the stability diagram. The LLD threshold, similar to the Sacherer analytical criterion, was also obtained from the integral form of the Lebedev equation [12] in Ref. [13] for a constant ImZ/k. However, as was shown recently in Ref. [6], the Sacherer stability diagrams can be obtained from the Lebedev equation using the low-frequency approximation.
Another approach was suggested by A. Hofmann and F. Pedersen [14] and used more recently, for example, for analysis of Landau damping in single [15] and multi-harmonic rf systems [16]. It is based on direct comparison of the coherent dipole frequency of the rigidbunch oscillations with the maximum incoherent frequency inside the bunch, both calculated for a constant reactive impedance ImZ/k and an elliptic particle distribution in longitudinal phase space (corresponding to the so-called "parabolic" bunches). Contrary to the Sacherer stability diagrams, here the synchrotron frequency spread is taken into account in the calculation of coherent frequencies. For the same conditions, the LLD threshold obtained by this method is only about 20-30% higher than the one from the stability diagram [9]. However, comparison with measurements in the LHC indicates that both these approaches may significantly underestimate the real threshold [6].
There are several ways to find the accurate solutions of the linearized Vlasov equation without neglecting the synchrotron frequency spread [12], [17], [18], [19]. In 1983, Y. H. Chin, K. Satoh, and K. Yokoya introduced the concept of van Kampen modes [20,21] for a description of bunch oscillations in the longitudinal plane [18]. Landau damping at low intensities was explained as decoherence of the van Kampen modes, while the LLD threshold is reached, once a van Kampen mode emerges from the band of the incoherent synchrotron frequencies. In 2010, A. Burov [22] used this concept together with an eigenvalue approach, suggested for analysis of single-bunch instabilities by K. Oide and K. Yokoya [19], to evaluate the LLD thresholds for different impedance models and rf potentials, also taking into account a potential well distortion. In particular, he found that the LLD threshold calculated for the constant ImZ/k is more than a factor of three lower than the one obtained from the Sacherer stability diagrams [23].
Even though in all previous studies a finite LLD threshold was obtained for a constant ImZ/k, we will show below that in most cases this threshold is zero unless some cutoff frequency is introduced (which effectively is also the case in all numerical calculations or particle simulations). The thresholds and the analytical criterion for the LLD are obtained by applying, for the first time, the concept of van Kampen modes to solve the Lebedev equation in its original matrix form [12]. These results are compared with those obtained using the Oide-Yokoya method, macro-particle simulations, and the beam measurements in the LHC.
The paper is organized in the following way. After this Introduction, the main equations of longitudinal motion and solutions for stationary particle distribution with intensity effects included, are presented in Section II. For completeness, in Section III, both the A similar analysis is also performed for a broad-band resonator model and the results are compared with macro-particle simulations using CERN Beam Longitudinal Dynamics code BLonD [24]. The influence of particle distribution on the threshold value is also investigated here. In Section V, the meaning of the LLD threshold for an accelerator operation is analyzed in detail by considering a bunch response to the most common, rigid-dipole, perturbation (a kick due to a phase/energy error or noise excitation). Using the van Kampen modes as a basis for an expansion of the rigid-dipole mode, the impact of this perturbation on the beam is evaluated semi-analytically and then compared with macro-particle simulations. In particular, it is demonstrated that the damping time of a perturbation below the LLD threshold and the amplitude of the residual bunch oscillations above this threshold are the most important characteristics of beam behavior. Moreover, the obtained results are also compared with available beam measurements in the LHC. The LLD thresholds for some interesting cases from a parameter space, which was not covered in the previous sections, are briefly discussed in Section VI. This includes the reactive impedance with the opposite sign (capacitive above transition or inductive below) and a special class of particle distribution functions, both leading to a finite LLD threshold even for a constant ImZ/k impedance. The main conclusions of this work are presented in Section VII. In the end, Appendixes A and B contain some important derivations.

II. LONGITUDINAL MOTION WITH INTENSITY EFFECTS
We start with a stationary situation, while time-dependent perturbations are treated in the following section. We chose the coordinate system with respect to the synchronous particle with the energy E 0 and the phase φ s0 that corresponds to the case of a single rf system without taking into account intensity effects. Then, the main equations of longitudinal particle motion can be written using conjugate variables {∆E/(hω 0 ), φ}, d dt where ∆E and φ are respectively the energy and phase deviations of the particle, ω 0 = 2πf 0 , f 0 is the revolution frequency, β is the particle velocity normalized by the speed of light, h is the harmonic number, q is the electrical charge of the particle, δE 0 is the energy gain per turn of the synchronous particle, η = 1/γ 2 tr − 1/γ 2 is the slip factor, γ is the relativistic Lorentz factor, and γ tr is the Lorentz factor at transition energy. These equations can be found in standard textbooks (e.g. [25]).
The present work focuses on the case of a single bunch in a single rf system, while the derivations can be easily adapted to any combination of rf waves. In the presence of intensity effects, the total voltage V t (φ) = V rf (φ) + V ind (φ), in addition to the rf voltage V rf , contains contributions from the beam-induced fields described as the beam-induced voltage V ind .
Then, V rf (φ) = V 0 sin(φ s0 + φ), with V 0 being the rf voltage amplitude, and the synchronous phase φ s0 = arcsin (δE 0 /qV 0 ) for η < 0 or φ s0 = π − arcsin (δE 0 /qV 0 ) for η > 0. The stationary beam-induced voltage can be written in the form (see Appendix A), where N p is the total number of particles in the bunch, Z k = Z(kω 0 ) is the longitudinal impedance at frequency kω 0 , and λ k is the Fourier harmonic of the normalized line density, The line density is defined by the distribution function , and the following normalization is imposed It is convenient for further analysis to use another set of variables (E, ψ), which correspond respectively to the energy and phase of the synchrotron oscillations, Here ω s0 = 2πf s0 is the angular frequency of small-amplitude synchrotron oscillations in a bare rf potential with The total potential can be obtained from the total voltage V t with the synchronous phase shift due to intensity effects ∆φ s that satisfies the relation The dependence of the synchrotron frequency on the energy of synchrotron oscillations ω s (E) = 2π/T s (E) in Eq. (7) can be found from the period of oscillations where φ min (E) and φ max (E) are the minimum and maximum phases of the particle with the energy of synchrotron oscillation E which satisfy the relation E = U t (φ).
The stationary particle distribution function F is a function of only the energy of synchrotron oscillations E with the line density It depends on the total potential that in calculations can be found using an iterative procedure, similar to the one used in Ref. [22] (see Appendix B). Below we will mostly consider distributions belonging to a binomial family, which covers a wide range of realistic bunch shapes, from flat (µ = −1/2) to Gaussian (µ → ∞), where the normalization factor In this case, the integration in Eq. (11) can be performed analytically yielding where Γ is the gamma function.
For µ → ∞ a bunch has a Gaussian line density and the corresponding bunch length τ 4σ is typically defined as four times the Root-Mean-Square bunch length σ, i.e. τ 4σ = 4σ.
The bunch length τ 4σ can be easily calculated from the Full-Width Half-Maximum (FWHM) bunch length τ FWHM : In practice, the proton bunches are usually not Gaussian and have a finite full length τ full defined as However, we will also use definition (15) below as it has important features related to the LLD threshold (see Sec. IV).

III. PERTURBATION FORMALISM
For analysis of beam stability, we need to consider the perturbationsF,λ andṼ ind , respectively, to the equilibrium distribution function F, line density λ, and induced voltage V ind , introduced in the previous section. If these deviations grow with time, the beam is unstable.
The linearized Vlasov equation in (E, ψ) variables is where by definition dψ/dt = ω s (E), while dE/dt can be found from a combination of the first and second equations of particle motion in the presence of perturbation Multiplying both sides of Eq. (18) byφ and taking into account Eq. (6) one obtains where the perturbed induced potentialŨ ind (φ, t) is defined similarly to Eq. (9) Finally, the linearized Vlasov equation can be written in the form Several approaches allow finding general solutions of this equation without neglecting the synchrotron frequency spread. In the following subsection we derive, the equation that was proposed by A. N. Lebedev in 1968 [12], but since then it was never used in its original matrix form for evaluation of LLD for the single-bunch case.

A. Lebedev equation
Here we derive the matrix equation first obtained in Ref. [12]. The detailed derivation with slightly different variables can be also found in [26]. The solution of Eq. (21) must be periodic in ψ, and can be represented as the sum of the harmonics e imψ (m = 0).
Presenting the perturbations for a given Ω asF(E, ψ, t) =F(E, ψ, Ω)e −iΩt andŨ ind (E, ψ, t) = U ind (E, ψ, Ω)e −iΩt , we havẽ with Then the solution of Eq. (21) is The perturbed induced voltage is related to the perturbed line density as (see Appendix A) where Z k (Ω) = Z(kω 0 + Ω). Then, using Eq. (20), the perturbed induced potential can be presented in the form and its Fourier harmonicsŨ ind,m , defined by Eq. (23), can be obtained as Here we introduce the dimensionless "intensity" parameter ζ, Some constant value ImZ/k, used for the normalization of the impedance, is well defined for the case of constant inductive impedance considered in the present work, while it can be arbitrarily chosen for other impedance models.
The function I mk (E), introduced for the first time in [12], plays an important role The last expression comes from the fact that φ(E, −ψ) = φ(E, ψ). Note, that I mk depends on ζ because φ(E, ψ) is modified via potential well distortion. It is easy to see also the following properties In the case of symmetric potential well, we also have I * mk = (−1) m I mk . An alternative way to evaluate I mk is to perform in Eq. (29) integration by parts where the phase of synchrotron oscillations is given by Eq. (7).
For perturbation of the line density at frequency Ω the harmonicsλ k (Ω) are related toF(E, ψ, Ω) as Here, the transformation of variables dφdφ = ω 2 s0 dψ dE/ω s (E) was used. At the next step, inserting Eq. (24) into Eq. (32), an infinite system of equations for harmonics of the line density perturbation can be obtained where the beam transfer matrices [26] are defined as After substitution of F from Eq. (12) the element G pk has a form Note that Eq. (33), referred to below as the Lebedev equation, is the general equation since no approximations were used so far. The elements G pk depend on intensity parameter ζ as they are found after the stationary problem is solved. The solution for Ω exists if the determinant of the following matrix is zero

B. Oide-Yokoya method
In this section, we re-derive, for completeness, another matrix equation that allows finding the bunch modes. According to this method [19], called later the Oide-Yokoya, the perturbed distribution functionF is expanded as Inserting this expansion into Eq. (21) [27], multiplying both sides by cos mψ or sin mψ and integrating over ψ, the following system of equations can be obtained where the symmetry ofŨ ind visible in Eq. (26) was also used. Inserting expansion (37) into Eq. (32), the perturbed line density harmonics can be presented in the form Combining Eqs. (38-40), one finally gets In the original work [19], the functions C m (E, Ω) were defined as a combination of the step- where E n is the n-th mesh point on the energy grid (in our calculations, it is assumed to be uniform with the total number of points N E , i.e. ∆E n = E max /N E ). Thus, we obtain an eigenvalue problem of a linear algebra: with the matrix elements δ nn is the Kronecker delta, m max is the maximum value of the azimuthal mode number, and k max is the maximum value of the revolution harmonic number. It was shown in Ref. [28] that approximating the integration over energy in Eq. (41) by a sum leads to the same eigenvalue problem. In general, all azimuthal modes and all frequency harmonics need to be included in calculations, i.e. m max → ∞ and k max → ∞. In practice, a truncation is used that might affect the results, which will be discussed in the following section. Note, that in the matrix element M nmn m dependence of impedance Z k on Ω was neglected, so that an eigenvalue problem can be solved.

IV. LOSS OF LANDAU DAMPING
In this section, we will consider the LLD thresholds in a general case and in application to the LHC, whose main parameters are listed in Table I. For the LHC examples, with the beam being above transition energy, we will focus on a non-accelerating case with φ s0 = π, and use the rf voltage that corresponds to the operational value at the flat bottom.  [20,21] has found that they have a continuous and a discrete part. In our case, at low intensities (ζ ≈ 0), all van Kampen modes belong to the continuous spectrum, Ω = mω s (E), and the corresponding eigenfunctions are singular.
Landau damping results then from the phase mixing of these modes, which do not represent the collective motion of the particles. Above the threshold the discrete van Kampen modes can emerge from the continuous spectrum, implying that Landau damping is lost [18]. These modes are described by regular functions and by definition they lie outside ω s (E). The Oide-Yokoya method [19] used to calculate numerically the van Kampen modes allowed to determine the LLD threshold in different situations [22].
A. The LLD threshold for constant inductive impedance

Analytic criteria
Here we first derive an analytic expression for the LLD threshold using the Lebedev equation (33). Then it will be compared with the results of the semi-analytic calculations using both matrix equations described in Section III as well as with the macro-particle simulations. In what follows we will focus on the case of ηImZ/n > 0 (space charge below transition or inductance above transition), while the opposite sign is briefly discussed in Sec. VI. I mk (0) = 0, the integral (35) defining the elements G pk converges for all p and k.
To find the solution of Eq. (36), we can use the following property of the matrix where tr(X) is the trace of an arbitrary square matrix X, and ε is the small parameter ε 1. Expansion of the exponent up to the first order of ε yields, with the identity matrix I. Thus, we get a general expression for the LLD threshold Naturally, the parameter ε ∝ ζ, while its dependence on the bunch length will be deduced below. Thus, potential well distortion has to be neglected in calculation of G kk as it would be a higher-order term of the parameter ε.
The element G kk can be calculated analytically if we consider short bunches (E max 1) in a single rf system.
If only the first azimuthal mode (m = 1) is taking into account, from Eq. (35) we can obtain Here we used a new variable x = E/E max and defined the maximum phase amplitude Below, we will consider the distributions of the binomial family g(x) = (1 − x 2 ) µ for which the normalization factor A N = φ 2 max /(2µ+2) can be found using Eq. (13), once the synchrotron frequency spread is neglected. Then G kk is where p F q (a 1 , .., a p ; b 1 , ..., b q ; z) is the generalized Hypergeometric function and y = kφ max /h. Considering the constant inductive impedance Z k = ik ImZ/k, the sum in Eq. (46) can also be analytically evaluated by approximating it with an integral which diverges, except the case of µ = 0 (a water-bag distribution). The derivative of this distribution function F ∝ δ(E − E max ). Thus the integrals in Eq. (34) can be evaluated analytically (without any approximations), yielding that each element G pk (Ω) is proportional to ω s (E max ). Since, ω s (E max ) = 0 due to the potential well distortion (see Appendix B), the determinant D(Ω, ζ) is never zero and the LLD threshold does not exist. For other values of µ, once we have truncated the sum at some k max , the LLD threshold is or in terms of intensity where y max = k max φ max /h and we introduced the function Examples of this function for two different µ values are shown in Fig. 1.
For µ = 1/2, the Sacherer formalism proposes the LLD threshold (in our notations) [10] Similarly, the Hofmann-Pedersen approach [14] applied for short bunches yields, The same expressions can be obtained from Eq. (50) for y max ≈ 3.32 and y max ≈ 2.95, respectively, (see also Fig. 1). This implies that to have these thresholds, the impedance each time has to be truncated at the frequency f c = f 0 k max ≈ 1/τ full , inversely proportional to the full bunch length τ full = 2φ max /ω rf . This is a low-frequency approximation while the constant inductive impedance Z k /k does not decay for k max − → ∞. correspond effectively to that used in Sacherer [9] and Hofmann-Pedersen [14] criteria are shown with squares.
Since the generalized Hypergeometric function 2 F 3 approaches zero for y → ∞ a simple expression for the LLD threshold can be obtained Thus, we see that the threshold is inversely proportional to the cutoff frequency and the fifth power in the dependence on the bunch length is replaced by the fourth. Our analytic formula (ε 1) is justified for impedance with k max h/8 (for µ ≥ 1/2).

Calculations using code MELODY
The analytic threshold (50) can now be compared with semi-analytical results obtained using code MELODY (Matrix Equations for LOngitudinal beam DYnamics calculations) [29]. In calculations based on the Lebedev equation (36), the determinant D(Ω, ζ) is numerically evaluated and the threshold ζ th corresponds to D(Ω, ζ th ) = 0. This is illustrated in Fig. 2 (left) for the LHC parameters and µ = 2 observed in measurements [7].
In Oide-Yokoya method, the eigenvalues of the matrix in Eq. (43) are calculated as a function of the intensity parameter ζ. To find the threshold, the difference between the maximum eigenfrequency and the maximum incoherent frequency is evaluated. The threshold corresponds to the zero of this function, as in an example in Fig. 2 (right).
In Fig. 3 we compare the LLD threshold as a function of the full bunch length calculated for two different cutoff frequencies using analytic Eq. (50) and code MELODY. Indeed, the analytic expression derived from the Lebedev equation agrees well with the corresponding exact solution evaluated semi-analytically. As expected there is some discrepancy for larger bunch lengths since the analytic threshold was derived in short-bunch approximation while still taking into account the synchrotron frequency spread. Moreover, one can see that numerical results obtained from MELODY using the Oide-Yokoya method and the Lebedev equation agree almost perfectly, with the maximum 2% relative error in the covered bunchlength range. We also observe that dependence on the bunch length is even slightly smaller than in the fourth power.
The semi-analytic results confirm that for a given bunch length, the threshold reduces as a function of the cutoff frequency (Fig. 4). The overall agreement between the analytical criterion (50) and the semi-analytic results is very good.
We also evaluated the LLD thresholds as a function of the scaled bunch length τ 4σ for a fixed cutoff frequency and different powers µ of the binomial distribution. The analytic prediction (50) slightly underestimates the LLD threshold evaluated semi-analytically, see  for an arbitrary impedance Z k /k if we define the following effective impedance where k eff is the effective cutoff frequency, above which ImZ k /k is always negative, i.e.
ImZ k /k < 0 for k > k eff . Then, the LLD threshold becomes

and again in terms of intensity
Here, in comparison to Eq. (51), ImZ/k is substituted by (ImZ/k) eff . Its applicability to an arbitrary impedance model is a subject of further studies, while in the present work we will focus on the broad-band resonator impedance where R is the shunt impedance, Q is the quality factor, and f r is the resonant frequency. To simplify comparison with the model of the constant inductive impedance, we choose Q = 1, R = (ImZ/k)Qf r /f 0 , with the same ImZ/k = 0.07 Ω as in the previous section while an impact of resonant frequency f r on the threshold is studied. Figure 6 shows dependence of the LLD threshold on the full bunch length for two values of the resonant frequencies f r , similarly to Fig. 3. The analytical threshold given by Eq. (55) agrees well with the results of semi-analytic calculations. Dependence of the LLD threshold on the resonant frequency is also well reproduced for a given bunch length (Fig. 7), as for the case of the truncated constant inductive impedance, see Fig. 4. Note that these LLD thresholds were found taking into account only one azimuthal mode m. This is sufficient for large f r τ full values as the LLD intensity thresholds are lower for them and the contribution of higher-order azimuthal modes can be neglected. For smaller f r τ full , the contribution of higher-order azimuthal modes reduces the LLD threshold by a few percents.
The results of semi-analytic calculations were compared with the macro-particle simulations using the CERN Beam Longitudinal Dynamics code (BLonD) [24] . We looked at the spectrum of the turn-by-turn data of different moments of the particle distribution depending on the azimuthal mode of interest (e.g., the mean value is sufficient to evaluate the frequencies of the dipole mode).
The bunch with 2×10 6 macro-particles was generated, matched with intensity effects and then tracked for a sufficient number of turns to have a proper frequency resolution of the spectrum (typically up to 10 6 turns). To evaluate induced voltage, we used at least 256 slices per rf bucket, that was sufficient to avoid uncontrolled emittance blowup due to numerical noise. Since at the LLD threshold, the mode frequency equals the maximum incoherent frequency it cannot be properly observed with a finite number of turns. However, above the   the van Kampen mode is seen as a strong peak above the incoherent band in Fig. 8.
Applying the Oide-Yokoya method in semi-analytic calculations by MELODY, the coherent modes above the threshold can be found in a standard way as solutions of the eigenvalue problem. On the other hand, to find them using the Lebedev equation, we fix the intensity parameter ζ and vary Ω >ω s until D(Ω, ζ) = 0. A good agreement between macro-particle simulations and semi-analytic calculations can be observed: as predicted the LLD threshold is lower for a higher resonant frequency. This is also the case for a higher-order radial mode of the dipole mode seen as the second emerged mode in Fig 8 (right).
As it will be shown below, even though the LLD threshold is lower for impedance with higher cutoff or resonant frequencies, it does not necessarily mean that an initial perturbation will have a stronger impact on the beam at the threshold.

V. IMPACT ON THE BEAM
A rigid-bunch dipole perturbation (a kick) could be a result of the injection phase or/and energy errors, but also of a phase noise in the rf system. Evaluation of the beam response to a rigid-dipole perturbation is a common way to study the LLD in simulations and measurements. N. G. van Kampen has shown that an arbitrary perturbation of the initial distribution behaves as a superposition of the waves that are solutions of the equation equivalent to Eq. (41) [20]. Below we will use expansion of a rigid-dipole perturbation on the basis of van Kampen modes to obtain analytically its time evolution. It was also used to describe the transverse oscillations of the colliding bunches [30]. Note that a similar approach can be applied for the analysis of higher-order perturbations (e.g., quadrupolar oscillations due to bucket mismatch).

A. Rigid-dipole kick as superposition of van Kampen modes
As a result of a kick, the whole bunch has a phase offset ∆φ with the respect to its synchronous phase. Since λ(φ + ∆φ) = λ(φ) + ∆φ dλ/dφ, the rigid-bunch dipole mode can be described via the derivative of the line density dλ/dφ. Its spectral harmonics (dλ/dφ) k can be obtained using Eq. (4) Dependence of harmonics on kI 0k means that it is not a pure dipole (m = 1) mode. The elements of the eigenvector of this mode can be approximated as C rd 1 (E) ∝ √ EF (E) and C rd m (E) = 0 for m > 1 [27]. We can define the harmonics of the approximate rigid-dipole mode,λ rd k , similar to Eq. (40), After finding eigenvalues and eigenvectors of matrix (43), this approximate mode can be expanded using the van Kampen modes as a basis where α n m is the expansion coefficient, Ω m n is the van Kampen mode with the index n that belongs to azimuthal mode m . We impose the following normalization of the eigenvectors Since the basis is not orthonormal, to find α nm , we construct a matrix from eigenvectors K mnm n = C m (E n , Ω m n ), evaluate its inverse matrix K −1 , and finally obtain B. Bunch offset evolution after a rigid-dipole kick Once the solution of the eigenvalue problem (43) and the mode expansion coefficients a mn from Eq. (61) are found, the bunch offset evolution can be described as where the factor κ = max [λ(φ + ∆φ) − λ(φ)] / max λ rd (φ) depends on the initial offset ∆φ. In what follows we will consider the broad-band resonator impedance model, while similar results can also be obtained for the truncated constant inductive impedance.
Considering a small initial kick ∆φ ω rf τ 4σ that does not cause a significant emittance blowup, we can evaluate analytically, using Eq. (62), the bunch evolution with time. The results for different intensities are shown in Fig. 9. The bunch oscillations are damped after the initial fast decoherence when intensity is lower than the LLD threshold. There are residual oscillations above the LLD intensity threshold and their amplitudes depend on intensity. At even higher intensities, one typically observes a beating.
We can also see that the intensity dependency of the damping time is similar for different resonator frequencies (left plot in Fig. 10). The damping time is short for N p N p,th , while it approaches infinity when the LLD threshold is reached. It means that Landau damping can be effectively lost even below the threshold depending on the machine cycle length.
Moreover, at the LLD threshold and above, the residual oscillation amplitude after a kick depends on the resonant frequency. Comparing the average residual oscillation amplitudes, we see that, at a given intensity, it is smaller for a higher f r . We also observe that amplitude dependence on intensity saturates for f r τ 4σ 1. These findings become more clear by analyzing the mode patterns in the frequency and time domains.
Comparison of coherent mode spectra for different values of bunch length and resonant frequency in Fig. 11 (left) shows that at the LLD threshold, the mode is not sensitive to details of particle distribution, but rather depends on impedance and its resonant frequency. This, in fact, explains why the LLD threshold significantly increases for f r τ full → 1 as can be seen in Fig. 7. For sufficiently low f r , the characteristic width of the mode becomes larger than the bunch length, so the mode cannot exist inside the bunch. We also observe that the mode is "wider" in the frequency domain for a higher f r , which would correspond to a more localized mode in the time domain. This explains the fact that beam response to the kick is stronger for lower f r . On contrary, if we compare the same modes well above the LLD threshold (at N p = 4 × 10 11 ppb, for example), they become very similar for the same bunch lengths but different f r (see right plot in Fig. 11). Typically, the coherent modes above the LLD threshold widen in the time domain as intensity increases (see examples in Fig. 12), so the beam has a stronger response to a rigid-dipole offset as shown in Fig. 10.

C. Comparison with simulations and LHC measurements
Below we will compare the results obtained so far from semi-analytic calculations using code MELODY with macro-particle simulations by code BLonD. In simulations, the amplitude of oscillations from the turn-by-turn bunch offset was extracted using the Hilbert transform. An example of comparison is shown in Fig. 13 for f r = 4 GHz. One can see a good agreement above the threshold, even though we used the approximate rigid-dipole mode in the analytic mode expansion. The exact phase of beating slightly differs, but the amplitudes are very similar. Below the threshold, where we expect a slow damping around the threshold, in the simulations we observe the effect of noise even when 20 × 10 6 macro-   [8]. The rf voltage is 10 MV, as during the measurements in the LHC. We also denote the bunch parameters for those the LLD was observed [7]. One can obtain from the fit that the residual oscillation amplitude scales as τ 6 4σ , while the LLD threshold is proportional to τ 4 4σ . The parameters from the measurements are in the region where Landau damping is lost, but residual amplitude is below 20% of the kick amplitude, in agreement with observations. In the future, however, the LHC impedance model at high frequencies has to be revised to have more accurate predictions of the LLD threshold for the HL-LHC beam intensity.

VI. OTHER CASES: DISCUSSION
The present work was dealing so far with the binomial particle distributions and the inductive impedance above transition energy (ηImZ/k > 0), which is typical for high energy colliders. We have seen that the LLD threshold is decreasing with cutoff or resonant frequency of impedance. Here we will consider two cases, when the LLD threshold does not depend on these frequencies.  [7] are shown by a red box.

A. Flat distributions
Previous studies have pointed out that the LLD threshold is sensitive to the smallargument steepness of the distribution function [23]. Based on this fact, the dedicated measurements were performed at the Tevatron that allowed suppressing the bunch oscillations due to LLD [31]. Our studies show that the LLD threshold is indeed higher for the distribution function with g (E = 0) = 0, used, for example, in Ref. [32] g where 0 < a < 1. For this case, the diagonal matrix elements G kk (48) are For a very large argument y → ∞, they simply become and approach zero as 1/y. This is the main difference in comparison to the particle distributions of the binomial family, where the elements G kk saturate at some constant level.
Calculating the sum of G kk to evaluate the LLD threshold using Eq. (46), one gets a weak (logarithmic) dependence on the cutoff frequency. However, the semi-analytic calculations with MELODY show that the LLD threshold does not depend on the cutoff frequency, as can be seen in Fig. 17. The discrepancy can be understood after numerical evaluation of a small parameter ε = ζ/φ 4 max , which, in this case, is significantly larger than for the binomial distribution and thus expansion (45) cannot be fully justified. Nevertheless, we propose the following analytic LLD threshold, based on the asymptotic behavior of G kk in Eq. (65), It agrees well with the results from semi-analytic calculations (see Fig. 17). In particular, the LLD threshold for distribution (63) with a = 1 is by orders of magnitude higher than the one for the binomial distribution with µ = 2 shown in Fig. 3. We have also found that once the LLD threshold is reached, the response to a rigid-dipole perturbation can be stronger than for a bunch with the binomial distribution. Note, in operation this distribution was obtained as the result of the specific rf manipulations (e.g., rf phase modulation) [31,33].

B. Constant inductive impedance below transition energy
Here, we will show that the LLD threshold does not depend on the cutoff frequency for the case of the inductive impedance below transition or the space charge above transition (i.e. η ImZ/k < 0). In this case, at the threshold, the van Kampen mode emerges below the minimum incoherent frequency, and we can calculate elements G kk for Ω = ω s0 (1 − E max /8), for µ > 1. For a very large argument y → ∞, they can be expressed as and also approach zero, similar to the case of flat bunches. Thus we expect that the LLD threshold does not depend on the cutoff frequency. Indeed, this can be seen in Fig. 18, where the semi-analytic calculations using MELODY are shown for different values of µ. There is a strong dependence on the tails of distribution and for the moment only fitted thresholds are proposed: |ζ th | ≈ 0.034φ 5 max for µ = 1.5, and |ζ th | ≈ 0.067φ 5 max for µ = 2.

VII. CONCLUSIONS
Loss of Landau damping (LLD) in the longitudinal plane can be an important performance limitation of the existing and future storage rings. In the present paper, the criterion of the emerged van Kampen mode was used to determine the LLD thresholds for different beam and machine parameters. In particular, we were able to derive a general analytic expression for the LLD threshold of the dipole oscillations in a single rf system exploiting the Lebedev equation.
Contrary to the previous studies, we have found that for a particle distribution of the binomial family, a constant inductive impedance ImZ/k above transition (the LHC case) or capacitive (space charge) below leads to a zero LLD threshold. In fact, it becomes inversely shown for the binomial distribution that even though the LLD threshold is lower for higher cutoff/resonant frequencies, an impact on the beam at the LLD threshold of impedance with higher frequencies is smaller, also with a strong dependence on the impedance model. On the other hand, for intensities well above the LLD threshold, the beam response is mostly defined by the bunch parameters. Again, these results were confirmed by simulations that required to use tens of millions of macro-particles and sufficient slicing in the induced voltage calculations to cope with the numerical noise. We also have found that the amplitude of bunch oscillations after a phase offset has a strong dependence on the bunch length (in the sixth power). Finally, our calculations are in good agreement with available measurements in the LHC.
There is also a special class of particle distribution functions, with zero derivative in the center of the bunch, which leads to a high and finite LLD threshold for a constant inductive impedance ImZ/k above transition energy (and capacitive below). Moreover, we also shortly discussed the case of a constant inductive impedance ImZ/k below transition energy (space impedance above) that results in the finite LLD threshold for the binomial distributions. A similar analysis can be also applied for higher-order azimuthal modes and different forms of the bare rf potential (double rf system). This could be a subject of further studies.
Below we will derive the induced voltage for the stationary and the time-dependent cases.
The stationary part of the induced voltage for a single bunch with intensity N p is related to the normalized line density λ(φ) as [25] V where the normalization πh −πh dφλ(φ) = 1 was imposed. Because of causality, sum can also be extended for k < 0. Inserting the relation between wake function and impedance (A2) one gets Using here the Dirac comb relation we can present the induced voltage as where Z k = Z(kω 0 ) and the harmonic of the line density Assuming the perturbation of the line density in the formλ(φ, Ω, t) =λ(φ, Ω)e −iΩt , the induced voltage can be expressed using the convolution of line density with wake function, similarly to Eq. (A3), [25] V ind (φ, t) = −qN p In the present work we consider the case when Ω ≈ mω s0 ω 0 for relatively low azimuthal modes m = 1, 2, 3, ..., so the phase factor e ±iΩφ/ω rf in Eqs. (A11, A12) is neglected for the rest of derivations.
The solution of this system of equations, if it exists, can be found for a sufficiently small convergence parameter > 0.
In the case of the constant inductive impedance Z k = ikImZ/k, Eq. (B1) can be simplified for the binomial distribution by using Eq. (14), to obtain an implicit form of the total Here we used the dimensionless parameter The total potential can be obtained analytically for µ = 1/2 as in Hofmann-Pedersen approach [14], but also for µ = 0, 1 and 3/2, for example, where either quadratic or cubic equation needs to be solved. In general, the first derivative can easily be found as where The second derivative that defines the small-amplitude synchrotron frequency is All higher-order derivatives can be calculated recursively.
The case of µ = 0 is of particular interest. From the above equations one can see that all derivatives vanish at the maximum and minimum particle excursions in the bunch, φ min and φ max , respectively. In means that the synchrotron period (10) is infinite for that trajectory and thus ω s (E max ) = 0.