%
% !BIB TS-program = biber
% !BIB program = biber
% !TeX spellcheck = en_US
% !TeX program = pdflatex
%
% support document for version: v 3.00  Feb 2026   Volker RW Schaa
%   
%
\documentclass[%              % paper size definition not needed as JACoW page size is set
               %boxit,        % show page boundaries to check whether paper is inside correct margins
               %refpage       % references placed on separate page
               %biblatex,     % BibLaTeX is to be used
               %nospread,     % flushend option: do not fill with whitespace to balance columns
               %hyphens,      % allow \url to hyphenate at "-" (hyphens)
               %luatex,       % use LuaLaTeX to process the file
               ]{jacow}

\usepackage[english]{babel}

\listfiles

%%
%% Lengths for the spaces in the title
%%   \setlength\titleblockstartskip{..}  %before title, default 3pt
%%   \setlength\titleblockmiddleskip{..} %between title + author, default 1em
%%   \setlength\titleblockendskip{..}    %afterauthor, default 1em

\newcommand{\B}{\mathbf B}

\begin{document}

%%%%% Title - Author block generated from indico contribution details %%%%%
%%%%%%%%%%%%%%%%%%%% Do not modify  %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% If you need to modify the content of this block, do it by modifying the  %%%%%
\title{MATCHING OF CRYSTALLINE BEAMS INTO THE BNL BOOSTER\thanks{Work supported by Brookhaven Science Associates, LLC under Contract No. DE-SC0012704 with the U.S. Department of Energy.}}

\author{S.J.~Brooks\thanks{sbrooks@bnl.gov}, Brookhaven National Laboratory, Upton, NY, USA}

\maketitle

%%%%% End of Title - Author block %%%%%

%%%%% You can modify the abstract to add LaTeX commands. However, please make sure that it is similar to the one in indico %%%%%
\begin{abstract}
A method has been developed to match large three-dimensional Coulomb crystals into existing accelerators.  Simulations show bunch temperatures below 1 Kelvin persisting for thousands of turns in the BNL Booster and for tens of thousands of turns in an idealised ring (both with no cooling).  Focussing forces balance space charge to give zero net phase advance and the crystal structure can be preserved for tens of milliseconds.  Magnetic focussing confines all three dimensions of the bunch without RF, since the bunch rotation in the ring plane is not at the ring revolution frequency.  A zero temperature, uniformly-filled ellipsoid of charge gives a 15-parameter model that can be used for matching and linear stability analysis.
\end{abstract}

\section{Booster Field Model}
The calculations that follow require a 3D Maxwellian field model of the booster that is free from discontinuities near the beam orbit.  The basic superperiod lattice from \cite{CourantBoosterLattice} was implemented, which contains identical bending dipoles 2.4\,m long and F and D quadrupoles 0.50375\,m long.  

In a local coordinate frame where (0,0,0) is the centre of the dipole arc, $z$ is the beam direction and $y$ is vertical, the dipole field is given by calculating
\begin{eqnarray*}
%R&=&L/\theta \\
t_y&=&\tan(y/f) \\
t_{z\pm}&=&\tanh((\pm z\cos\tfrac\theta2+(\tfrac{L}{\theta}-x)\sin\tfrac\theta2)/f) \\
t_{r\pm}&=&\tfrac12t_{z\pm}(1+t_y^2)/(1+t_y^2t_{z\pm}^2) \\
t_{i\pm}&=&\tfrac12t_y(1-t_{z\pm}^2)/(1+t_y^2t_{z\pm}^2) \\
\B_1 &=& (-(t_{i+}+t_{i-})\sin\tfrac\theta2,t_{r+}+t_{r-},(t_{i+}-t_{i-})\cos\tfrac\theta2),
\end{eqnarray*}
where $L,\theta$ are the dipole length and angle and $f$=0.15\,m determines the fringe field length.  The fringe field is shaped like $\tfrac12+\tfrac12\tanh(z/f)$ and the formulae for $t_{r,i\pm}$ come from the complex components of $\tanh(z+iy)$.  $\B_1$ is the field of a unit strength dipole, in general it is scaled to $\B=B_\mathrm{dip}\B_1$.

A quadrupole with soft ends does not have an analytic formula, so the series expansion in \cite{BrooksOffAxis} is used with a fringe field function again shaped like $\tfrac12+\tfrac12\tanh(z/f)$.  For the quadrupoles, $f$=0.05\,m and the series is calculated up to and including the $r^{10}$ term.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{WEP4324_f1}\qquad%
   \caption{Booster field map with successive zooms (top, bottom right), the small magnet is a quadrupole.}
   \label{fig:boosterfieldmap}
\end{figure}

Placing these field models and adding them together, with appropriate transverse extent limits that the beam does not encounter, gives the field illustrated in Fig.~\ref{fig:boosterfieldmap}.

\section{Matching Method}
Repulsive space charge forces dominate the dynamics of Coulomb crystals, while velocity spread and emittance are near zero.  Thus, the bunch is modelled as a uniform ellipsoid of charge with zero temperature, parametrised as shown in Fig.~\ref{fig:ellipsoidmodel}.  The ellipsoid is defined by its centroid $\mathbf s=\langle(\mathbf x,\mathbf v)\rangle$, covariance matrix $C=\langle\delta\mathbf x\delta\mathbf x^T\rangle$ and velocity distribution $V$ where $\delta\mathbf v=V\delta\mathbf x$.  Space charge forces are linear within a uniform charged ellipsoid \cite{LinearA,LinearB,LinearC}, so this form is preserved if external forces are also linear.  The centroid $\mathbf s$ evolves as a single particle while the 15 parameters $(C,V)$ have time derivatives calculated in \cite{firstordermodel}.  This model does not capture higher-order vibrational motion or phonons.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{WEP4324_f2}
   \caption{Parameters of the ellipsoid bunch model.}
   \label{fig:ellipsoidmodel}
\end{figure}

A matched solution for a ring would have $(\mathbf s,C,V)_\mathrm{in}=(\mathbf s,C,V)_\mathrm{out}$ over a single turn.  In this paper $(\mathbf s,C,V)$ are integrated with 4$^\mathrm{th}$ order Runge--Kutta steps and $\delta t = (8\,\mathrm{mm})/v_\mathrm{bunch}$.  The centroid dynamics are independent, so the closed orbit with $\mathbf s_\mathrm{in}=\mathbf s_\mathrm{out}$ can be found fairly standardly using multi-variable Newton iteration, taking care to find $t_\mathrm{out}$ that stops the particle exactly at the end of the turn.  

A similar iteration with 15 variables can also solve $(C,V)_\mathrm{in}=(C,V)_\mathrm{out}$ and the first derivatives of the mapping are used to obtain an eigenvalue stability analysis \cite{BrooksCOOL25}.  There are, however, some computational and physics issues.

If numerical derivatives are used, step sizes in each variable must be chosen carefully as $C$ and $V$ have different units.  A slightly enhanced Newton's method \cite{BrooksSVD,BrooksOptimiser} is used here, which also evaluates the mean square error, so units for the parameter and goal vectors are chosen to make $C$ and $V$ similar magnitude.

A second issue is that for some lattices, $V$ tends to infinity after a few cells because one of the ellipsoid's axes tries to invert.  Unintuitively, the model passes through infinite charge density at this time.  This can occur because although shrinking a 3D charge distribution to a 0D point or 1D line takes infinite energy, a 2D surface can be obtained with finite energy.  In reality, this would lead to semi-random pairwise scattering between the ions in the crystal, so is not desired.  The solution is to reduce the focussing in the lattice until the ellipsoid no longer inverts (changing the initial size for a given charge may also help).  This leads to surprisingly low levels of focussing: in the case of the booster, even the weak focussing from the sector dipoles was too much in the horizontal plane and the D quadrupoles had to be powered slightly while keeping the F quadrupoles turned off!

A third issue is that $(C,V)$ contain an odd number of variables (15).  A stable mapping is similar to a rotation and rotations in an odd number of dimensions have an invariant axis.  This does not happen for the six-variable optimisation to find $\mathbf s$, since rotations in an even number of dimensions must only leave a point invariant.  The location on the invariant axis is a free parameter, meaning the matched ellipsoid solution for $(C,V)$ is non-unique.  In the case of the booster, this parameter was observed to correspond to the angular momentum of the bunch in the vertical axis and was set to zero with an additional constraint $V_{xz}=V_{zx}$ (it was not clear a bunch could be easily spun during injection).  Somewhat intuitively, the matched ellipsoids with high angular momentum had a more oblate shape.

As well as rings, variations on this matching method could match desired ellipsoid parameters at the start and ends of transfer lines, if enough magnet parameters are optimisable.

% Table of `nominal' booster magnet settings, beam energy, species, machine circumference, rev freq, bare tunes
\begin{table}[!hbt]
   \centering
   \caption{Ring Parameters and Magnet Settings}
   \begin{tabular}{lcc}
       \toprule
%\textbf{Parameter} & \textbf{Value} & \textbf{Unit} \\
%       \midrule
Ion species & Ca$^+$ & \\
Kinetic energy per mass & 3.2332 & MeV/u \\
Main dipole field $B_\mathrm{dip}$ & 0.75518 & T \\
F quadrupole gradient & 0 & T/m \\
D quadrupole gradient & 0.2255 & T/m \\
Machine circumference & 201.78 & m \\
Revolution frequency & 123.47 & kHz \\
Zero current tunes $Q_{0x,y}$ & 1.0100, 1.0979 & cycles \\
       \bottomrule
   \end{tabular}
   \label{tab:ringparams}
\end{table}

% Table of matched ellipsoid parameters, bunch population etc.
\begin{table}[!hbt]
   \centering
   \caption{Matched bunch ellipsoid parameters.  Note that $C_{ii}=\sigma_i^2$ and other omitted entries are zero.}
   \begin{tabular}{lcc}
       \toprule
%\textbf{Parameter} & \textbf{Value} & \textbf{Unit} \\
%       \midrule
Number of ions $N$ & 2560 &  \\
Bunch size $\sigma_{x,y,z}$ & 211.3, 52.9, 207.5 & $\mu$m \\
$C_{xz}$ & -130 & $\mu$m$^2$ \\
$V_{xx}$ & -115911 & s$^{-1}$ \\
$V_{yy}$ & 118810 & s$^{-1}$ \\
$V_{zz}$ & 36 & s$^{-1}$ \\
$V_{xz}=V_{zx}$ & -391197 & s$^{-1}$ \\
       \bottomrule
   \end{tabular}
   \label{tab:bunchparams}
\end{table}

The ring parameters used in this paper are given in Table~\ref{tab:ringparams} and the bunch shape at the start of the lattice, derived using the above matching method, is given in Table~\ref{tab:bunchparams}.  The number of ions is chosen semi-arbitrarily to allow pairwise space charge calculations but the dynamics scale with $\sigma\propto N^{1/3}$.

\section{Multiparticle Simulation}
Paul type ion traps have been built at several accelerator labs \cite{SPODions,IBEX,IBEXions} and the addition of laser Doppler cooling \cite{Metcalf} can produce Coulomb crystals with a solid structure and temperatures in the milli-Kelvin range \cite{SPODpic}.  A simulation of such a cooled trap, using the same code as \cite{BrooksCOOL25}, produced an input distribution with temperature 5.8\,mK and aspect ratio as required by the matched ellipsoid.  A linear transformation was applied to give the matched $C$ and $V$ values.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{WEP4324_f3}
   \caption{Input bunch distribution (left) and bunch after 45\,ms in the booster (right).}
   \label{fig:trackbunch}
\end{figure}

Implementing the booster field model on the GPU instead of the Paul trap allowed the booster ring to be tracked.  $\delta t$ was set the same as for the ellipsoid model.  Figure~\ref{fig:trackbunch} shows the bunch before and after a 45\,ms (about 5500 turns) simulation, with only a small amount of emittance growth visible.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{WEP4324_f4}
   \caption{Temperature evolution during tracking.}
   \label{fig:tracktemperature}
\end{figure}

Figure~\ref{fig:tracktemperature} shows how temperature changes during the simulation, with an average heating rate of 0.5291\,K/s.  The temperature is defined by
$ T_n = \frac{m}{3k_B} \langle|\mathbf v - \mathbf p_n(\mathbf x)|^2\rangle $
where $\mathbf p_n(\mathbf x)$ is the $n^\mathrm{th}$ order polynomial least squares fit to $\mathbf v(\mathbf x)$.  Thus, $T_1$ subtracts both bunch average motion and linear lattice motion before calculating the temperature, whereas $T_2$ also removes a quadratic component, which makes a negligible difference here.  There is an initial jump to around 12\,mK within the first turn, which may be due to the equilibrium lattice structure being different between the Paul trap and the booster, followed by a heating rate that increases with temperature, which is expected when in the solid phase \cite{OkamotoIBS}.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{WEP4324_f5}
   \caption{Temperature variation with quadrupole strength.}
   \label{fig:dquadstrength}
\end{figure}

Figure~\ref{fig:dquadstrength} shows the result of varying the D quadrupole strength on the temperatures $T_{1,2}$ in multiparticle simulations at $t$=1\,ms.  Below $\sim$0.19\,T/m, there are no matched solutions of the model and above $\sim$0.27\,T/m the solution has unstable eigenvalues.  The temperature rise shown is due to higher order effects, likely crystal shear heating from the bunch aspect ratio changing, with minimal heating occuring close to where $\sigma_z/\sigma_x\simeq1$ and exponential growth either side.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{WEP4324_f6}
   \caption{Covariance matrix evolution through one superperiod ($\frac16$ turn) expressed in the local frame of the bunch.  Entries not plotted are approximately zero.  Uniform ellipsoid model predictions are shown with dashed lines.}
   \label{fig:covsp}
\end{figure}

The evolution of the bunch shape through one superperiod of the booster is shown in Fig.~\ref{fig:covsp}.  Here, the lattice structure of eight cells with missing dipoles in the third and sixth cells is visible.  Bunch sizes $\sigma_{x,y,z}$ only change slightly, with the most noticeable dynamics in the shear $C_{xz}$, which accumulates negatively in the drifts and reverses direction in the dipoles to give an overall periodic function.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{WEP4324_f7}
   \caption{Velocity matrix evolution in one superperiod.}
   \label{fig:velsp}
\end{figure}

The nonzero $V$ matrix elements are plotted in Fig.~\ref{fig:velsp}.  If $C$ is roughly analogous to the beta functions, $V$ corresponds to alpha.  The $V_{xx}$ and $V_{yy}$ transverse elements change rapidly in opposite directions in quadrupoles, with dipole weak focussing also visible making $V_{xx}$ decrease.  The zero vertical angular momentum condition $V_{xz}=V_{zx}$ stops being true in dipoles, since canonical rather than kinetic angular momentum is conserved and they differ in a magnetic field.

The bunch rotates $\sim$0.495 times in the $x$--$z$ plane on each turn, in the same direction as the ring revolution.

\section{Conclusion}
A regime has been found where milli-Kelvin bunch temperatures can persist in the BNL Booster ring.  The timescales approach that of a full machine acceleration cycle, so future work will investigate acceleration.  Stable parameters for the ring have been found up to 117.1\,MeV/u in the model and await confirmation with particle tracking.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{WEP4324_f8}
   \caption{Solitons forming from a non-crystalline beam.}
   \label{fig:solitons}
\end{figure}

The persistence of the bunch without bunching RF is surprising and likely relates to a Coriolis-like force from the ring rotation, which turns repulsive space charge into a rotational force.  The form of the matched bunch, flattened vertically with nearly circular shape in the horizontal-longitudinal plane, is strongly reminiscent of the `vortex effect' observed in cyclotrons near the space charge limit \cite{BaartmanSolitons}.  As Coulomb crystals are at the space charge limit themselves, this may be the same phenomenon.  To test this theory, a longer non-crystalline bunch of temperature 10\,K was placed in the booster lattice with the same settings, resulting in Fig.~\ref{fig:solitons}.  The bunch breaks up into persistent circular bunches resembling the last figure in \cite{BaartmanSolitons}.


\begin{thebibliography}{99}   	% Use for  10-99  references
%\begin{thebibliography}{9} 		% Use for 1-9 references
	
\bibitem{CourantBoosterLattice} E. Courant and Z. Parsa, ``The Booster lattice'', Booster Technical Note No. 1/BNL--105049-2014-IR/BNL-105049-2014-TECH (1986).  \doi{10.2172/1150383}
		
\bibitem{BrooksOffAxis} S.J. Brooks, ``Off-Axis Magnetic Fields Extrapolated from On-Axis Multipoles'', BNL-223622-2022-TECH (2013).  \doi{10.2172/1895080}

\bibitem{LinearA} Wei Cai, ``Potential Field of a Uniformly Charged Ellipsoid'', Department of Mechanical Engineering, Stanford University (May 28, 2007).  Available from \url{http://micro.stanford.edu/~caiwei/me340a/A_Ellipsoid_Potential.pdf}

\bibitem{LinearB} O.D. Kellogg, ``Foundations of Potential Theory'', Dover, New York (1953).

\bibitem{LinearC} N.M. Ferrers, ``An Elementary Treatise on Spherical Harmonics and Subjects Connected with Them'', Macmillan and Company, London (1877).  \\ Available from \url{https://archive.org/details/elementarytreati00ferruoft/page/108/mode/2up}

\bibitem{firstordermodel} S.J. Brooks, ``First Order Transport of a Cold Uniform Ellipsoid of Charge'', note BNL-229092-2025-TECH (2025).  \doi{10.2172/3000537}

\bibitem{BrooksCOOL25} S.J. Brooks, ``Proposed Ultralow-Emittance Beam Source for High-Luminosity Hadron Colliders'', in \emph{Proc. COOL'25}, Stony Brook, NY (Oct.~2025).  Available from \url{https://stephenbrooks.org/ap/report/2025-10/WEB3.pdf}  \doi{10.18429/JACoW-COOL2025-WEB3}
		
\bibitem{BrooksSVD} S.J. Brooks, ``Bounded Approximate Solutions of Linear Systems using SVD'', BNL-223624-2022-TECH (2015).  \doi{10.2172/1895088}
\bibitem{BrooksOptimiser} S.J. Brooks, ``Higher-Order Corrections to Optimisers based on Newton’s Method'', arXiv:2307.03820 (2023).


\bibitem{SPODions} R. Takai, H. Enokizono, K. Ito, Y. Mizuno, K. Okabe and H. Okamoto, ``Development of a Compact Plasma Trap for Experimental Beam Physics'', \emph{Japan. J. Appl. Phys.} \textbf{45}, No. 6A pp.5332--5343 (2006).

\bibitem{IBEX} S. L. Sheehy, E. J. Carr, L. K. Martin, K. Budzik, D. J. Kelliher, S. Machida and C. R. Prior, ``Commissioning and First Results of the IBEX Linear Paul Trap'', in \emph{Proc. IPAC 2017}, Copenhagen, Denmark, May~2017, pp.4481--4484.\\ \url{doi:10.18429/JACoW-IPAC2017-THPVA027}
	
\bibitem{IBEXions} L. K. Martin, S. Machida, D. J. Kelliher and S. L. Sheehy, ``A study of coherent and incoherent resonances in high intensity beams using a linear Paul trap'', New J. Phys. \textbf{21} 053023 (2019).

\bibitem{Metcalf} H. J. Metcalf and P. van der Straten, ``Laser Cooling and Trapping of Neutral Atoms'', in The Optics Encyclopedia, Vol. 38, pp.847--853 (Wiley-VCH Verlag GmbH \& Co. KGaA, Weinheim, Germany, 2007).

\bibitem{SPODpic} K. Izawa, K. Ito, H. Higaki and H. Okamoto, ``Controlled Extraction of Ultracold Ions from a Linear Paul Trap for Nanobeam Production'', \emph{J. Phys. Soc. Japan} Vol. \textbf{79}, No. 12, p.124502 (2010).


\bibitem{OkamotoIBS} H. Okamoto, H. Sugimoto and Y. Yuri, ``Coulomb Coupling and Heating of Charged Particle Beams in the Presence of Dispersion'', J. Plasma Fusion Res. SERIES, Vol. \textbf{8}, pp.950--954 (2009).


\bibitem{BaartmanSolitons} R. Baartman, ``Space Charge Limit in Separated Turn Cyclotrons'', in \emph{Proc. Cyclotrons'13}, Vancouver, Canada, pp. 305--309 (2013).

\end{thebibliography}

\end{document}