columns=fullflexible, %this is to allow simple copy-paste of the code to a zgoubi input file

keywordstyle=\color{red},

commentstyle=\color{listinggray},

otherkeywords={?}

%morecomment=[l]{!\ }% Comment only with space after !

}

\usepackage[section]{placeins}% to be able to use \FloatBarrier

\usepackage[capitalise]{cleveref}%to use \cref

\usepackage{longtable}

\addtolength{\textwidth}{1.5 cm}

\addtolength{\hoffset}{-0.8 cm}

%\addtolength{\topmargin}{-2 cm}

%\addtolength{\textheight}{3. cm}

%opening

\title{Zgoubi, would you do for me\\TRIUMF's High Resolution Separator?}

\author{Thomas Planche}

%\footnote{tplanche@triumf.ca}\\~\\ TRIUMF\thanks{This work has been supported by the Natural Sciences and Engineering Research Council of Canada. TRIUMF also receives federal funding via a contribution agreement through the National Research Council of Canada.}}

\begin{document}

\maketitle

\tableofcontents

\section{Basic Layout}\label{sec:layout}

The optical system we will consider in this tutorial is designed to separate rare isotopes with mass/charge differences of only one part in 20\,000 in beams with transverse emittances of at least 3\,$\mu$m\footnote{Un-normalized emittance. Note that to achieve such resolution an energy spread of the order of 1\,eV or less is required.}~\cite{TRI-DN-16-09}.

It is composed of:

\begin{itemize}

\item a source slit;

\item followed by an 80\,cm drift;

\item followed by two identical 90\,deg.~magnetic dipoles (bending radius~=~120\,cm, edge angle~$\approx$~26.5\,deg.) separated by a 2$\times$80\,cm drift;

\caption{TRIUMF high resolution separator basic layout (from source slit to image slit).}\label{fig:marco}

\end{figure}

\FloatBarrier

The basic parameters of this high resolution separator (HRS) were determined using the linear optics code {\tt TRANSOPTR}~\cite{optrOnGitlab}.

To achieve in practice the desired resolving power, it is essential to compensate non-linear aberration. For this purpose an electrostatic multipole corrector is added halfway between the two dipoles. The detailed design of the most critical components -- the dipoles and the multipole corrector -- was accomplished using a 3D finite element code ({\tt OPERA}) and non-linear optics codes such as {\tt COSY-INFINITY} and {\tt zgoubi}.

In this note I will go through the steps of the design work done using {\tt zgoubi}.

\section{Define your {\tt'OBJECT'}}

Let's start with something simple and try to track one single particle. Let's say we want to track a single 60\,keV $\rm{^{238}U}^+$ ion. At first we will track it through magnetic elements only, so the knowledge of the magnetic rigidity $B\rho$ is sufficient. As a reminder the magnetic rigidity is given by:

\begin{equation}

B\rho=\frac{p}{q}\,,

\end{equation}

where $q$ is the charge of the particle, and $p$ its momentum given by:

\begin{equation}

p^2c^2=E^2-m^2c^4\,,

\end{equation}

where $E$ is the particle's total energy, $m$ its mass, and $c$ the speed of light. For non-relativistic particles, like our 60\,keV uranium ion, the momentum can also be calculated from:

\begin{equation}

p=\sqrt{2mqV}\,,

\end{equation}

where $V$ is the beam potential (60\,kV in our case).%\footnote{$qV$ is often called `kinetic' energy. It is a true kinetic energy only in the non-relativistic approximation. In the relativistic case the beam potential is given by $(\gamma-1)mc^2$ and cannot be seen as a kinetic energy term, whatever way you look at it~\cite{2014PHYS-1}.}.

If you look into {\tt zgoubi} user's guide~\cite{meot2012zgoubi} you will find you need to call use the keyword {\tt 'OBJECT'} (or {\tt 'OBJET'} if you like to talk to {\tt zgoubi} in French). There are several ways (KOBJ=1 to 6) to have the 'object' you will track be one single particle. Let's try to use KOBJ=2:

\begin{lstlisting}[label=list:OBJET]

zgoubi HRS tuto: DIPOLE !zgoubi requires input files to be given a title on line 1

'OBJET'

?? !BORO[kg.cm]: Brho of 60 keV 238U+ in kG.cm

2 !KOBJ=2: initial coordinates must be entered explicitly

1 1 !total number of particles; number of distinct momenta

1 !1 or -9 (-9 disables the tracking of this particle)

\end{lstlisting}

{\it Assignment: copy-paste these lines into a text file that will serve as input to {\tt zgoubi}. Replace the red double question mark `{\color{red}??}' by the appropriate value of {\tt BRHO}.}

\caption{Parameters used to define the field map and geometrical boundaries, from {\tt zgoubi} user's guide~\cite{meot2012zgoubi}.}\label{fig:DIPOLE}

\end{figure}

We will first try to simulate our separator using {\tt zgoubi}'s built-in {\tt 'DIPOLE'} model. \cref{fig:DIPOLE} reproduces a figure from {\tt zgoubi} user's guide~\cite{meot2012zgoubi} that is very useful to understand the input values to {\tt 'DIPOLE'}.

To mimic our separator dipole, we will need to determine the parameters such that the beam is bent 90\,deg.~and passes through the center of the dipole at a radius of 120\,cm. To do that it is convenient to stop the tracking right in the middle of the dipole:

\FloatBarrier

\begin{lstlisting}[label=list:halfDIPOLE]

'DIPOLE'

2 !IL= 0 no outpot, = 2 output trajectory to zgoubi.plt, etc.

65 120. !AT[deg.], RM[cm]

65 ?? 0. 0. 0. !ACENT[deg.];B0[kG]; N; B; GX

7.0 0. !ENTRANCE FIELD BOUNDARY: fringe field extend[cm]; unused

Note that the tracking starts 20\,deg. before the effective edge of the magnet to leave room for the fringe field\footnote{About fringe fields: in most cases it is unimportant to know the precise shape of the field field. For dipoles (and quadrupoles~\cite{baartman1997intrinsic}) the lowest order aberration (second order for dipoles, third order for quads) is practically independent of the details of the fringe field. The next order aberration depends only on the fringe field extend. For more details read Ref.~\cite{baartman2001end}. It is only when you look at quite-high-order aberrations that the shape of the field fall off starts to matter.

For our dipole I chose a fringe field as simple as possible: fringe field extend~=~magnet gap (7\,cm) and a shape defined by a single Enge coefficient: {\tt C1}=1.8 (all the others are set to 0).} to fall off to (practically) zero, which leads to {\tt AT}=90/2+20=65\,degree.

{\it Assignment: copy-paste these lines right bellow the {\tt 'OBJECT'} definition (\cref{list:OBJET}). Guess the values of {\tt B0, RE}, and {\tt TE} such that the particle arrives at the magnet center:

\begin{itemize}

\item bent 45\, deg.;

\item at {\tt Y}=120\,cm (radial coordinate in the local reference frame of {\tt 'DIPOLE'});

\item with {\tt T}=0\,mrad (angle w.r.t.~the reference trajectory).

\end{itemize}

To check how good your guess was look into the output of {\tt 'FAISCNL'} (in this case: in the file `atDIPOLEcenter.fai'). You can also plot the step-by-step output into zgoubi.plt (triggered by IL=2). If you use {\tt gnuplot}, try something like:}

\begin{verbatim}

set xlabel "Theta/deg."

set ylabel "R/cm"

plot "zgoubi.plt" using ($22*180./pi):($10)

\end{verbatim}

\section{{\tt 'FIT'} dipole parameters}

To accurately (and rapidly) determined the values of {\tt B0}, and {\tt RE} we will use the {\tt 'FIT'} keyword. Add this right before the {\tt 'END'} statement:

{\it Assignment: Find values {\tt B0}, and {\tt RE} from fit. run {\tt zgoubi} with those fitted values and plot the resulting particle trajectory (from zgoubi.plt). You should get plots like in~\cref{fig:DIPOLEtraj}. }

\caption{Trajectory output into zgoubi.plt plotted in cylindrical (left) and Cartesian (right) coordinates. }\label{fig:DIPOLEtraj}

\end{figure}\FloatBarrier

Side note: looking for the a parameter number can be sometime quite frustrating, especially for an element like {\tt 'DIPOLE'} that has a large number of input parameters. What works best for me it to look for it by trial-and-error, watching the `INITIAL' value column in {\tt zgoubi}'s console output.

\section{{\tt'MATRIX'} to fit edge angles}

The dipole edge angles are set so that the phase advance from the source slit to the image slit is exactly 180\,deg. (in the horizontal direction). One could fit the transfer matrix coefficient $m_{12}$ and $m_{21}$ to be zero at the image slit. One can also fit the transfer matrix coefficient $m_{11}$ and $m_{22}$ to be zero half-way through the HRS, since the phase advance from the source slit should there be 180/2=90\,degree. I choose to follow the second approach.

{\it Assignment: we have only simulated one half of a magnet so far. Let's set up the other half: in {\tt 'DIPOLE'} set {\tt AT}=65$\times$2=130\,deg., and define the exit face parameters so that the exit face is the mirror image of the entrance face. Be careful with the signs! They are NOT consistent with the standard {\tt TRANSPORT} definition (see~\cref{fig:DIPOLE}) }

Before and after the dipole we should also have 80\,cm drift length (see~\cref{sec:layout}). The thing is that we already have a drift length built into the dipole definition: the 2$\times$20 extra degree we left for the field to fall off. So the drift length we have to add is $80-120\tan(20\,{\rm deg.)}$. So add this before and after your dipole definition:

\begin{lstlisting}[label=list:DRIFT]

'DRIFT'

36.32357 ! drift length=80-120*tan(20/deg.)

\end{lstlisting}

Now to we get {\tt zgoubi} to calculate transfer matrix coefficient: use the keyword {\tt 'MATRIX'} together with an {\tt 'OBJECT'} defined using {\tt KOBJ}=5 (see user's guide~\cite{meot2012zgoubi}).

The fit should converge to an value for the edge angle close to 26.56\,deg.\footnote{26.56\,deg. is the number that the linear optics code {\tt TRANSOPTR} finds}.

\section{Multiparticle simulation}

I got from my good friend Jim (who simulated the HRS using {\tt COSY-INFINITY}) a file containing the coordinates of 3\,000 randomly generated particles. You can download the file from here: \href{http://beamphys.triumf.ca/~tplanche/text/designs/HRS/zgoubi/DIPOLE/noEspread.beam}{link}. Each line of this file contains is (using {\tt zgoubi}'s notation): {\tt Y,T,Z,P,S,DP}. This format is what {\tt 'OBJECT'} can read when {\tt KOBJ}=3.01:

Jim warned me that in this file the vertical emittance is 6\,$\mu$m (while the horizontal one is 3\,$\mu$m), and that the aspect ration between vertical size and vertical angle was wrong by a factor of 4: to fix the aspect ratio and bring the vertical emittance down to 3\,$\mu$m I use the scaling factors on line 5 of the {\tt 'OBJECT'} input.

{\it Assignment: write an input file to simulate the entire HRS from source slit to image slit. Track the 3\,000 particles in noEspread.beam and same their final coordinates to a file. Now change the mass of the particle by a factor 1/20\,000 (i.e. change its {\tt BRHO} by a factor 1/10\,000), track gaian and save the final coordinates into another file. Plot the horizontal phase space coordinates from both files, you should get a plot like~\cref{fig:twomassesInit}:}

\caption{Two isobars with a relative mass difference of 1/20\,000 through the HRS, before correcting non-linear aberration.}\label{fig:twomassesInit}

\end{figure}

What you see in~\cref{fig:twomassesInit} is that the two masses completely overlap. This is (mostly) the effect of second order (sextupole) aberration. To achieve the designed resolution, we must correct the aberration.

\section{Dipole edge curvature}

To correct the second order aberration let's try to add curvature to the edges of our dipoles. Once again, we would like to adjust the edge curvature using a {\tt 'FIT'}, but what should the objective of that fit be?