On a model for weak turbulence in synchrotrons

In this work we formulate a theory of weak turbulence for the longitudinal degree of freedom in a stored, coasting beam. We employ a perturbative approach to the fully nonlinear Vlasov equation that is well suited to a set of well separated eigenfrequencies and study the case where the beam is marginally stable to longitudinal oscillations. We derive an averaged equation that describes the nonlinear flow of the fluctuation power among the multiplicity of harmonic modes, characterized by the properties of the machine impedance. In addition, we introduce the Schottky noise of the beam that acts as a thermal reservoir with which the wake-driven oscillations must come into equilibrium. We find steady-state solutions of the nonlinear equation set that indicate significant enhancement of the Schottky spectrum is possible.


Introduction
An outstanding issue in the operation of synchrotrons concerns the nature of equilibrium fluctuations of marginally stable beams. Observations have shown that in stored hadron beams the steady state fluctuation levels are often orders of magnitude larger than expected from traditional thermodynamic calculations [1,2]. While such fluctuations have commonly been considered benign, their presence does preclude the use of Schottky noise as a viable diagnostic and it greatly complicates the implementation of beam cooling techniques [2].
It has long been suspected that such fluctuations can arise from wakefield driven oscillations that are in a delicate equilibrium with dissipative mechanisms. A feasible linear model for such fluctuations in an unbunched coasting beam has been derived in the form of a dielectric function [3]. The equilibrium fluctuation levels for this linear case of interest are given by the fluctuation-dissipation theorem which relates the equilibrium fluctuations to the dielectric function and the Schottky voltage, which in the frequency domain is given by Here, U n is the potential associated with the nth longitudinal mode, i.e. the nth spatial Fourier harmonic of the machine circumference, U s n is the Schottky voltage, which will be derived in a later section, D n ( ) is the well-known dispersion relation for longitudinal modes, which also serves as an equivalent dielectric function for the longitudinal degree of freedom, and G n is the susceptibility. They are given by The energy distribution function, f 0 , is assumed normalized to 1, N is the total number of particles and ω 0 is the synchronous particle revolution frequency. Other quantities above are defined in the next section. The key consequence of this relationship is the fact that fluctuation levels increase strongly as conditions of marginal stability are approached, namely as D n ( ) → 0. We note that as in [3] we will restrict ourselves in this work to the longitudinal degree of freedom of an unbunched beam, which greatly simplifies the analytical procedure. This restriction is not fundamental, however, and our analysis can be extended to other degrees of freedom and to bunched beams. As realized here in the frequency domain, collective effects that are driven by wakefields are represented by an impedance Z n which relates the induced longitudinal voltage at the nth spatial harmonic to its corresponding component of the perturbed beam current. In this case, the machine impedance Z n determines the beam intensity where longitudinal modes become unstable and as the above expression suggests, the fluctuation levels rise sharply as the stability condition is approached, i.e., when the linear growth rate approaches zero. This expression has been confirmed by experimental measurement and is well-accepted in the field of accelerators.
However, as high-energy accelerators are pushed to ever higher intensities, namely closer to marginally stable conditions, the linear response intrinsic in the notion of the dielectric function employed in equation (1) must necessarily break down. The first indication of such effects was already reported a decade ago [4], along with anecdotal evidence that stretches back to the early years of synchrotron development. Namely, as intensity is increased, decidedly nonlinear 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT frequency mixing has been observed in a number of machines. In particular, in the now defunct Main Ring at Fermilab, a series of experiments revealed an orderly cascade of harmonics from a marginally stable resonator impedance. Such phenomena were identified as a form of three-wave mixing and the basic mechanisms outlined in [1]. However the details of this interaction were not studied quantitatively, and the application of the three-wave model to actual devices remains poorly understood.
In the field of plasma physics, an effort has been ongoing for many decades to understand the nature of turbulent states that can occur in systems where collective effects play a significant role [5]. The mathematical formalism grows increasingly complex as the degree of nonlinearity increases, and many such systems can only be effectively modelled numerically. However, in the case of a storage ring, the strict periodicity of the oscillations coupled with the wide separation of eigenfrequencies leads to the possibility that linear modes are weakly coupled through nonlinear interactions and the system behaviour can be determined perturbatively. This is the domain of the weak turbulence theory developed in plasmas and the argument can be made that stored beams satisfy the weak nonlinearity requirement better than many other dynamical systems.
In this work we have developed a comprehensive formulation of weak turbulence in a synchrotron, specifically for longitudinal fluctuations of unbunched beams. It is our purpose here to apply the theory of weak turbulence as formulated in plasma physics to the case of a stored beam. The similarities of the two domains are manifest, but specific aspects of a storage ring, namely its strict periodicity, the strongly decoupled degrees of freedom, and the complex nature of the wakefield or restoring force, lead to unique features in the mathematical treatment. It is worthwhile to note that a series of pioneering papers by Tzenov, e.g. [6], has appeared using related methods to study the nonlinear evolution of a coasting beam, however our approach differs somewhat in both expansion procedure and in our focus on the fluctuation spectrum.
In particular, we are interested in the behaviour of the state of marginal stability where nonlinear effects can be most likely expected to play a role. Many practical devices are commonly run under these conditions. This work extends earlier treatments [3] to include nonlinear interactions and provides a quantitative framework for evaluating experimental observations. It is based on a random phase average of the Vlasov equation including its fully nonlinear properties that leads to coupling among frequencies. The particular feature of a storage ring that facilitates this process is the countably infinite number of eigenfrequencies that are separated by integral multiples of the revolution frequency. This leads to a rich field of coupled modes.
The analysis proceeds by developing second-order expressions for the perturbed distribution function, followed by phase averaging which yields the evolution of the fluctuation spectrum on the slow time scale. A key feature of our work is the inclusion of Schottky noise as the baseline fluctuation spectrum, namely as the thermal reservoir with which wake-driven fluctuations come to equilibrium. This baseline spectrum results as a consequence of the fluctuation-dissipation theorem and we have effectively modified its steady-state form with next-order nonlinear perturbations. In the following sections we will use a Lorentzian equilibrium solution which allows calculation of analytical results. This choice is not a fundamental restriction, but will permit us to find explicit expressions for the nonlinear coupling. We first review the derivation of the fluctuation-dissipation theorem for a storage ring, followed by a formal expansion of the Vlasov equation including nonlinear corrections and the effects of Schottky noise. The resulting equation is introduced into a quadratic form that effectively gives us a separation of fast and slow time scales, and to lowest order yields the fluctuation-dissipation theorem. We then average the equations, giving rise to the so-called kinetic wave equation, effectively a flow equation for the 4 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT photonic fluid that is driven by the wakefields. We then proceed to find steady-state solutions of this set of equations for impedance models relevant to storage rings.

The fluctuation-dissipation theorem for a storage ring
In a stored beam in equilibrium, fluctuations in all macroscopic quantities will arise under the influence of the collective forces, i.e. wakefields that are generated by the beam. In the linear limit, these fluctuations arise from the dielectric response of the beam to small perturbations. Such perturbations, in turn, are generated by the graininess of the beam itself, the so-called Schottky noise of the beam. In the case of fluctuations in equilibrium, a relation can be derived between the fluctuation level and the dielectric function known as the fluctuation-dissipation theorem.
Restricting our attention to longitudinal fluctuations of an unbunched, or coasting, beam, we write the equations of motion in terms of the longitudinal phase-space variables (θ, ), where θ is the position around the ring and = E − E 0 represents the energy deviation from the ring synchronous energy. where I 0 = (eω 0 /2π)N, Z n is the impedance for the nth mode and f n (t, ) is the nth Fourier mode of f . In these canonical longitudinal coordinates for a storage ring, the Vlasov equation is written as also k 0 = −ηω 0 /β 2 0 E 0 is the frequency dispersion, where η represents the machine slip factor, β 0 is the Lorentz beta factor, and E 0 is the synchronous energy. Here we introduce the Schottky voltage U s (θ, t), and its Fourier coefficients U s n (t), as the potential that is induced from the Schottky current through the machine impedance at a given harmonic n. It is manifested as a forcing term in the Vlasov equation (5) above, consistent with a fixed level of shot noise intrinsic to the equilibrium beam distribution. f = f 0 ( ) represents the stationary solution to the Vlasov equation without Schottky noise. In this work, we assume f 0 ( ) in equation (6) below is independent of time. This amounts to foregoing the quasilinear diffusion theory whereby f 0 ( , t) is permitted to vary on a slow time scale, which has been developed previously for a storage ring elsewhere [7]. We shall assume here that such a time scale is long compared to fluctuation equilibration times. However, this assumption is not universally valid and a general model would include both the quasilinear diffusion theory and this fluctuation model solved in a simultaneous manner. The reason for this can be seen by anticipating our final result. Since we base our nonlinear coupling on the three-wave interaction, the rate of change of our fluctuation spectrum is amplitude dependent, scaling as |E n | 2 |E m | 2 , as shown in [5]. Quasilinear diffusion, however, scales as |E n | 2 and for any given beam distribution and machine impedance an amplitude range can be found where the quasilinear theory is relatively insignificant. However, it becomes dominant at sufficiently small fluctuation levels or when the distribution function possesses sharp gradients, hence a general model should include both limits. We will focus only on the threewave coupling process here with the caveat that the quasilinear theory needs to be examined and possibly included as well in any given situation.
To proceed, we expand the distribution function around an equilibrium in the usual way and Fourier analyse in θ: Then, the Vlasov equation becomes Commonly, the Vlasov equation is further analysed at this stage by taking a Laplace transform with respect to t which introduces initial conditions and the associated phenomenon of Landau damping. However, since we are only interested in the long-time behaviour of the system when all transients have died away, we choose to use the equivalent Fourier transform method taking care to ensure that we choose the proper integration contour so that the result convergences to a causal, and finite, result. Namely we evoke the standard practice of adding a small imaginary part to the frequency which suppresses any interactions as t → −∞. With this caveat in mind, we obtain Equation (8) is the basic iteration expression of our method, where we have included λ as an ordering parameter. To derive the linear solution from equation (8), we drop the nonlinear term, multiply by (−eω 0 Z n ) and integrate over to obtain Now the Schottky voltage is essentially the Schottky current times the machine impedance. As shown in [8] the Schottky voltage is given bŷ where θ 0a is a uniformly distributed random angle variable and ω a is another random variable representing the Lorentzian randomly distributed revolution frequency. The ensemble average is to be taken over the index a. Since the Schottky current is randomly distributed around the ring, the expectation value of this current, and hence the potential, is zero. However, if we construct the square modulus of the potential using the above equation and its conjugate, the ensemble average is finite and the result is the fluctuation-dissipation theorem as formulated in equation (1). This result for a storage ring was first described in [3]. It represents the linear response of the beam to the small, random perturbations that occur due to the graininess of the beam itself.
In equilibrium, the beam exhibits fluctuations that carry a signature of the frequency dependence of the machine impedance. For the purposes of our work, this baseline frequency spectrum will form the ground state around which we will further expand the Vlasov equation to include nonlinear mixing of harmonics. As such, this spectrum serves as a reservoir of thermal oscillations with which the full nonlinear fluctuation spectrum must be in equilibrium.

Derivation of the kinetic wave equation
In this section we derive a phase-averaged equation that describes the flow of power among the family of linear eigenmodes of the system. We restrict our analysis to the case of weak turbulence where the nonlinear interaction in the frequency domain can be limited to a set of well-separated eigenmodes. We first develop the nonlinear coupling equations from the Vlasov equation by implementing an iteration process [9], then perform a perturbation expansion of the resulting system, treating the mode amplitudes as statistically independent quantities. This leads to a kinetic wave equation that describes the differential flow of power across the frequency spectrum as the system approaches steady state. For the purposes of this work we are only interested in the equilibration of fluctuation levels that is accomplished on a phase randomization time, typically much shorter than changes in the distribution function can occur. We will begin by carrying out two steps of the iteration process as described in [9] using equation (8). We rewrite (8) by replacing the dummy variables m and by l and . We then replace n by n − m and by − and differentiate with respect to to obtain Substituting equation (10) into the right-hand side of equation (8) we have the following result 7 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT We apply the iteration process one more time because we wish to include the quadruple potential that represents the most general form of the three wave interaction. Since we only keep terms through O(λ 3 ), we simply replacef k ( ,˜ ) on the right-hand side by −i 2π We then multiply by (−eω 0 Z n ) and integrate over to obtain where we have dropped O(λ 3 ) terms.
The coefficient V( , − ) is the wave-wave coupling coefficient and W( , − , − − ) is the wave-particle-wave coupling coefficient known as induced scattering [5]. They are defined by where It can be shown that the V and W coefficients obey the following symmetry: In what follows, we take where we have chosen a Lorentzian distribution to facilitate the analysis that follows, while still including Landau damping by virtue of the energy spread intrinsic in f 0 . The dispersion relation is then

Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
This has a pole in the lower half plane and its zeros are given by n = nω 0 + ω n sgn(nk 0 ) nk 0 α, We note that for small enough beam current the zeros are in the lower half plane. In order to separate the fast and slow time scales, we introduce the quadratic expression which we will evaluate separately using both sides of equation (12). The form of this expression is chosen so that the fluctuation-dissipation theorem results when nonlinear terms are neglected. Moreover, this particular asymmetric form is chosen so that the leading term is of the form D n |U n | 2 , which will permit us to further expand D n ( ) to extract the slow time scale. Each factor individually represents the deviation from the fluctuation-dissipation theorem. This implies we are expanding about the Schottky spectrum as the equilibrium solution. We will also average over the fast time scale on both sides of equation (12) after applying the quadratic form. This will give us the desired form for the kinetic wave equation.
On the left-hand side of (12), we will first expand the dispersion relation in a Taylor expansion around an eigenfrequency, n , which introduces the slow-time scale variation of the mode amplitudes: where the prime indicates a derivative with respect to frequency. Note that D n ( n ) = D * n ( * n ) = 0. By taking the expectation value of the Schottky voltage U s n , the cross terms vanish in expression (17), which then becomes In the above equation, F( n ) represents the inverse Fourier transform of the square of the ensemble-averaged Schottky fluctuations, namely 9 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT make use of the following symmetry relations the imaginary part of equation (18) can be rewritten as We proceed to solve for the nonlinear corrections to this expression by means of a perturbation expansion, assuming sufficiently weak coupling between the modes. Anticipating that the total potential is a sum of the response to Schottky excitation and independently driven fluctuations at each of the eigenfrequencies, we expand the potential aŝ The goal at this stage is to find an equation which governs the amplitudes of the phaseaveraged potentials U (1) n . To this end we average the resulting equations by invoking random phases among the eigenmodes. This assumption is justified provided we focus on time scales longer than the inverse frequency width of the beam. Our assumption of randomly phased eigenmodes is invoked by considering each potential U n to be a statistically independent variable with each phase φ(n) randomly distributed in (0, 2π). Thus, when taking an ensemble average of terms with odd powers of U s n , we obtain for example We also assume n n such that n ≈ * n , which is true for marginally stable, well separated eigenmodes, and we will use the symmetries U −m = U * m and −m = − m which force the double sums from equation (12) to collapse to only two terms.
The full kinetic wave equation is thus given as d dt The function F( n ) can be evaluated explicitly by performing the contour integral in equation (19). In the case far from the stability boundary, the pole at nω 0 ± ink 0 α dominates and we obtain an expression for the Schottky power that is peaked at nω 0 , similar to that in [8]. But for marginally stable conditions, namely when i n ≈ 0, the pole closest to the real axis dominates and the Schottky power is shifted to the eigenfrequency, as shown in [3]. It is worthwhile to note that the frequency domain version of this coefficient differs from the time domain expression because we have moved one factor of D n ( ) to the left-hand side of equation (1) before inverting the Fourier transform. The expansion of the dispersion relation around a marginally stable eigenfrequency then permits us to extract the slow time variation of the fluctuation spectrum required for the kinetic wave equation.
For our purposes, we are mainly interested in solving the steady-state equation given by the time-independent part The complex coefficients V 1 , V 2 , V 3 , W in equations (21) and (22) above are calculated using the residue theorem and may be explicitly evaluated owing to our choice of f 0 ( ). Details of these calculations are given in the appendix.

Model results
For purposes of illustration we consider the case of a resonator impedance with a resonance located in the high mode number region. Moreover, we adopt a low-Q value that is consistent with the impedance commonly associated with the beampipe cutoff frequency, i.e. where the peak resonant impedance is chosen so that the growth rate is near zero at resonance. The resonator impedance has the following expression where R represents the resonator cavity's resistance, Q is the quality factor, and ω R is the resonance frequency. Such a choice facilitates the coupling across a band of eigenfrequencies, which will illustrate the physics. However, as we will discuss later, the non-smooth distribution functions often found in hadron storage rings can lead to similar results. In this case, we have found that the coefficient V 2 is several orders of magnitude larger than the rest of the coefficients in equation (22) as shown in figures 1-4. This situation results from the fact that this coefficient exhibits resonant behaviour as the eigenfrequencies approach integral multiples of each other, and specifically when the imaginary part of these frequencies is such that the pole in the expression for V , equation (13), occurs in the upper half plane, as is shown in the appendix. This implies that coupling of two modes at high n with a mode at lower n is the dominant process. A useful simplification, therefore, is to keep only the V 2 term in the nonlinear coupling. We may then readily find the steady-state solution by an iterative solution  of the nonlinear equation set where we use the Schottky voltage to determine the initial guess according to equation (1). A more generally relevant physical case occurs when there is a critical impedance, i.e., marginally stable conditions, at multiple frequencies. Such a situation might occur, for instance, when a resonator contributes an impedance at high mode number superimposed on the resistive wall impedance which is dominant at low mode number. In the above, Z 0 ≈ 377 is the impedance of vacuum, b is the beam pipe radius, µ 0 is the permeability of vacuum, 1/σ c is the material resistivity, and the square root expression represents the skin depth.   This term represents the wave-particle-wave interaction, or induced scattering, and is formally of higher order than the three-wave interaction depicted in the above coefficients. The largest value in this plot is significantly smaller than other coefficients consistent with the ordering assumed in this work.
The results for this model impedance are shown in figures 5 and 6 where we show the spectrum for a case where both high and low mode numbers have marginally stable bands. The enhancement of the effective Schottky spectrum can be significant, as shown, which is due to the effects of three-wave coupling. We note that there is also a slight reduction of the threshold for linear stability for this system through the interaction of the nonlinear terms.

Discussion and conclusions
It is evident from the form of the kinetic wave equation (21) that marginally stable systems permit abundant mixing of frequencies in the steady state, namely on a time scale that is long with respect to the inverse frequency width of the beam. For most synchrotrons, the stored beam lifetime is many orders of magnitude larger than this time scale, which justifies both our expansion procedure, as well as the steady-state solution. The details of the resulting fluctuation spectrum are sensitively dependent on both the frequency dependence of the machine impedance and on The deviation from the Schottky spectrum occurs at the lowest mode number which is excited by three-wave interactions. The ordinate is evaluated for a specific beam power monitor, but is essentially a relative scale. The enhancement above the Schottky spectrum occurs as a result of multiple three-wave interactions described by the kinetic wave equation. It is interesting to note that jumps in the spectrum occur at locations where specific mode couplings onset. The ordinate is evaluated for a specific beam power monitor, but is essentially a relative scale. the details of the energy distribution f 0 . For instance, we find that the nonlinear susceptibility depends on the product of the impedances at the eigenfrequencies in question and also scales as α −4 , to lowest order, where α is the energy width of the beam distribution. This implies that small deviations in the energy distribution can give rise to nonlinear effects even at modest impedances and beam intensities.
For instance, in the proton storage ring (PSR), an energy distribution with a sharp notch has been observed [10]. The complexity of this distribution was attributed to high-order transverse behaviour that, through the machine chromaticity, caused specific particle momenta to be lost at transversely resonant tunes. Similar behaviour has been seen in other hadron synchrotrons [1]. Once such a distribution is created, each of the notches in the distribution function creates large nonlinear susceptibilities that can give rise to broadband frequency coupling. This model can, in principle, qualitatively explain the observations, though without detailed knowledge of the energy distribution, quantitative comparison is difficult.
It is worthwhile to remind the reader that quasilinear diffusion, namely the change in the beam distribution itself due to the presence of fluctuations, is likely to occur simultaneously with the nonlinear wave-wave interactions treated in this work. A general solution to the problem of saturated beam states would then include a simultaneous solution of both the Kinetic Wave Equation and the quasilinear equation for f 0 ( , t), as given in [7].
We note also that some of the terms in the nonlinear expansion are always negligible, such as V 1 and W . These terms represent higher-order interactions that require specific impedances that are not found among the types we have studied. Finding solutions of the kinetic wave equation in synchrotrons that have complex impedance functions may not be possible in general because of the lack of detailed information. However, it is conceivable that in machines with sufficiently regular behaviour that quite predictable spectra can be determined, and hence used to determine machine parameters.
Finally, we would like to comment how this model may be extended to other situations, including the transverse degree of freedom or bunched beams. In principle, there is no restriction to the application of this theory provided the associated linear eigenmodes are well-separated and satisfy the three-wave coupling conditions. As such, even complex situations involving coupled transverse and longitudinal oscillations (e.g. synchro-betatron coupling) may be handled.
The key issue regarding the applicability of this model to any of these cases is whether a sufficient multiplicity of modes exists that will randomize the mode phases on a sufficiently short time scale.

Appendix. Discussion of coupling coefficients
We present here an exposition on how the coupling coefficients V 1 , V 2 , V 3 and W in equation (22) are evaluated. The choice of a Lorentzian distribution permits us to find exact analytical expressions for these coefficients.
Coefficients V 1 and W represent the coupling strengths between modes m, n and n − m. They can be written as To evaluate the integral in the expression above we will use residue theory. We will choose the contour of integration in such a way that we approach the Landau pole from a direction that is in agreement with causality. The integrand above is a rational function in the complex plane and has a pole of order one at q 1 = n /nk 0 α, a pole of order two at q 2 = ( n − m )/(n − m)k 0 α, and poles of order three at ± i. Two cases arise in the evaluation of the integral: the case when all the poles are distinct and the case when the first two poles enumerated coalesce. The results for these two cases follow.