Homework 8 Discussion

Some General Comments:

• grade: 20 possible
  + 20 pts each problem (10 each, if “re-done”)
  + ave score: 18/20

• note that if a gradient field “defocuses”, then it decreases the tune; if it focuses, then in increases the tune

For our Final Exercise Set (prior to our “Final Project”), select ONE of the three options below. Only one option will be accepted for grading:

1. Option 1: Dynamic Aperture Calculation (below)

2. Option 2: Beam-Beam Interaction Calculation (below)

3. Option 3: Re-write and re-submit TWO of your previous homework problems (not complete assignments; rather, two problems from previous assignments). They must be problems for which you received less than half-credit (e.g., 3/10 points, for example). Simply copying the answers from the web site is not acceptable; a good-faith attempt must be made to write and explain the solutions in your own words, with a demonstration of your mastery of the methods you’ve learned during the course.

As usual, the assignment must be submitted as an html R notebook (yourfirstname8.nb.html) which will include all code (R and/or python) produced by the author.

Dynamic Aperture Calculation

L  = 15
mu = 80.3*pi/180
F  = L/2/sin(mu/2)
Ncells = 26
nu = mu*Ncells/2/pi
C = 2*L*Ncells

Consider a FODO-style synchrotron, made up of 26 FODO cells, with a lens spacing of 15 m and with a cell phase advance of 80.3 degrees. Presume equal bending through angle \(\theta\) between each of the quadrupoles of the FODO system.

1. Give formulas for and compute the maximum and minimum values for \(\beta\) and \(D\) (dispersion) for such a FODO cell.

Mdf =  matrix( c( 1, L, -1/F, 1-L/F), byrow = TRUE, ncol = 2)
Mfd =  matrix( c( 1, L,  1/F, 1+L/F), byrow = TRUE, ncol = 2)

ang = 2*pi/Ncells/2
Dmax = L*ang/sin(mu/2)^2*(1+1/2*sin(mu/2))
Dmin = L*ang/sin(mu/2)^2*(1-1/2*sin(mu/2))
betMax = 2*F*sqrt((1+L/2/F)/(1-L/2/F))
betMin = 2*F*sqrt((1-L/2/F)/(1+L/2/F))

2. To the extent that the natural chromaticity (in the absence of errors) of a FODO-style synchrotron is \(\xi \approx -\nu\), what would be the values \(\xi_x\) and \(\xi_y\) for our system?

xix = -Ncells*mu/2/pi
xiy = -Ncells*mu/2/pi  # from symmetry, they will be the same
xix; xiy

[1] -5.799444
[1] -5.799444

3. Suppose sextupoles are placed at the locations of the quadrupoles for chromaticity control, with strengths \(S_1\) at the focusing quadrupoles and \(S_2\) at the defocusing quadruoples. Write a matrix equation that shows how the chromaticity would be adjusted in our system for these two given strengths.

4. If the change desired is to set both of the chromaticity values to zero, what numerical strengths \(S_1\) and \(S_2\) are required?

dxix = mu/2/pi  #  need to make "positive" changes to each cell
dxiy = mu/2/pi  #   this correction in each cell would give full 
                #   correction for Ncells

A = matrix(c( betMax*Dmax,  betMin*Dmin,
             -betMin*Dmax, -betMax*Dmin), 
    byrow=TRUE, ncol=2)/4/pi
Sxt = solve(A,c(dxix,dxiy))
S1 = Sxt[1] # 1/m^2  here, S = B''L_sxt/(Brho)
S2 = Sxt[2] # 1/m^2
S1; S2

[1] 0.01238819
[1] -0.02417646

Ntrk = 5000
x0  = seq(-100,100, 2)
xp0 = 0
xmax = 100

5. With these values of sextupole strengths, we now have a nonlinear system, which can affect the available stable phase space area. For an initial condition of \(x = x_0\) and \(x'\) = 0, at some value of \(x_0\) the motion is likely to become unstable. Track particles through the FODO system, including the effects of all the sextupoles, for a variety of values for \(x_0\), up to 5000 turns, and determine at what value(s) of \(x_0\) the motion becomes unstable. (Let’s call a particle unstable if its position \(x\) exceeds 100 mm within 5000 turns.) The region of \(x\) which contains stable motion could be called the 5000-turn dynamic aperture.

Note: \(\Delta x' = -\frac12 B''\ell/(B\rho)\)…

Nlast = 0

i=0
while(i<length(x0)){
   i=i+1
   V = c(x0[i],xp0)

   j = 0
   while(j<Ntrk){
      j = j+1
      m = 0
      while(m<Ncells){  # once around the ring
         m=m+1
         V = Mfd %*% V
         V[2] = V[2] - 1/2*S2*V[1]^2/1e3  # convert to mr
         V = Mdf %*% V
         V[2] = V[2] - 1/2*S1*V[1]^2/1e3  # convert to mr
         if(abs(V[1]) > xmax) break
      }
      Nlast[i]=j
      if(abs(V[1]) > 100) break
   }
}
plot(x0,Nlast,typ="l")

min(x0[Nlast>1000]); max(x0[Nlast>1000])

[1] -42
[1] 40

Beam-Beam Interaction

N   = 1e11
Ec  = 1e12      # collision energy, eV
mc2 = 938e6     # eV
gam = Ec/mc2
eps = 3*pi*1e-6 # m
betStar = 0.3   # m
sigx = sqrt(betStar*eps/gam)*1000  # mm
nu   = 20.5809
Qe   = 1.6e-19       # Coul
epsilon0 = 8.85e-12  # F/m
r0  = Qe/(4*pi*epsilon0*mc2)  # note: one "e" in `mc2` cancels one "e" in `Qe`^2
xibb  = N*r0/4/eps

Suppose a collider is operating with \(N\) = \(10^{11}\) protons in each bunch at a collision energy of \(1000\) GeV, colliding with similar counter-rotating proton bunches. The rms, normalized transverse emittance of each bunch is \(\epsilon_N\) = \(3\) \(\pi\) mm-mrad, and the amplitude function at the collision point is \(\beta^*\) = 0.3 m. The betatron tune of the collider is \(\nu\) = 20.99091.

The beam-beam interaction generates a shift in the tune of the collider, by an amount \(\Delta\nu = \xi\), where \[ \xi \equiv \pm\frac{r_0N}{4\epsilon_N} \] and \(r_0\) is called the “classical radius” of the particle, \[ r_0 = \frac{e^2}{4\pi\epsilon_0mc^2}. \] Here, we are actually considering a “round” beam, but will only study motion in one direction — \(x\).

1. For the situation described above, should we use the “positive” sign or “negative” sign for our tune shift? Justify your answer.

             # -1 if defocuses, tune goes down;  +1 if focuses, tune goes up
signXi = -1  #    if pos x pos or neg x neg, then -1;  if pos x neg, then +1
             #        here, we have protons on protons:  -1
xibb = xibb*ifelse(signXi <0, -1, 1)

1. As a particle passes through the interaction point, what would be the effective focal length, in the linear approximation, of the focusing field created by the oncoming bunch, in terms of \(\xi\) and \(\beta^*\)? Provide a numerical estimate.

\[ \Delta\nu = \frac{1}{4\pi}\beta q ~~~~~ \rightarrow ~~~~~q_{eff} = \frac{1}{F_{eff}} = \frac{4\pi\xi}{\beta^*} \]

Feff = betStar/4/pi/xibb
Feff  # m

[1] -5.867843

2. This extra focusing field would change the value of the beta function at the collision point from its design value of \(\beta^*\) = 0.3 m. By how much would it change?

dqEff   = 4*pi*xibb/betStar
dbonbAmp = dqEff*betStar/2/sin(2*pi*nu)
dbonbStar = -dbonbAmp*cos(2*pi*nu)
dbonbStar

[1] 0.04588251

3. We have seen that the beam-beam force (and space charge force) for a round beam has a functional form of \[ F \sim {\cal F}(u) = \frac{1-e^{-u^2/2}}{u} \] in terms of \(u \equiv x/\sigma\). Using R or python, make a plot of this functional form and its deriviate, and determine the value of \(x/\sigma\) that produces the maximum force.

Fbb  = function(u){ (1-exp(-u^2/2))/u}
Fpbb = function(u){ ( (1+u^2)*exp(-u^2/2) - 1 )/u^2 }
uMax = uniroot(Fpbb,c(1,4))$root
uMax

[1] 1.58521

curve(Fbb,0,10)
abline(h=c(0),v=uMax,lty=c(1,2),col=c("black","red"))

curve(Fpbb,0,10)
abline(h=c(0),v=uMax,lty=c(1,2),col=c("black","red"))

4. To study the phase space of the beam-beam interaction, note that the change in angle of a particle passing through an oncoming bunch will be \[ \Delta x' = \frac{F\cdot\ell}{pv} = \frac{e^2N'\cdot\ell}{\pi\epsilon_0 \sigma \;\; pv}\cdot {\cal F}(u) \] Noting also that the interaction region will be half the length of one of the beam bunches (due to the fact that the bunches are traveling toward each other), then \(N'\ell = N/2\) and \[ \Delta x' = \frac{e^2 \cdot (\frac12 N)}{\pi\epsilon_0 \beta^2\gamma mc^2\sigma}\cdot{\cal F}(u) = 2\;\frac{e^2}{4\pi\epsilon_0 mc^2}\;\frac{N}{\beta^2\gamma\sigma}\;{\cal F} = \frac{2r_0N}{\gamma\sigma} {\cal F} ~~~(\beta = v/c\approx 1) \] where we let \(\beta = 1\), which is most common for colliding beams storage rings.

If we now use our usual phase space variables \(x\) and \(P \equiv \beta x' + \alpha x\), then when passing through the oncoming bunch a particle’s \(P\) variable will change by \(\Delta P = \beta^*\Delta x'\), where \(\beta^*\) is the value of the Courant-Snyder parameter at the collision point. But, let’s re-write everything in terms of the rms beam size: define \(u \equiv x/\sigma\) and \(v \equiv P/\sigma\). Then,

\[ \Delta v = \Delta (P/\sigma) = \beta^*\Delta x'/\sigma = \frac{2\beta^*r_0N}{\gamma\sigma^2} {\cal F} = \frac{2\pi r_0N}{1}\frac{\beta^*}{\gamma\pi\sigma^2} {\cal F} = 2\pi\frac{r_0N}{\epsilon_N} {\cal F}(u) \] where we recognize \(\epsilon_N\) as the rms normalized beam emittance. The coordinates of a particle thus will change according to \[ \Delta u = 0 \\ \Delta v= 8\pi\xi\;{\cal F}(u) = 8\pi\xi \cdot\frac{1-e^{-u^2/2}}{u} \] when passing through the oncoming bunch. For a single pass around the collider, the particle’s phase space would evolve according to

\[ \begin{pmatrix} u \\ v \end{pmatrix}_{n+1} = \begin{pmatrix} ~~\cos (2\pi\nu) & \sin (2\pi\nu) \\ -\sin (2\pi\nu) & \cos (2\pi\nu) \end{pmatrix} \begin{pmatrix} u \\ v+\Delta v(u) \end{pmatrix}_{n} \]

dv = function(x){8*pi*xibb*(1-exp(-x^2/2))/x}

Now: Plot the \((u,v)\) phase space in a single figure for values of \(u_0\) = 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5 and for \(v_0\) = 0, where each condition is tracked for 2000 turns.

cn = cos(2*pi*nu)
sn = sin(2*pi*nu)
plot(0,0,typ="n",xlim=c(-alim,alim),ylim=c(-alim,alim),asp=1,xlab="u",ylab="v")
symbols(0,0,circles=sqrt(6),add=TRUE,fg="red",inches = FALSE)

m=0
while(m<length(u0)){
   m=m+1
   u = u0[m]
   v = 0
   i=0
   while(i<Ntrk0){
      points(u,v,pch=".")
      i=i+1
      u1 =  u
      v1 =  v + dv(u)
      u  =  cn*u1 + sn*v1
      v  = -sn*u1 + cn*v1
   }
}

ktune = c(20.50291, 20.580902, 20.6672073, 20.7498, 20.800776, 20.99091)
u0 = c(2:10)/2
NTplts = 4000

5. Now, create a grid of six phase space plots, one for each of the following values of betatron tune: \(\nu\) = 20.50291, 20.580902, 20.6672073, 20.7498, 20.800776, 20.99091. For each plot, track a collection of particles with initial conditions \(u_0\) = 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, \(v_0\) = 0, and track each particle for 4000 turns. In each phase space plot, indicate with a red circle the area which would contain 95% of the particles in the beam (in the absence of the beam-beam kick itself).

par(mfrow=c(2,3))
kpl = 0
while(kpl<length(ktune)){
   kpl=kpl+1
   alim = 10
   nu = ktune[kpl]
   cn = cos(2*pi*nu)
   sn = sin(2*pi*nu)
   plot(0,0,typ="n",xlim=c(-alim,alim),ylim=c(-alim,alim),asp=1,xlab="u",ylab="v",
      main=paste("nu =",round(nu,4)))
   symbols(0,0,circles=sqrt(6),add=TRUE,fg="red",inches = FALSE,lwd=3)
   m=0
   while(m<length(u0)){
      m=m+1
      u = u0[m]
      v = 0
      i=0
      while(i<NTplts){
         points(u,v,pch=".")
         i=i+1
         u1 =  u
         v1 =  v + dv(u)
         u  =  cn*u1 + sn*v1
         v  = -sn*u1 + cn*v1
      }
   }
}

---

1. Northern Illinois University and Fermi National Accelerator Laboratory↩