Homework 3 Discussion

Some General Comments:

• grade: 30 possible
  + 10 pts each problem
  + ave score: 20/30

• still: last names in your file names – almost there!
  
• computers and determining angles: beta > 0 always!
  
• tune shift
  
• some still need practice: W vs E vs pc vs mc^2 vs beta vs gamma …
  
• MAD-X: mux = psi_x/2pi; Q1 (summary table) = periodic mux
• periodic CS parameters vs. non-periodic values

• creating matrices, arrays, data frames, etc. – always check what’s going on; make sure it’s what you want

Betatron Motion

          [,1]       [,2]
[1,] 1.0237788 -28.970214
[2,] 0.1007168  -1.873243

Suppose a beam of particles entering a synchrotron has a phase space distribution (consider \(x\) only, for now) that is parameterized by the Courant-Snyder variables \(\beta\) = 25 m, \(\alpha\) = -1 and rms emittance \(\epsilon\) = 3.75 \(\pi\) mm-mrad when the particles first pass by a particular observation point in the ring. If a matrix calculation were performed, corresponding to the magnetic elements about the circumference of this ring, we find that the one-turn matrix at this same observation point is given by \[ M_0 = \begin{pmatrix} 1.0238 & -28.970~{\rm m} \\ 0.10072{\rm /m}& -1.8732 \end{pmatrix} \] which transports the vector \(\vec{X} = (x,x')^T\) according to \(\vec{X}_{n+1} = M_0\vec{X}_n = M_0^n\vec{X}_0\) where \(n\) is the revolution number, or turn number.

1. Create a distribution of several thousand particles that has the initial Courant-Snyder parameterization given above. Plot the phase space distribution after \(n\) = 0, 1, 2, 3, 4 and 5 revolutions, as well as histograms of the distribution in \(x\) (corresponding to possible beam width measurements at our observation point) for the same values of \(n\).

beta = 25     # m
alfa = -1  
gamm = (1+alfa^2)/beta
eps  = 3.75   # pi mm-mrad
M0 = matrix(c(1.0238, -28.970, 0.10072, -1.8732),ncol=2,byrow=TRUE)

Npar = 5000
x0  =  rnorm(Npar,0,sqrt(eps*beta))  # mm
xp0 = (rnorm(Npar,0,sqrt(eps*beta)) - alfa*x0)/beta  # mrad
# save; will use these values later...

x  = x0
xp = xp0
par(mfrow=c(2,2))
i = 0
while(i<6){
   i=i+1
   plot(x,xp,xlim=c(-1,1)*40,ylim=c(-1,1)*5,pch=".")
   hist(x, xlim=c(-1,1)*40)
   j = 0
   while(j<Npar){
      j = j+1
      U = M0 %*% c(x[j],xp[j])
      x[j]  = U[1]
      xp[j] = U[2]
   }
}

2. Create a plot of how the beam distribution rms width varies with revolution number, for \(0<n<200\).

beta = 25     # m
alfa = -1  
gamm = (1+alfa^2)/beta
eps  = 3.75   # pi mm-mrad
M0 = matrix(c(1.0238, -28.970, 0.10072, -1.8732),ncol=2,byrow=TRUE)

Npar = 5000
x  =  x0  # mm
xp = (x0 - alfa*x0)/beta  # mrad

xrms = sd(x)
i = 0
while(i<200){
   i=i+1
   j = 0
   while(j<Npar){
      j = j+1
      U = M0 %*% c(x[j],xp[j])
      x[j]  = U[1]
      xp[j] = U[2]
   }
   xrms[i] = sd(x)
}
plot(xrms, ylim=c(0,15), typ="l")

3. Next, using the matrix \(M_0\), find the periodic Courant-Snyder parameters \(\beta_0\) and \(\alpha_0\) at the observation point, as well as the betatron tune \(\nu_x\) of the synchrotron.

mu0 = acos(sum(diag(M0))/2)
#  NOTE:  if M[1,2]<0, then sin(mu)<0 -- beta>0 always!
if( M0[1,2]<0 ) mu0 = 2*pi - mu0

beta0  =  M0[1,2]/sin(mu0)
alpha0 = (M0[1,1]-M0[2,2])/2/sin(mu0)
mu0; beta0; alpha0

[1] 4.273758
[1] 31.99924
[1] -1.599962

nu = mu0/2/pi   # modulo 2pi
nu

[1] 0.6801898

4. Can you relate the period of the rms width variation to the betatron tune?

# Let's do a trial and error exercise...
plot(xrms, typ="l", xlim=c(0,50) )
xplt = c(1:200)
y = mean(xrms)*(1+.5*sin(2*pi*(2*nu)*(xplt)))
points(y,col="blue")

[ Looks like xrms oscillates at 2x the tune.]

5. Create phase space plots of (\(\beta_0\;x'+\alpha_0\;x\)) vs. \(x\) of the distribution for \(n\) = 0, 1, 2, 3, 4, 5.

x = x0
p = alpha0*x0 + beta0*xp0


R0 = matrix(c(cos(mu0),sin(mu0),-sin(mu0),cos(mu0)),ncol=2,byrow=TRUE)
i = 0
par(mfrow=c(2,3))
while(i<6){
   i=i+1
   plot(x,p,xlim=c(-1,1)*80,asp=1,pch=".")
   j = 0
   while(j<Npar){
      j = j+1
      U = R0 %*% c(x[j],p[j])
      x[j] = U[1]
      p[j] = U[2]
   }
}

6. Change \(\alpha\) and \(\beta\) of the initial distribution to be equal to \(\alpha_0\) and \(\beta_0\), thus creating a new corresponding initial particle distribution; repeat the plots of part (e) above.

x0  =  rnorm(Npar,0,sqrt(eps*beta0))  # mm
xp0 = (rnorm(Npar,0,sqrt(eps*beta0)) - alfa0*x0)/beta0  # mrad

x = x0
p = alpha0*x0 + beta0*xp0


R0 = matrix(c(cos(mu0),sin(mu0),-sin(mu0),cos(mu0)),ncol=2,byrow=TRUE)
i = 0
par(mfrow=c(2,3))
while(i<6){
   i=i+1
   plot(x,p,xlim=c(-1,1)*80,asp=1,pch=".")
   j = 0
   while(j<Npar){
      j = j+1
      U = R0 %*% c(x[j],p[j])
      x[j] = U[1]
      p[j] = U[2]
   }
}

As a check:

x = x0  # mm
p = alpha0*x0 + beta0*xp0  # mm

xrms = sd(x)
i = 0
while(i<200){
   i=i+1
   j = 0
   while(j<Npar){
      j = j+1
      U = R0 %*% c(x[j],p[j])
      x[j] = U[1]
      p[j] = U[2]
   }
   xrms[i] = sd(x)
}
plot(xrms, ylim=c(0,15), typ="l")

[So, in this case the beam is well matched to the periodic CS parameters of the ring, and hence – while each particle oscillates in phase space – the overall width of the distribution does not change significantly turn-by-turn.]

The Standard FODO Cell

Note: a MAD-X file associated with this Homework Problem can be found here.

1. Provide a formula for the Courant-Snyder parameters \(\beta_{max,min}\) for a thin lens FODO cell in terms of the cell phase advance, \(\mu\), and the half-cell length \(L\).

\[ \beta_{\pm} = \frac{L}{\sin\mu/2}\sqrt{\frac{1\pm\sin\mu/2}{1\mp\sin\mu/2}} \]

2. Make a single plot of \(\beta_{max}/L\) and \(\beta_{min}/L\) vs. \(\mu\) for \(0 < \mu < 150^\circ\).

bmaxonL = function(x){sqrt((1+sin(x*pi/180/2)) / (1-sin(x*pi/180/2)))/sin(x*pi/180/2)}
bminonL = function(x){sqrt((1-sin(x*pi/180/2)) / (1+sin(x*pi/180/2)))/sin(x*pi/180/2)}
curve(bmaxonL,0,150,ylim=c(0,20),xlab="Cell Phase Advance [deg.]", ylab="betaMax,Min / L")
curve(bminonL,add=TRUE,lty=2,col="blue")

3. Compute \(\beta_{max}\) and \(\beta_{min}\) for the Tevatron FODO cell, where the lenses were separated by 30 m, using a cell phase advance of 75 degrees. What focal length does this correspond to?

betaMax = bmaxonL(75)*30
betaMin = bminonL(75)*30
F = 30/2/sin(75*pi/180/2)
betaMax;betaMin # meters

[1] 99.93074
[1] 24.3024

F  # meters

[1] 24.64019

4. Use MADX to perform the same calculation at the downstream end of each quadrupole in a thick lens system. Use Tevatron-like quadrupole magnets with effective lengths of 2 m. Note that the total cell length should be equivalent to that of the “thin lens” version and that the total cell phase advance should again be 75 degrees. What values for \(\alpha_x\) and \(\alpha_y\) do you find at these locations? What value of \(|K| = |B'|/B\rho\) is required in this calculation? What focal length does this value of \(K\) correspond to?

We generate a model as follows:

LQ = 2    ;  ! quad length
F  = 25   ;  ! quad focal length

L0 : drift, L=28 ;
QF : quadrupole, L=LQ, K1= 1/F/LQ ;
QD : quadrupole, L=LQ, K1=-1/F/LQ ;

Tcell : line=(L0, QD, L0, QF) ;

Using a MATCH command in MAD-X we can generate the appropriate 75 degree phase advance, and perform a TWISS calculation.

match, sequence=Tcell ;
    constraint, sequence=Tcell, range=#e, mux=75/360, muy=75/360 ;
     vary, name=QF->k1,  lower = -1/25, upper=1/25, step=5.0e-3 ;
     vary, name=QD->k1,  lower = -1/25, upper=1/25, step=5.0e-3 ;
    simplex, calls=5000, tolerance=1.0e-6;
endmatch ;

From the output of calculation, \((\alpha_x,\alpha_y)\) = (2.022264919, -0.5561859735) at the output of the focusing quad, and = (-0.556511509, 2.022216597) at the output of the defocusing quad. The output of the match operation yields for \(K_1\) of the quadrupoles:

qf->k1             =      0.02075938165 ;
qd->k1             =     -0.02076360524 ;
1/(qf->k1*qf->l)   =        24.08549582 ;

The last line is the effective focal length of the magnet.

5. Re-compute the values of \(\beta_{max}\) and \(\beta_{min}\) at the mid-point of each thick quad, using MADX. Also, what are the values of \(\alpha\) at these points?

Here, we create “half-length” quadrupoles and generate a new beam line:

K1fit = QF->K1 ;

QF2 : quadrupole, L=LQ/2, K1= K1fit ;
QD2 : quadrupole, L=LQ/2, K1=-K1fit ;

Tcell2 : line=(QF2, L0, 2*QD2, L0, QF2) ;

The output will now show the values of the CS parameters at the midpoints of the magnets; we see that \((\alpha_x,\alpha_y)\) = (0,0) (at least to the \(10^{-16}\)-level) at both midpoint locations for this symmetric system, as these points are where the values of \(\beta\) reach their maxima. Note also that these maximum betas have values:

	thin lens	MAD-X 2-m quads
\(\beta_{max}\)	99.93074	99.25046839
\(\beta_{min}\)	24.3024	24.46556152

illustrating that for this calculation, the thin lens approximation is good to better than 1%.

6. In the thin lens calculations, how would you determine CS parameters at the quad mid-points?

[ Change the matrix computation from… \[ \begin{pmatrix} 1 & L \\ 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ q & 1 \end{pmatrix} \begin{pmatrix} 1 & L \\ 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ -q & 1 \end{pmatrix} \] to \[ \begin{pmatrix} 1 & 0 \\ -q/2 & 1 \end{pmatrix} \begin{pmatrix} 1 & L \\ 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ q & 1 \end{pmatrix} \begin{pmatrix} 1 & L \\ 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ -q/2 & 1 \end{pmatrix} \] etc. ]

THE Booster Ring

Note: a MAD-X file associated with this Homework Problem can be found here.

The original Fermilab Booster synchrotron – still in operation today – was designed in the late 1960’s and has the following configuration:

• 24 identical periods
  
• each period contains two BF magnets and two BD magnets, each of which both bends and focuses (or defocuses) the beam
• each BF magnet and each BD magnet is 2.889612 m long (design value), and are rectangular, rather than sector, elements
  
• each BF magnet operates with a central field value of 0.72588 T, and has a field gradient of 1.60761 T/m
  
• each BD magnet operates with a central field value of 0.61727 T, and has a field gradient of -1.71101 T/m
  
• the fields above correspond to a final proton kinetic energy of 8 GeV
• each of the 24 periods has the following geometric configuration:
  • (0.6 m) – BF – (0.5 m) – BD – (6.0 m) – BD – (0.5 m) – BF – (0.6 m)

1. Verify that the above configuration should bend the beam through 360 degrees. What is the total circumference of the ring, and what fraction of the circumference contains bending?

# Total bending length x field = (BL)_total should equal 2pi*Brho
Lmag = 2.889612
Bf = 0.72588
Bd = 0.61727
BLtot = 24*(Bf + Bd)*2*Lmag  # T-m
mc2 = 0.938  # GeV
W0  = 8      # GeV/c
bg  = sqrt(((8+mc2)/mc2)^2-1)
Brho = mc2*bg*10/2.9979  # T-m

BLtot/Brho - 2*pi

[1] 0.0001016444

C = 24*(4*2.889612 + 2*0.5 + 2*0.6 + 6)
C

[1] 474.2028

pkf = (24*4*2.889612)/C
pkf

[1] 0.5849876

2. Create an appropriate MADX input file with a beam line that describes the Fermilab Booster. Compute and print out the values of the Courant-Snyder (Twiss) parameters at the end of each element according to the geometry described above.

Some relevant lines from the input file:

! The magnets:
BD: rbend, L= 2.889612, angle= 0.61727*2.889612/brho,  k1= -1.71101/brho ;
BF: rbend, L= 2.889612, angle= 0.72588*2.889612/brho,  k1=  1.60761/brho ;

! Drift spaces:
O  :   drift,  L= 0.6 ;
OO :   drift,  L= 0.5 ;
OOO:   drift,  L= 6.0 ;

! The layout:
BCEL:  line=(     O, BF, OO, BD, OOO, BD, OO, BF, O )   ;

and a table from the output:

* NAME                   S           BETX           BETY             DX               MUX              MUY
"BCEL$START"             0    33.67385206    5.271002857    3.204695719                 0                0
 "O"                   0.6    33.68454285    5.339301055    3.204695719    0.002835519438    0.01803901468
 "BF"          3.490214615     20.8840135    10.82233146    2.604892322     0.01891236161    0.08519977287
 "OO"          3.990214615    17.30501983    12.99977669     2.40508385     0.02309882393    0.09191082297
 "BD"           6.88026237    7.587521356    20.46123159    1.850080196     0.06707513392      0.118018933
 "OOO"         12.88026237    7.587521356    20.46123159    1.850080196      0.2122354719     0.1653852637
 "BD"          15.77031013    17.30501983    12.99977669     2.40508385      0.2562117819     0.1914933737
 "OO"          16.27031013     20.8840135    10.82233146    2.604892322      0.2603982442     0.1982044238
 "BF"          19.16052474    33.68454285    5.339301055    3.204695719      0.2764750864      0.265365182
 "O"           19.76052474    33.67385206    5.271002857    3.204695719      0.2793106058     0.2834041966
 "BCEL$END"    19.76052474    33.67385206    5.271002857    3.204695719      0.2793106058     0.2834041966

3. What are the betatron oscillation tunes \((\nu_x, \nu_y)\) of the synchrotron?

The output of the TWISS command will generate values for “Q1” and “Q2” which are the horizontal and vertical tunes (here, through one period); so:

value, 24*TABLE(SUMM,Q1);    ! note:  here used one period, hence x24
24*table( summ q1 ) =        6.703454539 ;
value, 24*TABLE(SUMM,Q2);
24*table( summ q2 ) =        6.801700719 ;
nux0 = 24*Table(summ,Q1) ;
nuy0 = 24*Table(summ,Q2) ;

Also, we can see that MUX and MUY for one period is given in the above table. These numbers are \(\Delta\psi/2\pi\) and so the tunes will be 24 times those values.

4. Plot the periodic amplitude functions \(\beta_x\), \(\beta_y\) and the dispersion function \(D_x\) throughout a single period.

Booster Period

5. Protons are injected into the Booster ring from a linear accelerator, when the proton kinetic energy is 400 MeV. Make a plot of the revolution time (in microseconds) for an ideal proton as a function of magnetic field corresponding to the ring’s range of operation.

bet  = function(x){ sqrt(1-(mc2/(mc2+x))^2) }
tprd = function(x){ C/2.9979e8/bet(x) }
curve(tprd(x)*1e6,0.4,8,
   xlab="Kinetic Energy [GeV]", ylab="Rev. Time [microsec.]")

ke = function(bg){mc2*(sqrt(1+bg^2) - 1)}
bgMin = sqrt(((mc2+0.4)/mc2)^2-1)
bgMax = sqrt(((mc2+8.0)/mc2)^2-1)
curve(tprd(ke(x*bgMax))*1e6,bgMin/bgMax,
   xlab="Fraction of maximum field", ylab="Rev. Time [microsec.]")

6. Suppose that a single BF magnet in the ring has its gradient altered by +5%. Compute the resulting change in the horizontal and vertical tunes of the Booster using appropriate modifications to the matrices found above. Estimate the change by using the simple tune shift formula and compare.

Here, we take a single BF magnet (BFerr) and increase its gradient by 5%, then build a new ring description that has one period with this magnet in it, and 23 periods with “normal” magnets :

BFerr: rbend, L= 2.889612, angle= 0.72588*2.889612/brho,  k1=  1.60761/brho*1.05 ;
BCELerr:  line=(     O, BF, OO, BD, OOO, BD, OO, BFerr, O )   ;
BRing:   line=( 11*(BCEL), BCELerr, 12*(BCEL) ) ;

We then re-compute the Twiss paramters and extract the new tunes (Q1 and Q2):

value, TABLE(SUMM,Q1);
table( summ q1 )   =         6.72132284  ;
value, TABLE(SUMM,Q2);
table( summ q2 )   =         6.797235611 ;
nux1 = Table(summ,Q1) ;   ! note:  here used complete ring
nuy1 = Table(summ,Q2) ;

We see that the new tunes differ from the original tunes, due to the error field, by:

nux1-nux0          =      0.01786830108  ;
nuy1-nuy0          =     -0.004465107912 ;

We now use our “tune shift formula” to do a quick check:

\[ \Delta\nu_{x,y} = \pm \frac{1}{4\pi}\beta_{x,y}\Delta K\ell \]

Bpf = 1.60761   # T/m for the BF magnet
betaxBF = 26    # from the MAD output; use average values within the BFerr magnet...
betayBF = 8.0 
dnux =  1/4/pi*betaxBF*(Bpf/Brho*0.05)*Lmag # "+" since   focuses in x
dnuy = -1/4/pi*betayBF*(Bpf/Brho*0.05)*Lmag # "-" since defocuses in y
dnux; dnuy

[1] 0.01620822
[1] -0.004987145

---

1. Northern Illinois University and Fermi National Accelerator Laboratory↩