This work was motivated by a novel application in particle physics, namely the dynamics of a generic circular particle accelerator. Since this application is expected to be unfamiliar to most readers of this journal, we offer the following description of a synchrotron [7"Cornell Electron Storage Ring", CLASSE: CESR. Cornell laboratory for accelerator-based sciences and education, . Available from: http://www.lepp.cornell.edu/Research/CESR/WebHome.html].

The synchrotron is a particle accelerator in which a “Particle” actually consists of a group of electrons called a “Bunch,” and the collection of all bunches is called a “Train.” We ignore the interactions of electrons inside each bunch and treat the entire bunch as a single particle.

Each bunch leaves an electrical disturbance behind it as it traverses around the synchrotron, and these wake fields are the main source of coupling in the model. The coupling is mediated by several sources, including ion coupling and the electron cloud effect.

Particle paths in the synchrotron are circle-like but are not exact circles. Since the synchrotron lacks a central force, the circle-like particle orbits are achieved through the use of about 100 electromagnets spread around the periphery Fig. (1).

Fig. (1) Two bunches moving clockwise along a polygonal path through the use of a system of electromagnets.

This means that the magnetic external forcing is periodic in rotation angle θ; assuming that the angular velocity of the bunch is constant with θ = ωt, the forcing is periodic in time as well. We can express this forcing function as a Fourier series, and we shall approximate this series by the first couple of terms in it, namely the constant term and the first cosine term.

We model each bunch as a scalar variable x_i(t), i = 1,...,n. Here, x_i is the vertical displacement above equilibrium of the center of mass of the i^th bunch. Each x_i is modeled as a damped parametrically-forced oscillator, and we write:

The nonlinear terms are included to create a more realistic model since most natural phenomena are nonlinear and linear models are a convenient approximation. The nonlinear parameter, γ, can be chosen to adjust the scale of the problem.

The damping terms are also included to create a more realistic model.

The coupling terms on the right-hand side may be viewed as the strength of the plasma interactions: the radiation from a bunch produces an electron cloud which travels behind the bunch. This radiation dissipates away very quickly though, and so this coupling has a short range and can only influence the dynamics of the next bunch in the train. The coupling strength is affected by both the spacing between bunches as well as the charge of each bunch, and α encapsulates both of those effects.

The dynamics of the first bunch are given by eqs. (8, 9).

As this system is 2-dimensional, the complete dynamics can be expressed in the phase plane A₁-B₁. The Numerical Results section contains various graphs demonstrating the full dynamics of the first bunch, but the rest of this section will focus purely on calculating the equilibrium points of the system.

Setting these equations equal to zero gives us the following equilibrium points:

Here, we use the standard notation to denote equilibrium points of the variables A_i, B_i.

To simplify matters, we will focus on the special case when δ₁ = 0 and μ = 0. Setting δ₁ = 0 lets us examine the point of resonance directly without considering points in the neighborhood of the resonance, and setting μ = 0 ignores the effects of damping which ends up having a negligible effect on the analysis of the steady state. The fixed points in this case are:

(12)

Here, we note our first major observation: Nontrivial real equilibrium points for A₁ only exist for γ < 0, and nontrivial real equilibrium points for B₁ only exist for γ > 0. However, the magnitude of the amplitude is the same in both cases, and the analysis of both cases is identical. Without the loss of generality, we will take γ > 0 and = 0.

The dynamics of the second bunch are given when i = 2 in eqs. (10, 11).

(13)

(14)

Note that these equations share a similar structure to eqs. (8 and 9) but with additional terms resulting from the one-way coupling. These additional terms mean the dynamics of the second bunch are 4-dimensional instead of 2-dimensional like the first bunch, and we cannot view the dynamics in a phase plane. Our analysis of the second bunch will focus only on the steady-state solutions and not on the general dynamics.

Since we are taking γ > 0, all equilibrium points for the first bunch require = 0; substituting = 0 into eqs. (13, 14) yields:

(15)

(16)

Note that + is a nonnegative value, and is only identically zero in the trivial solution (, ) = (0,0). Since we are interested in nontrivial solutions, we require that + be positive.

However, since γ > 0, the expression 3γ( + ) + 2 must be positive. Thus, the only way for (16) to equal zero is for = 0.

Substituting = 0 into (15) yields:

(17)

Since is already known to us from eq. (12), solving this cubic equation will give us the solution for .

This solution gives us the steady-state of B₂. The full dynamics of the system is 4-dimensional, involving A₁, B₁, A₂ and B₂, and cannot be expressed in a phase plane.

A graph of the relationship between and α can be seen in Fig. (2).

The dynamics of the third bunch are given when i = 3 in eqs. (10, 11).

(18)

(19)

Note that these equations share the exact same structure as eqs. (13, 14). Just like with the second bunch, the presence of additional terms means that the dynamics of the third bunch are 4-dimensional. This means that we can’t view the dynamics in a phase plane, and so our analysis of the third bunch will focus only on the steady-state solutions and not on the general dynamics.

Fig. (2) Plot of as a function of α, for γ = 1 and .

We know from the analysis of the second bunch that = 0 is necessary for all steady state solutions of the second bunch; substituting = 0 into eqs. (18, 19) yields:

(20)

(21)

Note that + is a nonnegative value, and is only identically zero in the trivial solution ( + ) = (0,0). Since, we are interested in nontrivial solutions, we require that + be positive.

However, since γ > 0, the expression 3γ( + ) + 2 must be positive. Thus, the only way for (21) to equal it to zero is for .

Substituting into (20) yields:

(22)

Since is already known to us from eq. (17), solving this cubic equation will give us the solution for .

The graph of the relationship between and α will look qualitatively similar to the relationship between and α Fig. (2) but they will differ quantitatively since the value of that is substituted into eq. (22) will in general be different than the value of that is substituted into eq. (17).

Still, given the similarity between eq. (22 and 17), it is natural to ask if this pattern continues for all later bunches. Indeed, this is the case, and so we will generalize the results found in these past two sections to the n^th bunch in the system, where n can be any integer.

Due to the nature of one-way coupling, the dynamics of all the bunches except the first one are identical, as can be seen in the results found for the second and third bunches. Thus, we can easily derive a formula for calculating the steady-state solution for the n^th bunch for any n > 1:

(23)

Unfortunately, calculating in practice requires calculating the amplitudes of all bunches in front of it, since the recursive relationship cannot be simplified into a formula dependent only on the first bunch. Part of the problem is that each step requires solving a cubic equation which, although solvable in principle, is a tangled mess. In practice, it is much easier to use numerical root solving methods to find the amplitudes for all the bunches.

The other, bigger problem is that cubic equations have three roots; if all three roots are real and distinct, then we have a multi-valued function. Fortunately, this only occurs for a range of α values, and we can easily determine this region through analytic means.

Note that in Fig. (2), there are two points with infinite slope: these are the points that divide the function into multi-valued regions and single-valued regions. If we flipped the two axes, then these two points change from having infinite slope to zero slopes see Fig. (3). Thus, we want to derive and find the values of α for which the derivative is equal to zero.

Fig. (3) Plot of α as a function of , for γ = 1 and .

We start by moving the α term in eq. (23) to the right hand side and differentiating both sides.

Here is a constant since we’re assuming all up to i = n - 1 have already been found.

Solving for , we obtain:

We set this equal to zero to find the local extrema:

Finally, we substitute this result into eq. (23) to find the corresponding α value:

(24)

One important feature of this result is that the range of α values for the n^th bunch depends on the bunch before it. As an example, for , the α range depends on ; multiple limit cycles are possible for when α[-3^-3/2,3^-3/2] ≈ (-0.2,0.2).

The significance of this result is that fixing α does not guarantee that all bunches will be either multi-valued or single valued; it is possible to choose α such that only has one limit cycle, but has three.

A natural question to ask at this point is whether there exist values of α that do guarantee single valuedness for all . Due to the relationship of eq. (24), a uniform bound on α requires a uniform bound on , so this question is equivalent to asking if there is a bound on how large can become.

It turns out that the are indeed bound by an upper bound, and this can be shown by examining the limit as n → ∞.

In particular, we care about the limit of the amplitude ||. One way for this limit to exist is if B* = = as n → ∞. In this case, each bunch has the same amplitude as the bunch before it, and each bunch is oscillating in phase with the bunch before it.

By setting B* = = in eq. (23), we obtain:

(25)

Thus, we find that the possible in-phase limits are:

(26)

As long as α ≥ - 1/2, these three limits are distinct; otherwise, no in-phase limit cycles are possible.

Another possibility is to examine the case when B* = = as n → ∞. In this case, each bunch has the same amplitude as the bunch before it, and each bunch is oscillating 180 degrees out of phase with the bunch before it.

By setting B* = = - in eq. (23), we obtain:

(27)

Thus, we find that the possible limits are:

(28)

As long as α ≤ 1/2, these three limits are distinct; otherwise, no out-of-phase limit cycles are possible.

Therefore, in the limit as n → ∞, there are three possible cases:

• For α < - 1/2, only out-of-phase limit cycles can exist.
• For -1/2 < α < 1/2, both out-of-phase and in-phase limit cycles can exist.
• For 1/2 < α, only in-phase limit cycles can exist.
  • Fig. (4) shows the upper bound for both types of limit cycles.

Fig. (4) The solid curve is the upper bound on amplitude of the in-phase limit cycle. The dashed-dot curve is the upper bound on the amplitude of the out-of-phase limit cycle.

To help visualize the dynamics of the first bunch, graphs showing the phase plane for eqs. (8 and 9) are provided below. All phase plane graphs were made using PPLANE [10 MATLAB v.7.5 by Mathworks; applet pplane7.m by Polking, J. C , Rice University, Department of Mathematics, Available from: http://math.rice.edu/∼dfield/dfpp.html].

Fig. (5) compares the cases when γ > 0 and γ < 0. Qualitatively, the graphs are the same; the only difference is that one has the equilibrium points on the A axis and the other has the equilibrium points on the B axis.

Fig. (6) shows the effect of including small values of δ₁ and μ. The inclusion of damping has broken the homoclinic orbit and all orbits are eventually attracted to one of the two stable equilibrium points. Since the two basins of attraction are intertwined, it can be hard in practice to determine which equilibrium point will be reached from a given initial condition.

Fig. (7) shows the effect of including larger values of δ₁ and μ. The equilibrium points have gone through a pitchfork bifurcation and there is now only one equilibrium point: the origin.

Since all bunches other than the first bunch are coupled to another bunch, we cannot express their dynamics in a phase plane. Instead, we will show the amplitudes through in a cobweb diagram for a fixed number n.

Each diagram contains the graphs of Eq. (23), with on the x-axis and on the y-axis, and the line y = x. The sequence starts with on the x-axis, and proceeds as follows:

• Move vertically to the line y = x.
• Move horizontally to the curve given by Eq. (23) [may be multi-valued].

Fig. (5) Phase plots of the A1-B1 dynamics. The left plot shows the dynamics for γ = 1 and the right plot shows the dynamics for γ = -1. Both plots have α = 0 and μ = 0.

Fig. (6) Phase plot of the A1-B1 dynamics. Parameter values are α = 0.1, μ = 0.1 and γ = 1.

Fig. (7) Phase plot of the A1-B1 dynamics. Parameter values are α = 1, μ = 1 and γ = 1.

This process is repeated n times, with each horizontal change determining the next value of .

For example, Fig. (8) shows the case when α > 1/2. The diagram begins at the point (0.8,0), representing the value . After moving up to the line y = x at (0.8,0.8), the process then moves to the right to (1.1,0.8), where 1.1 represents the next value, . After several iterations, the values of approach the point (1.2,1.2), where 1.2 is the size of the stable limit cycle in this region.

Fig. (8) Cobweb diagram for α = 0.6 and γ = 1.

Fig. (9) shows the case when α < -1/2. The diagram begins at the point (0.8,0), as before. This time, after moving up to the point (0.8,0.8), the process then moves to the left to (-1.2,0.8), where -1.2 represents the next value, . After several iterations, the values of alternate between two values, 1.3 and -1.3, reflecting the out-of-phase nature of the stable limit cycle in this region.

Fig. (9) Cobweb diagram for α = -0.6 and γ = 1.

For -1/2 < α < 1/2, there are three possible limit cycles. Since the output of eq. (23) is multi-valued in this region, there are many different cobweb diagrams for a given starting point , with the realized outcome determined by the initial condition in the original system (1, 2). The code used to generate the diagrams picks one of the three limit cycles at random, as this best represents the unpredictability of knowing the precise initial condition.

Fig. (10) shows the chaotic nature of this iteration map in the region -1/2 < α < 1/2. Since the behavior here is randomized, there is no clear pattern to be discerned here. However, we note that the attractor seems to be a fractal of some kind since there are gaps that are never reached.

Fig. (10) Cobweb diagram for α = 0.1 and γ = 1.

Fig. (11) Plots of vs α for α ≥ 0.3 and γ = 1. The top left graph shows n = 2, the top right graph shows n = 3, the bottom left graph shows n = 4, and the bottom right graph shows n = 5.

The most important question in this model is determining how long a train of bunches can be before the tail becomes unstable. We have shown that there is a theoretical upper bound to the limit cycle amplitude, but it remains to be seen how quickly the numerical sequence given by Eq. (23) approaches this limit.

Fig. (11) shows a sequence of graphs of vs α for 2 ≤ n ≤ 5 and α ≥ 0.3, which demonstrates the speed at which this sequence approaches the limit. However, we note that the sequence does not seem to converge pointwise at the same rate, and for larger α values, it takes longer to converge to the limit.

Fig. (12) shows the graph of vs α for a larger range of α values, including the multi-valued region. While the multi-valued region does not show convergence to either of the limit curves, the values do stay bounded between the maximum values of the two limit curves. Thus, even in this region we, can place an upper bound on how large the can grow.

Fig. (12) Plots of vs α for -2≤α≤2 and γ = 1.

Fig. (13) shows the graph of vs α. For α > 0.5, this figure is almost identical to Fig. (12), but for α < -0.5 the sign of the curve is now negative instead of positive. This reflects how the limit cycle in this region is out-of-phase, as is out of phase with . Since the magnitude of the limit curve is the same for both positive and negative branches, the sign of is of little concern.

Fig. (13) Plots of vs α for -2≤α≤2 and γ = 1.