Spin motion at and near orbital resonance in storage rings with Siberian Snakes. Part I: at orbital resonance

Here, and in a sequel, we invoke the invariant spin field to provide an in--depth study of spin motion at and near low order orbital resonances in a simple model for the effects of vertical betatron motion in a storage ring with Siberian Snakes. This leads to a clear understanding, within the model, of the behaviour of the beam polarisation at and near so--called snake resonances in proton storage rings.


Introduction
In earlier papers we and collaborators have emphasised the utility of the invariant spin field (ISF) and the amplitude dependent spin tune (ADST) for analysing spin motion in circular particle accelerators and storage rings [1,2,3,4,5,6,7,8,9,10]. In particular, under certain conditions, the ISF is unique up to a global sign and in that case it allows estimates to be made of the maximum equilibrium beam polarisation and the maximum time averaged beam polarisation in proton storage rings. Then, for example, for a given equilibrium distribution of particles in phase space, the maximum attainable polarisation at the chosen high energy can be estimated before embarking on extensive computer simulations of the effect on the polarisation of acceleration from low energy. Once a machine configuration has been found which appears to be acceptable at the chosen high energy, one then studies the effect of acceleration to assess whether the configuration is still acceptable. Acceleration can involve crossing many spin-orbit resonances and that can lead to a loss of polarisation. The latter problem can be partially solved by the inclusion in the ring of so-called Siberian Snakes [11,12], magnetic field configurations that cause the average spin precession rate on the design orbit to be independent of the nominal beam energy. Nevertheless, spin-orbit resonances can still occur but their identification then often requires a more careful definition of the spin precession rate than has been common among practitioners, involving the amplitude dependent spin tune. A full understanding also requires a careful definition of an adiabatic invariant for spin motion. In most of the numerical investigations described in [1,2,3,4,5,6,7,9,10], orbital resonance is avoided. Moreover, it is shown that the spin-orbit systems tend to avoid exact spin-orbit resonance. These and other matters are explained and illustrated in great detail in the sources cited above. In order to keep this paper to a reasonable length we will assume that the reader is familiar with that material.
Of course, the ISF and the concepts derived from it, may be of little help if the ISF is not unique. That can be the case if the orbital motion is resonant or if the system is on spin-orbit resonance [5,7]. Nevertheless, as we show below, special choices from sets of non-unique ISF's can be useful for investigating spin motion near some kinds of orbital resonance. Moreover, the ISF is still useful at rational vertical orbital tunes corresponding to the so-called odd order snake "resonances". At these tunes the Siberian Snakes apparently do not succeed in preventing loss of polarisation during acceleration [13,14,15,16,17,18]. However, with the exceptions of [19,20,5,21], discussions about spin motion at or near to these tunes have made no reference to the ISF. The treatment in [19,20] involved a mathematical approximation to the model used in this paper. Then in [5] it was pointed out for the first time that at these tunes the ISF is an irreducibly discontinuous function of the vertical orbital phase and that the discontinuities can be moved, thereby demonstrating non-uniqueness. In [21] the necessity of the discontinuities was disputed (see Section 3.4). In Section 2 we explain that exactly at these special tunes, the term "snake resonance" does not fit with our preferred definition of spin-orbit resonance. Nevertheless, for simplicity, we adopt the now traditional nomenclature. In [5] it was also made clear how non-uniqueness can occur at other rational tunes.
In this paper and in a sequel (called Part II) we extend the investigations in [5]. In the initial and pioneering work on snake resonances in [13,14,15], emphasis was placed on the significance of the so-called "perturbed spin tune", a measure of the angles of spin rotation around the real, unit length, eigenvectors of 1-turn SO(3) spin maps. See also [22]. However, these eigenvectors are usually not solutions of the Thomas-Bargmann-Michel-Telegdi (T-BMT) equation along the trajectories. Thus, while it is clear from calculations that the "perturbed spin tune" can show strong variations, we do not consider its behaviour to be relevant to the discussion [8].
In [13,14,15] spin motion was also analysed in terms of an essentially perturbative expansion of the p-turn SU(2) spin transfer matrix, T (p), and it was found that at snake-resonance tunes, |T 21 (p)| could increase without limit as the number of turns p increased. In so far as it relates to positions in tune space, this behaviour, which is an artifact of the perturbative approach, appears to be consistent with the snake resonance phenomenon. However, although an unlimited increase of a matrix element in a perturbative expression for a rotation matrix does suggest exceptional behaviour, it destroys the unitarity of the matrix, thereby demonstrating an invalid approximation and implying a consequent limitation of the predictive power of the calculation. For example, in the absence of other input, one might suppose that an unlimited growth of |T 21 | could infer that initially vertical spins are simply flipped. Alternatively, the growth might be a hint that the vertical component of the beam polarisation oscillates as spins rotate around a horizontal axis. Lastly, simulations reported in [10,23,24] demonstrate the effects of varying the rate of acceleration near snake-resonance tunes. The number of turns needed to traverse a given energy range depends on the energy gain per turn. Then, if weight is given to the perturbative treatment, the number of turns determines how large |T 21 | can become. The rate of acceleration is certainly important in the Froissart-Stora calculation [25] of the loss of polarisation when crossing spin-orbit resonances in rings without snakes, and the phenomenology is well understood. However, the simulations in [10,23,24] show that the dependence of the final polarisation on the acceleration rate can be complex and unexpected and that no clear picture emerges.
To summarise, in our opinion, although snake resonances have presented problems [17,18], the numerical and theoretical investigations made so far have provided no completely coherent picture of spin motion at and near snake-resonance tunes, either with or without acceleration.
These papers provide a new contribution towards such a picture, at least within our adopted simple model. We carry out our study against the background of our standard philosophy, namely that to detect exceptional behaviour, one should start spin-orbit tracking simulations with an equilibrium distribution of particles in phase space and with each spin parallel to the ISF vector corresponding to the position of the particle in phase space [7]. Then any unexpected behaviour is signalled by long term or turn-to-turn variations of the polarisation of the beam. This gives a much cleaner view of the situation than if one just begins in the common way with spins parallel to the direction of the ISF on the closed orbit. Accordingly, with the ISF at the centre of our discussion we show how, in the cases considered, the long term behaviour of spins can be inferred, at least qualitatively, from some features of the ISF. On low order orbital resonance, an ISF can be calculated almost trivially from the spin maps of a few turns.
For our purposes, and in order to allow direct comparison, it suffices just to consider a model used in earlier literature [13,14,15], namely a model with two Siberian Snakes. Since the ranges of the relevant parameters and the number of possible configurations is huge, this study, which is mainly numerical, is not exhaustive. We fully appreciate that storage rings do not run on low order orbital resonance, that spin-orbit resonances need not be well separated, that particles have three modes of oscillation and that particle motion in real rings can be nonintegrable. Nevertheless our study provides useful insights.
The paper is structured as follows. We continue in Section 2 by recalling the simple idealised and traditional model of spin motion for protons considered in [5,13,14,15] and specify the notation commonly used to describe it. Then in Section 3 we use the model to study spin motion exactly at orbital resonances including an odd order snake resonance and show how the chief features of spin motion can be guessed from the characteristics of the ISF. We summarise our studies in Section 4. Part II of this study completes the picture by addressing spin motion close to, but not at, an odd order snake-resonance tune. The numerical calculations were carried out with purpose-built spin-orbit tracking codes, with the spin-orbit tracking facilities in the code SPRINT [3,4] and with the SODOM-II algorithm [26] embedded in SPRINT.

Recapitulation -the single resonance model with two snakes
Spin motion in the electric and magnetic fields at the point z in the 6-dimensional phase space at beam energy E 0 and at the position s around the ring, is described by the T-BMT precession equation d S/ds = Ω( z; s, E 0 ) × S [27,28,1] where S is the spin expectation value ("the spin") in the rest frame of the particle and Ω( z; s, E 0 ) contains the electric and magnetic fields in the laboratory and depends on the beam energy E 0 . The ISF, whose value at ( z; s) is denoted byn( z; s), is a 3-vector field of unit length obeying the T-BMT equation along particle trajectories ( z(s); s) and fulfilling the periodicity conditionn( z; s + C) =n( z; s) where C is the circumference 1 . Thusn( M( z; s); s + C) =n( M ( z; s); s) = R 3×3 ( z; s)n( z; s) where M ( z; s) is the new position in phase space after one turn starting at z and s, and R 3×3 ( z; s) is the corresponding spin transfer matrix. For convenience we have suppressed the dependence of M , R andn on E 0 . In addition to the kinematical constraint |n| = 1, a complete definition of the ISF requires the specification of a constraint on its regularity with respect to z. For example, one could require thatn( z; s) is continuous in z. It is clear that such regularity conditions are needed since, for example, a piece-wise continuous ISF exists if a continuous one exists but not vice versa. See Section 3.4 and [5]. However, since the emphasis of the paper is on numerical results, we only occasionally dwell on the matter of regularity. We use the term "global uniqueness" if two ISF's can differ only by a sign. Thus in the case of global uniqueness, either exactly two ISF's, ±n, exist as in Section 3.1 or none, as in Section 3.4. We use the term "local uniqueness" if any two ISF's,n andn ′ are parallel, i.e.n ×n ′ = 0, so thatn andn ′ can differ only by a sign function. Of course global uniqueness implies local uniqueness but not vice versa. Since the issue of local uniqueness is beyond the scope of this paper, it will be addressed only briefly. If an ISF exists and parameters such as E 0 are constant, the scalar product J s = S ·n/| S| is invariant along a trajectory. For a turn-to-turn invariant particle distribution in phase space, a distribution of spins initially aligned along the ISF remains invariant from turn-to-turn, i.e., in "equilibrium". Moreover, for integrable orbital motion and away from both orbital resonances and spin-orbit resonances (see below), the average | n( z; s) | ofn over the phases on a torus is the maximum attainable time averaged beam polarisation P lim . Away from orbital resonances and spinorbit resonances the actual time averaged polarisation can be written as P lim P dyn where the P dyn = | J s | depends on the history of the beam [4]. For a turn-to-turn invariant particle distribution in phase space P lim = | n( z; s) | is also the maximum attainable equilibrium beam polarisation. This is reached when P dyn = 1.
Under appropriate conditions J s is an adiabatic invariant while system parameters such as the beam energy E 0 are slowly varied [3,9]. In factn then serves as a "template" for spin motion. Several examples of this are given in Section 3.
The ADST ν s ( J) at the amplitudes (actions) J, is the number of spin precessions around then per turn on a trajectory, viewed in a so-called uniform precession frame (UPF). See [7] for precise definitions for smooth systems, i.e., systems with continuously differentiable functions, and for an explanation of how a particular ADST is, in fact, a member of an equivalence class. Note that although the systems in this paper are not smooth in s due to the presence of point-like snakes (see below), their smoothness in z facilitates a close analogy with the smooth systems of [7].
In general, an ADST does not exist if the trajectory is on orbital resonance but on the other hand, one avoids running a machine on orbital resonances, at least those of low order. If an ADST exists, it depends only on J, hence the name ADST.
The ADST provides a way to quantify the degree of coherence between the spin and orbital motion and thereby predict how strongly the electric and magnetic fields along particle trajectories disturb spins. In particular, the spin motion can become very erratic close to the spin-orbit resonance condition ν s ( J) = k 0 + k 1 Q 1 + k 2 Q 2 + k 3 Q 3 where the Q's are orbital tunes and the k's are integers. Near these resonances the ISF can spread out so that P lim is very small. The spin tune on the design orbit ν 0 ≡ ν s ( 0) always exists and so doesn 0 (s) ≡n( 0; s).
In this paper we shall be concerned mainly with those orbital resonances where the Q's are rational. We write the fractional parts, [Q i ], of rational tunes Q i (i = 1, 2, 3) as a i /b i where the a i and b i are integers. Here and later the brackets [...] are used to signal the fractional part of a number. For rational [Q i ] a trajectory is periodic over c turns where c is the lowest common multiple of the b i . This opens the possibility that in this case the ISF at each ( z; s) can be obtained (up to a sign) as the unit length real eigenvector of the 3 × 3 orthogonal matrix representing the c-turn spin map (c.f. the calculation ofn 0 from the 1-turn spin map on the closed orbit). However, the corresponding eigentune cν c extracted from the complex eigenvalues Λ c = e ±2πicνc , depends in general on the synchrobetatron phases at the starting z. Thus in general ν c cannot be used to find a spin tune. Nevertheless if c is very large the dependence of ν c on the phases can be very weak so that it can approximate well the ADST of nearby irrational tunes. For non-resonant orbital tunes, the spin tune can be obtained using the SODOM-II algorithm [26] or from averaging the pseudo spin tune [3,4].
In perfectly aligned flat rings with no solenoids,n 0 is vertical and ν 0 can be chosen to be aγ 0 where γ 0 is the Lorentz factor on the closed orbit and a is the gyromagnetic anomaly of the particle. In the absence of skew quadrupoles, the primary disturbance to spin is then from the radial magnetic fields along vertical betatron trajectories. The disturbance can be very strong and the beam polarisation can be small near the condition aγ 0 = κ ≡ k 0 ± Q 2 where k 0 is an integer and mode 2 is vertical motion. This can be understood in terms of the "single resonance model" (SRM) whereby a rotating wave approximation is made in which the contribution to Ω from the radial field along a vertical betatron trajectory is dominated by the Fourier component at κ with resonance strength ǫ(J 2 ). The SRM can be solved exactly and the ISF is given by [30] is the difference between the vertical betatron phase and the phase of the Fourier component and (ê 1 ,ê 2 ,ê 3 ) are horizontal, vertical and longitudinal unit vectors. The tilt ofn away from the verticaln 0 is | arcsin(ǫ/λ)| so that it is 90 • at δ = 0 for non-zero ǫ. At large |δ|, the equilibrium polarisation directionsn(J 2 , φ 2 ; s), are almost parallel ton 0 (s) but as we see from the above formula, at δ = 0,n lies in the horizontal plane and P lim = 0. In this simple model ν s exists and is well defined near spin-orbit resonances for all Q 2 . In our calculations we choose the phase of the Fourier harmonic to be zero so that φ 2 represents the phase of the vertical betatron motion.
It is found both in practice and in simulation, that in the absence of special measures, acceleration of the beam through δ = 0 at practical rates can lead to loss of beam polarisation. This loss can be ascribed to a loss of invariance of J s and it can be quantified in terms of the Froissart-Stora formula [25]. Luckily, the loss of polarisation can be reduced by installing pairs of Siberian Snakes [11,12], magnet systems which rotate spins by π, independently of z, around a "snake axis" in the machine plane. For example, one puts two snakes at diametrically opposite points on the ring. Thenn 0 ·ê 2 = +1 in one half ring and −1 in the other. With the snake axes relatively at 90 • , the fractional part of ν 0 becomes 1/2 for all γ 0 . For calculations one often represents the snakes as elements of zero length ("point-like snakes"). Then if, in addition, the effect of vertical betatron motion is described by the SRM, and orbital resonances are avoided, at most J 2 , the fractional part of the ADST is 1/2 too, independently of γ 0 [10,31,5]. This is a special feature of this model. Thus for [Q 2 ] away from 1/2, the system is not at the first order spin-orbit resonance ν s (J 2 ) = [Q 2 ]. Therefore such resonances are not crossed during acceleration through δ = 0 and the polarisation can be preserved. This is confirmed by tracking simulations. However, simulations have shown also that the polarisation can still be lost if [Q 2 ] =ã 2 /2b 2 where here, and later,ã 2 andb 2 are odd positive integers withã 2 < 2b 2 [13,14,15]. This is the "snake resonance phenomenon" and it has also had practical consequences [13,14,15,17,18], especially for smallb 2 . Such a [Q 2 ] fits the condition 1/2 = (1 −ã 2 )/2 +b 2 [Q 2 ]. Since such tunes correspond to orbital resonance an ADST does not exist at most amplitudes. Then, according to our definition the system is not on a spin-orbit resonance ν s (J 2 ) = (1 −ã 2 )/2 +b 2 [Q 2 ]. However, for nearby irrational [Q 2 ] an ADST can exist, namely with the value 1/2. Then one can say that the system is close to spin-orbit resonance. This case is studied in Part II. Because the system is on orbital resonance and using the analogy with the smooth systems [7], even a smoothn need not be globally unique. Even if it were, there would be no guarantee that the maximum time averaged polarisation on a torus would be given by | n( z; s) |. We investigate these matters in the next section. Note that the rings in the Relativistic Heavy Ion Collider, RHIC [17,18] contain two snakes and that the RHIC team has avoided running near snake-resonance vertical tunes. Even away from the dangerous orbital tunes just mentioned, snake layouts should be chosen carefully. Methods for choosing layouts are discussed in [3,4].
Although one can describe spin motion in terms of orthogonal 3x3 matrices, here, we prefer to use SU(2) matrices. Correspondingly, the orientation of a spin is encoded in a two-component spinor 2 . We write the SU(2) matrices as where I is the 2 × 2 unit matrix,m is the unit vector along the effective rotation axis, ψ is the angle of rotation around that axis and the three components of σ are the Pauli matrices. The rotation is right handed when ψ > 0. Equation (1) can be re-written as where 3 i=0 r 2 i = 1. We call the real ordered quadruple (r 0 , r) a unit quaternion [32,4]. Spin maps are then concatenated using the multiplication rule where (a 0 , a), (b 0 , b) and (c 0 , c) are unit quaternions. The elements of the usual 3 × 3 matrices are given by R ij = (2r 2 0 −1)δ ij +2r i r j +2r 0 ǫ ijk r k where δ ij is the Kronecker symbol and ǫ ijk is the Levi-Civita symbol. Note that the R ij are homogeneous quadratic forms in the r i . This implies that R ij (r 0 , r) = R ij (−r 0 , − r) which simply reflects the fact that SU(2) covers SO(3) twice. In this paper, as in [5], we consider a system with two point-like snakes placed at diametrically opposite points on the ring. The snake axes are respectively at 0 • and 90 • to the longitudinal direction. The effect of vertical betatron motion is modelled by the SRM. The components of the unit quaternion for one turn starting with phase φ 0 2 just before the first (0 • ) snake are then As mentioned above, on orbital resonance, the vectorn can be obtained (up to a sign) as the eigenvector of unit length of the appropriate c-turn spin map. In terms of unit quaternions,n is simply the unit vector along the vector r (c) for the c-turn unit quaternion and we are free to choose the sign. It is clear from (4) that with small but non-zero ǫ/λ, the 1-turn spin map is close to a rotation by the angle π around an axis close to the vertical. This is expected on physical grounds too: at large |δ|, i.e., far from the parent resonance, or at small ǫ, the perturbation embodied in ǫ is relatively unimportant and the spins precess by an amount per turn similar to that on the design orbit. Then, the map for an odd number of turns is also close to a rotation by the angle π around the vertical but the map for an even number of turns is close to the identity. If λ is an even integer, the 1-turn spin map is always a rotation by the angle π around the vertical.

Polarisation in the model ring at rational [Q 2 ]
We now use our model to study and contrast the equilibrium beam polarisation, the time averaged beam polarisation and the beam polarisation surviving after acceleration, for the first members of the three classes of rational tunes just listed, namely for [Q 2 ] = 1/3, 1/4 and 1/6. We are primarily interested in [Q 2 ] at and near 1/6 but the other cases serve to familiarise the reader with the "normal" cases.

Off orbital resonance
To set the scene, and at variance with the title of this section, we first consider a case where the system is off orbital resonance and off spin-orbit resonance so that the smooth ISFn is globally unique. Thus figure 1 shows the components ofn for δ = 0 in the range 0 < [φ 2 /2π] ≤ 1 obtained by stroboscopic averaging [1,2,3,4] at the irrational tune Q 2 = 47+ √ 5−2 = 47.236067977 . . . 3 . In this and in all other figures in this paper, the spins are viewed just before the 0 • snake. Furthermore, for all calculations in this paper, the resonance strength, ǫ, is 0.4 and the integer k 0 is 1800, corresponding to a proton energy of about 970 GeV. These are the values used in [5] and we use them again here to allow comparisons to be made.
We remind the reader thatn is 2π-periodic in φ 2 . In principle, the stroboscopic averaging could have been carried out at each value of [φ 2 /2π] separately. However, away from orbital resonances one can cover a torus by simply findingn at some [φ 2 /2π], setting a spin parallel to thisn and then recording the spin components while transporting the spin for a large number of turns. Since J s is invariant along a trajectory we then have the components ofn all along the trajectory. This is the approach adopted for figure 1 and we see confirmation thatn is a single valued continuous function of [φ 2 /2π]. The average n ofn over φ 2 is vertical and P lim = 0.47. The ADST is 1/2. Figure 2 shows the beam polarisation, sampled every hundred turns for 10 6 turns, for an ensemble of particles distributed uniformly in the range 0 < [φ 2 /2π] ≤ 1 at δ = 0 when the  spins are all initially vertically upward. The horizontal components remain at zero but the vertical component oscillates, at least for millions of turns, between time independent maxima and minima with a time average of about 0.3. As expected, this is less than P lim . A constant polarisation equal to the maximum 0.47 could have been attained by setting the spins initially parallel to their respectiven vectors. See also figure 9 in [1]. Inspection of the turn-by-turn data reveals that the oscillations have a period of about four turns, as expected for a [Q 2 ] close to one quarter and an ADST of 1/2. In the simple SRM and at δ = 0 the analogous simulation would exhibit a beam polarisation oscillating between +1 and −1 as the spins precessed around the horizontaln at a rate λ = ǫ.  Figure 3 shows the components ofn for the parameters of figure 1 except with δ = 10.6, a value corresponding to a beam energy far from that of the parent resonance, with non-even λ but otherwise arbitrary. The vectorsn(φ 2 ) are almost vertical so that P lim is high, namely 0.998. Figure 4 shows the curve for P lim together with the beam polarisations, as ensembles are accelerated through δ = 0 at the rates of 100 MeV, 500 MeV and 1 GeV per turn (p.t.). The acceleration is simulated by incrementing δ by four equal amounts, namely just after each snake and at the mid-points of the two arcs. At the start, δ = −10.6 and the particles are distributed uniformly in [φ 2 /2π] with each spin initially set parallel to its correspondingn(φ 2 ), which is almost vertical. For protons, a rate of 100 MeV per turn corresponds to ∆ ≈ 0.19 for the change of aγ 0 per turn. For this rate the beam polarisation follows the curve for P lim vs. δ, dipping to the value 0.47 at δ = 0. Moreover, detailed inspection shows that at each δ the distribution of spins matches the ISF. This is a nice demonstration of the adiabatic invariance of J s in this case [9]. The invariance of J s is lost at the higher rates. Slightly different curves are obtained if the spins are set vertically upward at the start.
The rate of 100 MeV per turn corresponds to a value ǫ 2 /α ≈ 5.3 in the Froissart-Stora formula [25] where α = ∆/2π. The Froissart-Stora formula describes the final polarisation when a spin-orbit resonance is crossed in the SRM and for these parameters it would predict almost full spin flip, corresponding to adiabaticity. However, our model includes the snakes and there are therefore no first order spin-orbit resonances to cross. So the Froissart-Stora formula does not apply. Nevertheless for our model, the rate of 100 MeV per turn is adiabatic.

On orbital resonance: [Q 2 ] = 1/3
We now consider our first case of orbital resonance, namely with Q 2 = 47 + 1/3, corresponding to odd a 2 and b 2 . Figure 5 shows the components ofn at δ = 0 and ǫ = 0.4. These components are obtained by normalising to unity the r (3) corresponding to three turns in the range 0 < [φ 2 /2π] ≤ 1/3, namely 0 • to 120 • , and then transporting then for each [φ 2 /2π] in this range for two or more turns with the 1-turn spin map, thereby filling up the full phase range. Note that the curves are single valued functions of [φ 2 /2π] as required. The average | n | ofn over [φ 2 /2π] in figure 5 is 0.05 and n is vertical. While the smooth ISFn of figure 5 is globally unique, one looses global uniqueness if one allows discontinuities, as demonstrated in figure 6. There, we introduce changes of sign inn by hand at the arbitrarily chosen angles of 17.5 • and 90 • , while constructingn in the range 0 • to 120 • using r (3) . We then transport thisn for two or more turns as before. Naturally, the sign-discontinuities (often simply called "discontinuities" from now on) are transported too. In particular, we see that the transportedn is still a single valued function of [φ 2 /2π]. The average | n | in figure 6 is 0.164. It is clear that neithern nor | n | are unique. Of course, an unlimited number of discontinuities could be introduced in the same way. Then the curves would be smooth almost nowhere. Each of then obtained in this way would correspond to a permissible equilibrium spin distribution. Then obtained by stroboscopic averaging [1] over the whole range 0 < [φ 2 /2π] ≤ 1 can have discontinuities with positions that depend on the "seed" spin field used in the stroboscopic average but these discontinuities can be removed to give the curves in figure 5. Since these discontinuities are sign discontinuities, we do not exclude the possibility that the ISF is locally unique. However, this issue is beyond the scope of this paper since it would lead us into a discussion of regularity conditions. If the long term tracking simulation of figure 2 is repeated but with [Q 2 ] = 1/3, the vertical component of the beam polarisation oscillates quickly between about -0.3 and +0.8 for at least 5 · 10 6 turns with a time average of about 0.25. This is higher than the | n | in figure 5 but no significance can be attributed to this since | n | is not unique. Figure 7 shows the components of the smooth ISFn for the conditions of figure 5 but with δ = 10.6. | n | is high as expected, namely 0.997 since r (3) is close to vertical. The existence of then of figure 7, means that an ensemble of exactly vertical spins is close to a permissible equilibrium spin distribution. Figure 8 shows the beam polarisation for acceleration through δ = 0 from δ = −10.6 to δ = +10.6 at the rates of 50 MeV, 300 MeV and 1 GeV per turn for this Q 2 . At the start, the particles are distributed uniformly in [φ 2 /2π] and the spins are set parallel to the almost vertical n vectors of the smooth ISF. Up to an acceleration rate of 50 MeV per turn, J s is invariant, with the beam polarisation dipping down to 0.05 around δ = 0 and returning to a high value at the end. This is a demonstration that with the chosen smoothn, J s can be adiabatically invariant, although the proof in [9] does not guarantee this because the system is on orbital resonance. At the higher acceleration rates, the invariance is lost. By using stroboscopic averaging for irrational [Q 2 ] near 1/3 one finds ISFs similar to that in figure 5.

On orbital resonance: [Q 2 ] = 1/4
For our second case of orbital resonance we choose Q 2 = 47 + 1/4, corresponding to an odd a 2 and a b 2 which is twice an even integer. Figure 9 shows the components ofn at δ = 0 and ǫ = 0.4 obtained, in analogy with the previous case, from r (4) in the range 0 < [φ 2 /2π] ≤ 1/4 and from transporting thosen for three or more turns. In this case we see "stray" points at multiples of 45 • corresponding to the phases where the 4-turn map is the identity. For this figure we have imposed the constraint that the components are continuous in the range 0 • to 90 • , apart from the stray points. If we had not imposed smoothness, the components would have changed sign at 45 • and the resulting discontinuities would have been transported to the remainder of the phase range. So, for these parameters and for [Q 2 ] = 1/4,n can have discontinuities as in the case of any rational Q 2 . But in contrast to a case discussed below, these discontinuities can be suppressed. Then obtained by stroboscopic averaging over the whole range 0 < [φ 2 /2π] ≤ 1 is smooth as in figure 9. Of course, as in the case of [Q 2 ] = 1/3, we can also introduce an unlimited number of sign-discontinuities. The curves of figure 9 give | n | = 0.43. Note the . This is consistent with the prediction in [4, Section 4.8] that in mid-plane symmetric rings the ISF is well behaved close to the condition ν 0 = k 0 + 2k 2 Q 2 , (k 0 , k 2 ∈ Z). If the long term tracking simulation of figure 2 is repeated but with [Q 2 ] = 1/4, the vertical component of the beam polarisation oscillates quickly, initially between about -0.1 and +0.7. But these limits gradually change and become 0.1 and 0.4 respectively after 5 · 10 6 turns. The time average of about 0.25. This is lower than the | n | in figure 9 but no significance can be attributed to this since | n | is not unique. Figure 10 shows the components ofn for the conditions of figure 9 but with δ = 10.6. The average | n | is 0.99. Note that in contrast to the 3-turn map used for [Q 2 ] = 1/3, at large |δ| the 4-turn map is close to the identity. Nevertheless, r (4) is close to vertical.  At the start, the spins are set parallel to the almost verticaln vectors of the smooth ISF. Up to an acceleration rate of 100 MeV per turn, J s is invariant, with the beam polarisation dipping down to 0.43 around δ = 0 and returning to a high value at the end. This is again a demonstration that with the chosenn, J s can be adiabatically invariant although the system is on orbital resonance. At the higher acceleration rates, the invariance is lost.  We now come to the first of the two cases of primary interest for this study, namely the case when [Q 2 ] = 1/6, i.e., a case of a snake resonance. Again, the integer part of Q 2 is 47 and ǫ = 0.4. Figure 12 shows the components ofn at δ = 0 obtained by transporting for five or more turns then obtained from r (6) in the range 0 < [φ 2 /2π] ≤ 1/6. We see stray points at phases which are multiples of 30 • and 90 • corresponding to the phases where the 6-turn map is the identity. The vector r (6) has sign-discontinuities at these points but for this figure we have imposed the constraint that the components ofn are continuous in the range 0 • to 60 • , apart from the stray points. One sees thatn still has discontinuities, namely at phases which are multiples of 60 • . Thus, in spite of smoothingn in the initial range of 0 • to 60 • , discontinuities persist. They cannot be removed without creating a vector field which becomes double valued when it is transported turn-by-turn. However, the discontinuities can be moved. These effects explain the failure of the MILES algorithm forn at snake-resonance tunes in [21] where the need for discontinuities in this model is nevertheless disputed. It is clear that the curves in figs. 7 and 8 in [21] do not representn [5].
Stroboscopic averaging over the whole range 0 < [φ 2 /2π] ≤ 1 generates the curves of figure  12 directly i.e., without extra smoothing. The discontinuities ofn occur at phases where the raw stroboscopic average passes through zero. The passage through zero is smooth. So discontinuities inn do not imply discontinuities in the stroboscopic average.
Our numerical calculations show thatn has such discontinuities at snake-resonance tunes at most values of ǫ and that the minimum number of discontinuities is 2b 2 .
Of course, ifn is represented as the locus of points on the unit 2-sphere, one finds disjoint segments. The average | n | over [φ 2 /2π] in figure 12 is 0.13. An arbitrary number of extra discontinuities can be introduced by hand.
If the long term tracking simulation of figure 2 is repeated but with [Q 2 ] = 1/6, the polarisation oscillates quickly, but with constant upper and lower limits with a time average of about 0.1, at least up to 5 · 10 6 turns. Thus the time averaged polarisation does not vanish. This is illustrated in figure 13.  Figure 14 showsn obtained as for figure 12 but with δ = 10.6. Except when λ is an even integer this is typical of then at large |δ| (and also at small ǫ). The value of r (6) is very small and the 6-turn spin map is close to a rotation of 2π aroundn. The discontinuities persist but in contrast to the earlier examples, the vertical component ofn is close to zero and the horizontal components are piece-wise almost independent of [φ 2 /2π]. The average | n(φ 2 ) | is essentially zero. It would remain close to zero if sign-discontinuities were introduced by hand. Since the horizontal components ofn are piece-wise almost independent of [φ 2 /2π] but also different, and since the 1-turn spin map is a rotation of about π around an axis close to the vertical, it essentially changes their signs from turn to turn, causing the discontinuities. Such  figure 15 where we repeat the long term tracking simulation of figs. 2 and 13 but at δ = −10.6 and [Q 2 ] = 1/6. We now see that the polarisation falls, but slowly, over many tens of thousands of turns and subsequently oscillates around zero. Then the time averaged polarisation is close to | n(φ 2 ) | ≈ 0. Nevertheless, since the system is on orbital resonance, the theorem [3,4] on the maximum time averaged polarisation does not enforce this. Although the initial spin distribution is not        in equilibrium, it is not surprising that it takes about 10 5 turns before the polarisation reaches zero. This is due to the fact that at large |δ| the eigentune, 6ν 6 , of the 6-turn spin map is almost independent of [φ 2 /2π] and very close to an integer for this case. Since J s is invariant along a trajectory, we can view the motion of a spin as a precession at a fixed angle cos −1 ( S ·n/| S|) around itsn. In this case the angles are about 90 • . With eigentunes almost independent of [φ 2 /2π] and close to an integer, the projections of spins on the planes perpendicular to their respectiven's spread out (decoher) only slowly. Then, at the viewing position, the spins return almost to their original directions after six turns. For large |δ|, the 1-turn spin map corresponds to a rotation of about π around an axis close to the vertical. So, it is again no surprise that the polarisation in figure 15 takes many turns to reach zero. For even larger |δ| (e.g., over 100),n can be taken to be horizontal but the polarisation remains vertical and it takes many millions of turns for it to show signs of falling. There is no fall if λ is an even integer since then, the 6-turn map is the identity. Figure 16 shows the beam polarisation for acceleration through δ = 0 at the rates of 50 KeV, 10 MeV and 500 MeV per turn for [Q 2 ] = 1/6. At the start, the particles are uniformly distributed in [φ/2π] and the spins are set parallel to the almost horizontaln vectors of that ISF which deforms into the ISF's of figures 12 and 14. The initial beam polarisation is essentially zero. During acceleration at rates up to 50 KeV per turn, the beam polarisation rises to 0.13, corresponding to the | n | of figure 12, and then returns to around zero. A detailed inspection of the data shows that for a rate of 10 MeV per turn, the spins deviate slightly from their respectiven vectors at large |δ|. However, this effect is not apparent in the average over [φ 2 /2π] contained in the beam polarisation. This is again a demonstration that with the chosen n and the chosen layout of accelerating cavities, J s can be approximately invariant even for these discontinuous ISF's and that at the higher acceleration rates, the invariance is lost. The approximate invariance is confirmed in figure 17 which shows the corresponding behaviour of the phase average of J s , J s . In figure 17 we have suppressed data at δ's wheren is indeterminate because λ is an even integer. Figure 18 shows the beam polarisation as the simulation of figure 16 is repeated but with the spins initially vertically upward and for rates of 50 KeV and 10 MeV per turn and for 50 MeV per turn, where J s is still approximately invariant. For these rates of acceleration the angle between a spin and itsn remains around 90 • . Then the beam polarisation during acceleration depends just on the geometry of the ISF and on the history of the rate of decoherence of the projections of the spins on the planes perpendicular to then's. These rates depend, in turn, on the magnitude of 6ν 6 and its dependence on [φ 2 /2π]. We therefore expect that the final polarisation could depend sensitively on the magnitude of the rate of acceleration and on its time dependence. This is confirmed in figure 18 which shows that at a rate of 50 KeV per turn, the polarisation is effectively lost at positive δ but that at the much higher rate of 50 MeV per turn the final polarisation is around -0.4 at the end of the acceleration cycle. By now, the reader will have realised that the polarisation of -0.4 cannot represent an equilibrium state. This is confirmed in figure 19 where, after acceleration up to δ = 10.6, δ is frozen and the ensembles are tracked for a further 5 · 10 6 turns. Figure 19 shows that after some large oscillations the polarisation gradually decays to zero in a way and on a time scale familiar from figure 15. It also shows that although the polarisation can be small at the end of the acceleration (as in the case of 10 MeV/turn), the spin distribution is by no means isotropic but is such that the polarisation can return to a large value later. In fact after the 5 · 10 6 turns, the curves of spin vector versus [φ 2 /2π] are smooth curves for all three acceleration rates 4 . This suggests that contrary to conventional expectation, a complete loss of polarisation is not inevitable during acceleration exactly at a snake resonance with [Q 2 ] = 1/6, at least not within the confines of our model. This completes Part I of our investigation.

Summary and conclusion
In this paper we have presented and contrasted four scenarios for spin motion on and off orbital resonance within the confines of our simple model, and by this means we have developed a clean, elegant account of the special features of spin motion at a snake resonance. In all four cases S ·n is an invariant at low enough rates of acceleration. For the first three cases ([Q 2 ] = 0.236067977 . . . , 1/3, 1/4) the ISF is close to vertical at large |δ|, i.e., far away from the energy for the parent resonance, and the spin motion is unexceptional. For example, after acceleration from a large negative δ to a high positive δ, an initially vertical spin is still close to vertical. These cases serve to emphasise the exceptional form of the ISF when [Q 2 ] = 1/6. In this case, far away from the parent resonance, the ISF lies close to the horizontal plane. Then in contrast to the other three cases, an ensemble of particles with a uniform distribution of [φ 2 /2π] and with vertically upward spins, cannot be at spin equilibrium. The subsequent evolution of the beam polarisation depends on the chosen initial δ and is exemplified in figures 13 and 15. In particular, the polarisation oscillates at a rate depending on the proximity of the eigentune of the 6-turn spin map to an integer and on the extent of the variation of that eigentune with [φ 2 /2π]. Then at the energy of the parent resonance (δ = 0), the polarisation oscillates quickly and the time averaged polarisation is small but non-zero. At most large |δ|, the time averaged polarisation is zero but the polarisation oscillates slowly and it reaches zero for the first time only after many thousands of turns.
As soon as one sees that at most large |δ| the ISF for [Q 2 ] = 1/6 lies close to the horizontal plane, it is no surprise that in this case the time averaged beam polarisation can become small in the long term. Acceleration adds little to the story, except that within our model, after starting with an ensemble of vertical spins at δ = −10.6, the final polarisation depends on the rate at which one passes from the spin motion underlying figure 15 to the spin motion underlying figure 13 and then beyond to large positive δ. The key features of spin motion at [Q 2 ] = 1/6 are encoded in the ISF. We see no necessity to invoke the perturbed spin tune [14,15]. Instead, we appeal to the eigentune of the 6-turn spin map, a quantity with physical significance.
We emphasise that the main results presented here refer to a very special case, namely for our model right at [Q 2 ] = 1/6 and with ǫ = 0.4. As pointed out in [5], the ISF is extremely complicated for values of [Q 2 ] just below and just above 1/6. This is consistent with the prediction in [4, Section 4.8] that in mid-plane symmetric rings the ISF need not be well behaved close to the condition ν 0 = k 0 + (2k 2 + 1)Q 2 , (k 0 , k 2 ∈ Z). Thus in Part II of this study we extend our calculations to cover such values of [Q 2 ] and to larger values of ǫ. It will be shown there that although the ISF for [Q 2 ] = 1/6 has the special form described above, this is an exception and that the loss of polarisation during acceleration near to [Q 2 ] = 1/6 has a different origin. We also comment on the findings in [10,23,24].
The analysis should then be extended to real synchrobetatron motion with misalignments for a typical optic of a real ring and with the fields of real snakes. See, for example, [33].
Other snake-resonance tunes should also be covered. We note with interest that according to simulations for RHIC, the loss of polarisation during acceleration is less severe when the simulations are carried out with the magnetic fields of real snakes rather than with point-like snakes [34].