Accelerators vs. ergodicity

(Sethna, “Entropy, Order Parameters, and Complexity”, ex. XXX)
© 2021, James P. Sethna, all rights reserved.

Synchrotrons push charged particles around circular orbits at near the speed of light, stably over billions of orbits. Why haven’t the nonlinear fields used to steer the particles kicked them out of the beam? Like the cat map in Exercise 5.8 and the three-body problem of Exercise 4.4, physicists distill the dynamics using a Poincaré section into a two-dimensional area-preserving nonlinear map. The classic example, due to Chirikov, is the so-called “standard map”
   (Subscript[θ, n+1])=  (Subscript[θ, n])+ (Subscript[p, n+1]) =  (Subscript[θ, n])+ (Subscript[p, n])+K sin((Subscript[θ, n])) 
   (Subscript[p, n+1]) =  (Subscript[p, n])+ K sin((Subscript[θ, n]))
Here  p is is proportional to the squared amplitude of the deviation from center of beam and  theta represents the phase of the oscillation around that center. Large  K  represents a strong nonlinearity in the fields. For plotting, we shall take θ mod 2 π,  and occasionally also p mod 2 π.

We remember from Exercise 5.8 that area-preserving maps mimic the volume-preservation in phase space guaranteed by Liouville’s theorem.

(a) Show that the standard map preserves area, by showing that its Jacobian has determinant one.

vars = {θ,p};
ChirikovNoMod[K_][{θ_,p_}] := {...,...}
jacobian = Table[D[...,vj],{i,1,2},{vj,vars}]
Det[jacobian]

(b) Iterate the standard map a few thousand times for K=4 , starting with ((Subscript[θ, 0]), (Subscript[p, 0]))=(0.1,0.11)  (chosen arbitrarily), and plot (Subscript[p, n]) versus n for a thousand points n. Does the amplitude of the particle oscillations around the beam path remain bounded in amplitude as  n  gets large? Would designing an accelerator with this large an anharmonicity be a good idea?

pts = NestList[ChirikovNoMod[4], {...,...}, 1000];
p = Transpose[pts][[2]];
ListPlot[p, AxesLabel->{"iterations","p"}, PlotLabel->"K=4"]

The situation changes for smaller anharmonicity.

(c) Implement the standard map graphically, so that you can interactively select a variety of initial conditions and run each for a few thousand steps to explore the trajectories, using  mod 2 π,  for θ but not for  p . Check that you get trajectories for K=0.2 that are compatible with those shown in the figure in the text. Plot the behavior for K=0.7 , selecting a variety of initial conditions. Be sure to illustrate “oval” KAM tori around stable periodic orbits, “horizontal” KAM tori that span from θ=0 to 2π, and chaotic regions. (Remember that your trajectory is a cross section once per orbit around the accelerator. So a closed curve becomes a tube, which mathematically is called a torus. KAM stands for Kolmogorov, Arnol’d, and Moser, who proved that tori survive for certain irrational winding numbers, see part e.) 
Explore initial conditions spanning p in (0, 2 π) . Do the trajectories stay bounded in  p? Which keeps  p  from growing indefinitely -- the oval or horizontal tori? 
Keeping the magnets close to the beam is best -- until the beam starts hitting the magnets.  p  can be used to estimate the maximum deviation from the beam center. What range  Δp  would you use to design an accelerator at K=0.7?

Chirikov[K_][{θ_,p_}] := {Mod[...,2π],...}

ExploreStandardMap[K_,nPoints_:1000,initialPoint_:{0.1,0.11}] :=
DynamicModule[{pts={NestList[Chirikov[K],initialPoint,nPoints]}},ClickPane[Dynamic@ListPlot[pts,PlotRange->{{0,2π},All}, AxesLabel->{"θ","p"},PlotLabel ->Row[{"K=",K}]],(pts=Append[pts,NestList[Chirikov[K],#,nPoints]])&]]

ExploreStandardMap[0.7,10000]

The equilibrium state predicted from statistical mechanics would fill all of phase space uniformly. The map at  K=0.7  is not ergodic.

In Section 4.2 we defined that a time evolution was ergodic if and only if all the ergodic components of the space either have zero volume or have a volume equal to the space. Here ergodic components were sets R that remain invariant under the map (so if  ((Subscript[θ, 0]), (Subscript[p, 0])) is in R then ((Subscript[θ, n]), (Subscript[p, n])) is also in R).

(d) Identify a few ergodic components in your plot for K=0.7 in part (c). Are the KAM tori ergodic components? Do they have non-zero volume? Do they surround ergodic components with non-zero volume? Does the chaotic regions surrounding θ ~ p ~ 0  appear to be an ergodic component with non-zero volume? 

(e) Alter your standard map to also take  p mod 2 π , to keep the motion bounded. Plot and examine carefully the time evolution at K=4 . Do all initial conditions fill out the entire volume? If not, describe ergodic components that neither have zero volume nor fill the square.

ChirikovModp[K_][{θ_,p_}] := {Mod[...,2π],Mod[...,2π]}
ExploreStandardMapModp[K_,nPoints_:1000,initialPoint_:{0.1,0.11}] :=
DynamicModule[{pts={NestList[...[K],initialPoint,nPoints]}},ClickPane[Dynamic@ListPlot[pts,PlotRange->{{0,2π},All}, AxesLabel->{"θ","p"},PlotLabel ->Row[{"K=",K}]],(pts=Append[pts,NestList[...[K],#,nPoints]])&]]

ExploreStandardMapModp[4,10000]

We define the winding number of a particle trajectory as the average number of oscillations θ / 2π   per iteration of the map. The KAM theorem tells us that the tori whose winding numbers are difficult to approximate by rationals are the most stable.

(f) What is the winding number at  K=0,  for a trajectory starting at initial condition  ((Subscript[θ, 0]),(Subscript[p, 0]))? Write a routine to approximate the winding number for  K=0.7 and  θ = π, and plot it for  0< p <2π . By examining your plot of the different trajectories for K=0.7  from part (c), explain the plateaus you observe at rational winding numbers.

WindingNumber[K_,p0_,nPoints_, θ_:π] := Nest[ChirikovNoMod[K],{θ,p0},nPoints][[1]]/ ...
Plot[WindingNumber[0.5,p0,100],{p0,0,2π}, PlotLabel->"Winding number, K=0.7", AxesLabel->{"p","winding number"}]

The last ``horizontal’’ KAM torus for the standard map is destroyed by the nonlinearity at  K=0.971635… . Its winding number is the inverse Golden Mean,  ((Sqrt[5])-1) / 2 = 0.618033… , which is the most irrational number -- the number hardest to approximate by ratios of small integers. This transition is associated with self-similar behavior and scaling, and has been studied using the renormalization group methods we study in Chapter 12.

The following is graphics for the figure in the exercise.

R = 1;
A = 0.05;
k =2π+ (Sqrt[5]-1)/(4π);
center = ParametricPlot[{R Cos[ϕ], R Sin[ϕ]},{ϕ,0,2π}, PlotStyle->{Dashing[Large], Thickness[0.01]}, Axes->False];
orbit = ParametricPlot[{(R+A Cos[k ϕ])Cos[ϕ], (R+A Cos[k ϕ]) Sin[ϕ]},{ϕ,-0.1,2π+0.1}, PlotStyle->{Black, Thickness[0.01]}, Axes->False];
dots = ListPlot[Table[{R Cos[ϕ], R Sin[ϕ]}, {ϕ, π/(2 k), 2π, 2π/k}], Axes->False, AspectRatio->1, PlotStyle->{Green,PointSize[0.03]}];
section = Plot[0,{r,R-3A,R+3A}, PlotStyle->{Red, Thickness->0.02}];
AcceleratorOrbit = Show[center,orbit, dots,section]

SetDirectory["/Users/sethna/6562/Exams/Prelim/Fig"]
Export["AcceleratorOrbit.pdf",AcceleratorOrbit]
Export["AcceleratorOrbit.png",AcceleratorOrbit, ImageResolution->500]