# Multi-turn injection into a heavy-ion synchrotron in the presence of space
charge

S. Appel  O. Boine-Frankenheim  GSI, Helmholtzzentrum für Schwerionenforschung
GmbH, Planckstr. 1, 64291 Darmstadt, Germany  Technische Universität
Darmstadt, Schlossgartenstraße 8, 64289 Darmstadt, Germany

###### Abstract

For heavy-ion synchrotrons an efficient Multi-Turn Injection (MTI) from the
injector linac is crucial in order to reach the specified currents using the
available machine acceptance. The beam loss during the MTI must not exceed the
limits determined by machine protection and by the vacuum requirements.
Especially for low energy and intermediate charge state ions, the beam loss at
the injection septum can cause a degradation of the vacuum and a corresponding
reduction of the beam lifetime. In order to optimize the injection of intense
beams a very detailed simulation model was developed. Besides the closed orbit
bump, lattice errors, the position of the septum and other aperture limiting
components the transverse space charge force is included self-consistently.
The space charge force causes a characteristic shift of the optimum tunes and
a smoothing of the phase space density.

###### keywords:

Synchrotron , multi-turn injection , space charge , simulation

††journal: arXiv

##  1 Introduction

In heavy-ion synchrotrons the Multi-Turn-Injection (MTI) should fill up the
available transverse acceptance, avoiding the already occupied phase space
area [1, 2, 3]. In order to increase the space charge limit, heavy-ion
synchrotrons are operated with intermediate charge state ions. Therefore
stripping injection is not an option and the MTI has to respect Liouville’s
theorem for the chosen charge state. The present study focuses on the MTI into
the GSI heavy-ion synchrotron SIS18 [4]. The SIS18 is presently being upgraded
in order to increase the beam intensities for the FAIR (Facility for
Antiproton and Ion Research) project [5]. For the reference U28+ ions the
SIS18 intensities at injection energy (11.4 MeV/u) should reach the limit
determined by transverse space charge and error resonances (space charge
limit). Heavy-ion ion beams are injected from the GSI UNILAC and are stacked
along the horizontal plane using multi-turn injection [4]. The beam loss
should be well below 30% to avoid an increase of the dynamic vacuum pressure
as well as activation and damage. The transverse beam size of the stacked beam
should be within the machine acceptance ($\epsilon_{x}=150$ mm mrad,
$\epsilon_{y}=50$ mm mrad). Due to the elliptic beam pipe the vertical machine
acceptance is much smaller than the horizontal. Also the resulting rms
momentum spread ($\Delta p/p\leq 10^{-3}$) has to remain below a maximum value
determined by the available rf bucket area for fast ramping. In the SIS18 the
microbunches injected from the UNILAC with 36 MHz debunch, filament and form a
coasting beam within a few turns. During debunching the longitudinal space
charge force of the intense microbunches leads to an increases of the final
momentum spread [6]. In this study we will first discuss optimum multi-turn
injection schemes for low intensities. Afterwards the impact of transverse
space charge on the optimum MTI settings and the resulting phase space
distribution is discussed.

##  2 Multi-turn injection model

Figure 1: The septum and bumper magnets

In the SIS18 the injected beamlets are stacked in the horizontal phase space
until the machine acceptance is reached. To fulfill Liouville’s theorem, four
bumper magnets create a time variable closed orbit bump such that the
electrostatic injection septum deflects the next incoming beamlet into free
horizontal phase space and close to the formerly injected beamlets. Fig. 1
indicates the position of the four bumper magnets along the circumference. If
we require that the four bumper magnets produce no closed orbit distortion
outside of the injection region and at the injection position $s_{0}$, then
the amplitude $x_{0}$ as well as the angle $x_{0}^{\prime}$ of the closed
orbit are the degrees of freedom. The resulting angular kicks produced by the
bumper magnets are [1, 2]

$\displaystyle\varphi_{1}$
$\displaystyle=\frac{d_{02}x_{0}-b_{02}x^{\prime}_{0}}{b_{21}}\ ,$
$\displaystyle\varphi_{2}$
$\displaystyle=\frac{-d_{01}x_{0}+b_{01}x^{\prime}_{0}}{b_{21}},$ (1)
$\displaystyle\varphi_{3}$
$\displaystyle=\frac{-b_{41}\varphi_{1}-b_{42}\varphi_{2}}{b_{43}}\ ,$
$\displaystyle\varphi_{4}$
$\displaystyle=\frac{b_{31}\varphi_{1}+b_{32}\varphi_{2}}{b_{43}}$ (2)

with

$\displaystyle b_{lk}$
$\displaystyle=\sqrt{\beta_{l}\beta_{k}}\sin(\phi_{l}-\phi_{k}),$ (3)
$\displaystyle d_{lk}$
$\displaystyle=\sqrt{\frac{\beta_{k}}{\beta_{l}}}\sin(\phi_{l}-\phi_{k})-\frac{\alpha_{l}}{\beta_{l}}b_{lk}.$
(4)

Here $\alpha_{l}$ and $\beta_{l}$ are the horizontal lattice parameters and
$\phi_{l}$ is the phase advance at the point $s_{l}$.
To achieve high beam intensities the injected beamlets should be packed as
compact as possible. In normalized phase space coordinates the injected
beamlets as well as the beam pipe are approximately circular, therefore the
MTI packing problem is similar to the packing of ropes and cables. Assuming
that $n_{eff}$ beamlets with radius $a$ are packed into a given machine
acceptance (or container) with radius $R$ like shown in Fig. 2 (also named
hexagonal packing) and imagine hexagons which circumscribe those beamlets.
Then the dilution is defined as the area of the container divided by the area
of these hexagons [7]

$d=\frac{2\sqrt{3}}{\pi}\frac{R^{2}}{a^{2}}\frac{1}{n_{eff}}.$ (5)

Figure 2: Optimal packing of seven identical beamlets of radius $a$ into a
machine acceptance with radius $R$.

The hexagonal packing of circles onto an annular ring leads to a shell
structure and the dilution will be larger than one. In [7] it is shown that
for $n_{eff}=7$ and $n_{eff}=19$ the smallest dilution is achieved, namely
$d_{7}=1.16$ and $d_{19}=1.12$.
The incoming beamlets will have a linear $x_{i}$ and an angular
$x^{\prime}_{i}$ displacement with respect to the closed orbit and will
therefore undergo betatron oscillations. One turn later the injected beamlets
will pass the injection point again, but due to the betatron oscillation will
avoid the septum. If the orbit is not reduced fast enough, the beamlets will
hit the inner side of the septum after $n$ revolution turns, depending on the
betatron oscillation tune and get lost. Therefore the beam loss has its
maximum located at the resonance condition

$Q_{x}n=integer$ (6)

where $Q_{x}$ is the horizontal tune. During injection the number of injected
turns $n_{i}$ and the revolution turns $n$ are equal, after injection this
condition does not hold any longer because $n_{i}$ remains constant. If $\eta$
characterizes the ratio between the lost and injected particles it follows

$n_{eff}=n_{i}(1-\eta).$ (7)

For a loss free injection $\eta$ is zero and the effectively accumulated
beamlets are equal to the number of injected turns. The quality of the MTI is
determined by the number of accumulated beamlets and the associated loss. For
a given machine acceptance a high MTI quality implies a small dilution factor.
The beamlets injected at latter turns will have a larger betatron amplitude as
the orbit bump is reduced. This leads to the formation of

$m\sim\begin{cases}Q_{f}^{-1},&\text{if }0<Q_{f}\leq 0.5\\\
(1-Q_{f})^{-1},&\text{if }0.5<Q_{f}\leq 1\end{cases}$ (8)

arms determined by the fractional tune $Q_{f}$. In our later considerations we
will assume fractional tunes between $0.0$ and $0.5$.
To achive small dilution factors one would choose a fractional tune, which
leads to a large number of arms. However, one has to consider the restriction
given by the geometry and by the machine resonances. If the first beamlet is
injected exactly on the closed orbit; one can place at maximum six beamlets
with the same size around this beamlet as shown in Fig. 2. To avoid beam loss
due to the already occupied phase space area, the fractional tune should not
be smaller than $1/6$. But then one has to accept a larger dilution factor for
the next shell than given by the geometry optimum (six instead of twelve).
This could be improved by ramping the tune using pulsed quadrupoles. However,
this option will not be considered in the present study.

Figure 3: Snapshot of the MTI from a simulation showing the shell structure in
normalized phase space.

Fig. 3 shows a snapshot of a MTI simulation in normalized horizontal phase
space for a fractional tune of $Q_{f}=0.17$ resulting from the present
arguments as optimal tune for the injection. The linear reduction of the bump
per turn

$\Delta x\approx Q_{f}(2a+d_{s})$ (9)

in the simulation has been chosen such, that the incoming beamlets miss the
inner side of the septum with their edges. The septum thickness $d_{s}$ can be
neglected for the SIS18 injection since the beamlet radius is in the range of
$5$ and $10$ mm and the septum thickniss is $0.1$ mm. The vertical solid line
in the figure indicates the septum. The newly injected beamlet is located to
the right of the septum and will then oscillate towards the stacked beamlets
to the left of the septum. The formation of six arms is clearly visible as
well as the resulting shell structure. In the simulation the decoherence due
to the chromaticity $\xi$ has been supressed to demonstrate more clearly the
influence of the coherent betatron oscillations of the injected beamlets. The
decoherence would fill the voids between the injected beamlets after many
turns.
The MTI quality can be further optimised by adjusting the start value of the
closed orbit bump, its ramping curve and the position of the incoming beamlet.
For an optimal filling of the phase space (e.g. small dilution factor) Ref.
[3] propose to map the injected emittances onto upright ellipses

$\frac{\alpha_{0}}{\beta_{0}}=-\frac{x^{\prime}_{i}-x^{\prime}_{0}}{x_{i}-x_{0}}$
(10)

and mismatch the lattice function

$\frac{\beta_{0}}{\beta_{i}}=\left(\frac{\epsilon}{\epsilon_{i}}\right)^{1/3}$
(11)

in order to adapt the curvature of the incoming beamlet to the ring acceptance
curvature. The subscript $i$ represents the parameter for the incoming beamlet
and the subscript $0$ for the ring. $\epsilon$ is the available ring
acceptance and $\epsilon_{i}$ is the emittance of the injected beamlet. The
beam center of the incoming beamlet should approximately placed as such, that
the edge of the incoming beamlet touches the septum at the outside (see Fig.
3)

$x_{i}\approx a+x_{s}+d_{s}+x_{st}$ (12)

with $x_{s}$ as the septum position and $x_{st}$ as the possible steering
offset. The angular center position of the incoming beamlet is limited by the
septum deflection angle $x^{\prime}_{s}$ and the possible steering of the
linac beam $x^{\prime}_{st}$

$x^{\prime}_{i}\approx x^{\prime}_{s}+x^{\prime}_{st}.$ (13)

For the orbit bump reduction also other ramps are proposed besides a linear
ramp to better adapt the orbit bump to the free phase space regions. An
exponential ramp is given through

$x(n)=x_{0}\exp(-(2a+d_{s})\tau n)$ (14)

for the orbit bump variation per turn. By varing $\tau$ one can possibly find
an efficient ramping curve which provides an optimal filling of the phase
space accompanied by low beam loss. To further improve the MTI quality also
the steering parameters $x_{st}$ and $x^{\prime}_{st}$ should be varied as
long the mapping condition into upright ellipses is fulfilled.

##  3 MTI simulation model

The MTI model described in the previous section has been implemented in the
particle tracking codes PATRIC (see e.g. [8]) and pyORBIT (see e.g [9, 10]).
Both codes use space charge solvers together with a Particle-In-Cell (PIC)
numerical scheme. The simulation model includes the closed orbit bump, lattice
errors and the position of the septum. To consider loss on machine acceptance
we placed an acceptance collimator vis-$\grave{a}$-vis to the septum. The
simulation model also contains the possibility to obtain the closed orbit bump
from the SIS18 control program approximation. For technical reasons the
present SIS18 control program uses an approximation for the calculation of the
angular kick from the bumper magnets, since the evaluation of the horizontal
lattice parameters as well as the phase advance at the bumper magnet position
and at the injection point of a given horizontal tune are not possible. To
find a function depence on horizontal tune for a fixed amplitude ($x_{0}=10$
mm) and angle ($x_{0}^{\prime}=1$ mrad) the angular kicks for several
horizontal tunes given by Eq. 1 were calculated and quadratic functions were
fitted on these results. By normalizing this function with the bump amplitude
of $x_{0}=10$ mm mrad one can adjust the four functions on each bump amplitude
by multiplying them with the desired bump amplitude [11]. The new control
system does not have such limitations and therefore this bump calculation will
be used only for the comparison with current perform SIS18 experiments.
In PATRIC the Poisson equation is solved on a static 2D transverse grid
together with an FFT solver and momentum kicks corresponding to the local
space charge field strength are applied. For the macro-particle tracking
PATRIC uses the linear transfer matrix computed by the lattice program MADX.
pyORBIT is a Python/C++ implementation of the ORBIT (Objective Ring Beam
Injection and Tracking) code. The code is freely available at [9]. In pyORBIT
the transverse space charge potential for the particle distribution is
evaluated on a dynamic 2D grid using an FFT solver and the kicks on the macro-
particles are obtained from second order interpolation [10]. In pyORBIT the
particles are tracked with a TEAPOT (Thin-Element Accelerator Program for
Optics and Tracking [12]) method. pyORBIT allows to build the TEAPOT lattice
from a MAD input file by analyzing the input file. We implemented a python
routine to include the model for the closed orbit bump.

##  4 Comparison of MTI simulations and experiment

Figure 4: Simulated and measured MTI loss as a function of the horizontal
tune. The loss maxima given by Eq. 6 are indicated as vertical lines. The
measurement results were provided by [13].

Accurate predictions of the MTI efficiency require a careful validation of the
machine model with experiments. The MTI efficiency depends on various machine
and beam parameters. Often, some important MTI parameters are not precisely
known or can change from cycle to cycle. Therefore, an accurate measurement of
the horizontal tune is most importent. In the SIS18 the machine tune is
obtained with very high precision from the transverse Schottky spectrum.
During a dedicated MTI experiment with $U^{28+}$ ions [13], besides the tunes
also the injected emittance was measured ($\epsilon_{i}=6.5\pm 0.7$ mm mrad)
in the transfer channel to the SIS18 a few meters before the injection point.
The comparison between measurements and simulations is shown in Fig. 4 for low
currents. The measurements as well as the simulations with both codes show
beam loss maxima located at the resonant tunes. In the simulations the machine
injection settings have been used (orbit amplitude of $x_{0}$ = 85 mm, bump
reducing of $\Delta x_{0}$ = 2.5 mm per turn and the number of injected turns
of $n_{i}=21$). In the SIS18 it is not possible to measure the slope and the
position of the incoming beamlet at the septum. The shown agreement between
measurements and simulations was obtained by setting the center of the
injected beamlets to $x_{i}$ = 90 mm (approximately the position of septum
plus twice the beamlet radius) and its slope to $x^{\prime}_{i}$ = 7.9 mrad
(the range of the linac steering angle is 1-3 mrad) in the simulations.
The experiment provided by [13] indicate that the fraction tune of
$Q_{f}=0.17$ is not always optimal for an efficiently injection into the
SIS18. This discrepency can be explained by a mismatch of the beamlet slope in
the experiment. In the simulation an even lower beam loss around $Q_{f}=0.17$
could be achieved by changing the beamlet slope. As shown in Fig. 5 for a
beamlet slope of $x^{\prime}_{0}=7.0$ mrad for tunes around 4.17 the loss is
reduced by $30\%$. Also a shift of the maxima is observable. In other
dedicated MTI experiments a lower loss for tunes around 4.17 could be obained
[13].

Figure 5: MTI loss as a function of the horizontal tune for two different
beamlet slopes.

##  5 Improvement of the SIS18 MTI quality

In order to find an efficient ramping curve we performed several simulations
by varying the parameter $\tau$ for different initial emittances (show in Fig.
6). The chosen emittance values are in the range of the expected emittances
for the planned UNILAC injector upgrade. In the simulation the injection is
stopped if the rms emittance of the stacked beam exceeds the machine
acceptance. This automatically maximizes the number of effectively accumulated
beamlets for a given initial emittance. The chosen injection stop condition
results in a different number of injected turns $n_{i}$ in each simulation. To
take into account any possible subsequent beam loss the simulation continues
for 200 turns. As indicated in Fig. 6 $\tau=3$ gives a good compromise between
a compact packing of the beamlets and loss. For all emittances a MTI qualtiy
with $\eta\approx 0.15$ loss and dilution factor $d\approx 2$ could be
achieved. If higher loss during injection is acceptable a slower reduction
(e.g. a smaller $\tau$) could be selected to reach a larger number of
efficiently accumulated beamlets. The speed of the orbit reduction is
reflected in the dilution factor. For a slower reduction the dilution could be
reduced to $d=1.5$ at the cost of a larger loss. The simulation results
indicate that for smaller emittances the dilution is slightly larger then for
larger emittances. The achieved dilution factors are similar to the dilution
factors of $d=1.4-2$ presented by [14] for a conventional one plane multi-turn
injection.
For the SIS18 an optimum steering offset of the incoming beamlet was found for
$x_{st}=-0.2a$ and $x^{\prime}_{st}=0.5$ mrad as well as
$x^{\prime}_{i}-x^{\prime}_{0}=0.25$ mrad and $x_{i}-x_{0}=2.8$ mm. The ratio
of the lattice parameters at the injection point in the SIS18 for the
horizontal tune of $Q_{f}=4.17$ is $\frac{\alpha_{0}}{\beta_{0}}=-0.09$ and
the SIS18 septum deflection angle is $x^{\prime}_{s}=6$ mrad.

Figure 6: The number of effectively stored turns and the beam loss as a
function of the ramping rate $\tau$ for four different emittances.

##  6 Space charge effects

The PATRIC and pyORBIT codes were used to study the effects of the transverse
space charge force on the multi-turn injection. The incoherent tune shift for
a round, incoming beamlet with a homogeneous density is [15]

$\Delta Q_{sc}=-\frac{r_{p}Z^{2}N_{I}}{2\pi\beta^{2}\gamma^{3}A\epsilon_{i}}.$
(15)

In case of the stacked beam, the emittance $\epsilon_{i}$ in the upper
equation should be replaced by

$\frac{1}{2}\left(\epsilon_{x}+\sqrt{\epsilon_{x}\epsilon_{y}\tfrac{Q_{x}}{Q_{y}}}\right)$
(16)

and multiplied with the form factor $g_{f}$ to consider the non-homogeneous
density profile of the stacked beam. For a homogeneous profile $g_{f}$ is $1$
and for a Gaussian distribution $g_{f}$ is $2$. $\epsilon$ is four times the
rms emittance of the stacked beam. Suppose the vertical emittance remains
constant ($\epsilon_{y}=\epsilon_{i}$) than the horizontal space charge tune
shift for the stacked beam is

$\Delta Q_{sc}^{s}=\Delta
Q_{sc}\frac{2g_{f}}{d+\sqrt{\tfrac{d}{n_{eff}}\tfrac{Q_{x}}{Q_{y}}}}.$ (17)

During multi-turn-injection, the increase of the intensity and of the
emittance usually compensate each other and the horizontal space charge tune
shift remains below the initial one

$\Delta Q_{sc}^{s}\lessapprox\Delta Q_{sc}.$ (18)

For the SIS18 this condition is fulfilled in all our simulation examples.
Other injection schemes, e.g. employing two planes with a tilted septum, could
possibly lead to larger space charge parameters for the the stacked beam than
for the incoming beamlet.

Figure 7: Simulation result obtained with PyORBIT and PATRIC for the MTI
efficiency as a function of the horizontal tune with and without space charge
effects.

Fig. 7 shows the effect of space charge on the MTI efficiency. One can observe
that the maxima and minima of the MTI efficiency are shifted to the right by
$\Delta Q_{sc}$ under the influence of space charge. The observed shift of the
maxima/minima can be described very well through

$(Q_{x}+\left|\Delta Q_{sc}\right|)n=integer.$ (19)

In the simulations the space charge tune shift for an injected beamlet was
$\Delta Q_{sc}\approx-0.04$.
The simulations indicate that the presently chosen high working point
$Q_{x}=4.17$ in the SIS18 is a good candidate in order to reach the FAIR
design intensities. We do expect that the MTI efficiency for high beam
currents (with space charge) will be close to the one achieved for low
currents (low space charge) shown in Fig. 6.

Figure 8: The snapshot of an MTI simulation with space charge (same simulation
parameters as in Fig. 3)

During and right after the injection process the stacked beam distribution
will be non-stationary. The space charge potential energy will be larger than
the potential energy of the corresponding matched equilibrium distribution.
Therefore we can expect that the nonlinear space charge force will drive the
beam distribution closer towards an equilibrium [16]. In addition also
chromaticity or beam mismatch causes a filamentation of the distribution
between the injected beamlets. Fig. 8 shows a snapshot of the broadening of
the beamlets and a filamentation driven by space charge during the multi-turn
injection. Besides the much higher intensity, the simulation parameters are
similar to the ones chosen for Fig. 3. Especially close to the center, with
space charge, the individual beamlets cannot be distinguished. At a latter
stage a phase space distribution with a few larger holes can develop. In the
center of Fig. 8 the initial stage of such a hole can already be observed. For
the chosen beam intensities the evolution of the beam distribution is
dominated by the effect of space charge. Other effects, like the decoherence
due to chromaticity, can be neglected.

##  7 The role of the vertical emittance

Figure 9: Increase of the space charge tune shifts during injection if the
vertical emittance remains constant. The vertical line indicates the
horizontal space charge tune shift of the incoming beamlet.

In our simulation model the vertical emittance remains constant during and
also after the horizontal MTI. The corresponding evolution of the vertical and
horizontal space charge tune shifts are shown in Fig. 9 for the example case
of an initial vertical emittance of $\epsilon_{y}=15$ mm mrad. During the
injection the beam loss as well as the dilution varies turn by turn. Presently
the vertical emittance in the SIS18 increases after injection by
$\frac{\Delta\epsilon_{y}}{\epsilon_{y}}\approx 0.2$. In [17] the vertical
emittance was obtained from the vertical beam profiles (for an
${}^{40}Ar^{18+}$ beam). The initial UNILAC emittance was measured by [13] as
$\approx 15$ mm mrad. Analytical estimates indicate that the emittance growth
in the vertical plane is caused by injection steering errors and optical
mismatch, as we will show in this section. The growth of the vertical
emittance affects the horizontal space charge tune shift. More importantly the
vertical emittance determines the (vertical) space charge limit at $\Delta
Q_{y}\approx-0.5$ in the SIS18, which is beyond the scope of the present
study. In an optimal scenario the vertical emittance should increase in a
controlled way up to the vertical machine acceptance, which is at $A_{y}=50$
mm mrad in the SIS18 [5].
During the horizontal MTI the vertical emittance increases in a non-controlled
way due to steering errors and optical mismatch [18, 19]. The emittance
increases due to a radial displacement of the beam center after an injection
steering error

$\frac{\Delta\epsilon_{y}}{\epsilon_{y}}=\frac{\Delta
r^{2}}{2\epsilon_{y}\beta_{y}},$ (20)

if the radial displacement is

$\Delta r=\sqrt{\Delta y^{2}+(\alpha_{y}\Delta y+\beta_{y}\Delta
y^{\prime})^{2}}.$ (21)

Assuming the steering errors are $\Delta y=0.5$ mm and $\Delta y^{\prime}=0.5$
mrad the emittance growth is $\frac{\Delta\epsilon_{y}}{\epsilon_{y}}=0.1$.
The emittance growth is very sensitive to the slope error due to the
multiplication with the beta function. A doubling of the slope error leads to
$\frac{\Delta\epsilon_{y}}{\epsilon_{y}}=0.4$. For a mismatch of the lattice
functions the emittance growth is

$\frac{\Delta\epsilon_{y}}{\epsilon_{y}}=F-1$ (22)

with

$F=\frac{1}{2}\left[\frac{\beta_{0}}{\beta_{i}}+\frac{\beta_{i}}{\beta_{0}}+\left(\alpha_{0}-\alpha_{i}\frac{\beta_{0}}{\beta_{i}}\right)^{2}\frac{\beta_{i}}{\beta_{0}}\right].$
(23)

Assuming the lattice functions are known with an error of 20$\%$ we obtain an
emittance growth of $\frac{\Delta\epsilon_{y}}{\epsilon_{y}}=0.2$.
A mismatch in the dispersion gives the following change in emittance

$\frac{\Delta\epsilon_{y}}{\epsilon_{y}}=\frac{\Delta
D^{2}}{2\epsilon_{y}\beta_{y}}\sigma_{p}^{2}$ (24)

if the radial dispersion displacement is

$\Delta D=\sqrt{\Delta D_{y}^{2}+(\alpha_{y}\Delta D_{y}+\beta_{y}\Delta
D_{y}^{\prime})^{2}}.$ (25)

Since the injection kicker bends the incoming beamlet vertically there is a
systematic vertical dispersion mismatch which results in
$\frac{\Delta\epsilon_{y}}{\epsilon_{y}}=5\times 10^{-5}$.
The total vertical emittance increase after injection due to the
aforementioned effects should not be larger than
$\frac{\Delta\epsilon_{y}}{\epsilon_{y}}=0.3$, which is larger than the
observed $\frac{\Delta\epsilon_{y}}{\epsilon_{y}}\approx 0.2$.

##  8 Conclusions and Outlook

We studied the effect of transverse space charge on the multi-turn injection
(MTI) using the example of the SIS18 at GSI. Two different simulation codes,
PATRIC and pyORBIT, where adapted to the problem. The MTI quality (intensity
and beam loss) depends very strongly on the horizontal tune. The beam loss
maxima are located at the lower order resonances. An optimal fractional tune
is identified at 1/6. The transverse space charge force causes a shift of the
beam loss maxima and of the optimum tune. We find that this shift corresponds
to the space charge tune shift of the incoming beamlets. For the SIS18
conditions we find an exponential decrease of the orbit bump results in good
compromise between compact stacking and beam loss. More or less independently
from the injected emittance the beam loss could be reduced to $\eta=0.15$ due
to a rather compact packing (e.g. $d=2$). A direct comparison of the beam
distribution with and without space charge indicates that the nonlinear space
charge force between the injected beamlets drives the final beam distribution
towards a more homogeneous density with a lower total potential energy. Future
studies will focus on the application of the skew quadrupoles in the SIS18 to
possibly increase the number of injected turns through a coupling with the
vertical plane. Another option to increase the MTI efficiency is to ramp the
tune during injection.

## Acknowledgement

The authors thank Y. El-Hayek for his measurement results and D. Ondreka for
useful discussions about the injection into the SIS18.

## References

## References

  * [1] C. Gardner, Multiturn injection of heavy ions into the booster, Tech. rep., Brookhaven National Laboratory, Upton, New York (1991).
  * [2] K. Wille, The Physics of Particle Accelerators: An Introduction, Oxford University Press on Demand, 2000\.
URL http://books.google.de/books?id=k1STXLP4lkgC

  * [3] C. R. Prior, Multiturn Injection into Accumulators for Heavy Ion Inertial Fusion, in: 5th European Particle Accelerator Conference, 1996, p. 3.
  * [4] B. Franczak, SIS 18 Parameterliste, 10.09.1987 (1987).
  * [5] P. Spiller, R. Bär, U. Blell, O. Boine-Frankenheim, L. Bozyk, K. Blasche, J. Falenski, E. Fischer, E. Floch, G. Franchetti, B. Franczak, P. Forck, O. Gumenyuk, I. Hofmann, P. Hülsmann, M. Kauschke, M. Kirk, H. Klingbeil, H. G. König, A. Kovalenko, P. Kowina, A. Krämer, D. Krämer, M. Kumm, U. Laier, M. Mehler, J. P. Meier, G. Moritz, P. Moritz, C. Mühle, K.-P. Ningel, C. Omet, I. Pschorn, N. Pyka, H. Ramakers, P. Schnizer, G. Schreiber, C. Schroeder, M. Schwickert, Y. Shim, M. Sitko, B. Skoczen, S. Sorge, P. Spiller, J. Stadlmann, A. Stafiniak, K. Sugita, B. Weckenmann, H. Welker, FAIR Technical Design Report December 2008, Tech. Rep. December (2008).
  * [6] S. Appel, O. Boine-Frankenheim, Microbunch dynamics and multistream instability in a heavy-ion synchrotron, Phys. Rev. ST Accel. Beams 054201 (2012) 1–10.  doi:10.1103/PhysRevSTAB.15.054201.
  * [7] S. Kravitz, Packing cylinders into cylindrical containers, Mathematics magazine 40 (2) (1967) 65–71.
URL http://www.jstor.org/stable/10.2307/2688509

  * [8] O. Boine-Frankenheim, V. Kornilov, S. Paret, Measurement and simulation of transverse Schottky noise with space charge, Physical Review Special Topics Accelerators and Beams 11 (7) (2008) 1–7.  doi:10.1103/PhysRevSTAB.11.074202.
  * [9] S. Cousineau, A. Shishlo, J. Holmes, PyORBIT.
URL https://code.google.com/p/py-orbit/

  * [10] A. Shishlo, S. Cousineau, V. Danilov, J. Galambos, S. Henderson, J. Holmes, M. Plum, O. Ridge, THE ORBIT SIMULATION CODE : BENCHMARKING AND APPLICATIONS *, in: Proc of ICAP 2006, 2006, pp. 53–58.
  * [11] D. Ondreka, GSI, System planning, privat communication.
  * [12] L. Schachinger, R. Talman, TEAPOT : A THIN-ELEMENT ACCELERATOR PROGRAM FOR OPTICS AND TRACKING 22 (1987) 35–56.
  * [13] Y. E. Hayek, Minimierung der systematischen Anfangsverluste im SIS18, Ph.D. thesis (2013).
  * [14] A. Chao, Handbook of accelerator physics and engineering, 1999.
  * [15] K. Schindl, Space Charge, in: D. Brant (Ed.), Proc, CAS, Intermediate Course on Accelerator Physics, CERN, 2006.
  * [16] M. Reiser, Theory and Design of Charged Particle Beams, Wiley Series in Beam Physics and Accelerator Technology, Wiley-VCH, 2008.
  * [17] S. Paret, V. Kornilov, O. Boine-Frankenheim, T. Weiland, Transverse Schottky and beam transfer function measurements in space charge affected coasting ion beams, Phys. Rev. ST Accel. Beams 13 (2) (2010) 22802.  doi:10.1103/PhysRevSTAB.13.022802.
  * [18] D. Möhl, Sources of emittance growth, in: D. Brandt (Ed.), Proceedings of the CAS. Intermediate accelerator physics, CERN, Zeuthern, Germany, 2005, pp. 245–268.
  * [19] W. Fischer, W. Mackay, S. Peggs, J. Wei, Emittance growth in RHIC during injection, Tech. rep., Brookhaven National Laboratory, Upton New York (1996).


This section decribes methods used for the three main stages of the SpineNet
pipeline (VB detection, VB labelling and IVV grading) from an algorithmic
standpoint. These stages are illustrated in Figure
LABEL:fig:spinenet_overview. Broadly speaking, whole spine scans and smaller
scans (e.g. lumbar, cervical) are processed identically, with only minor
differences in how the patch splitting is done in the VB detection stage.
These differences are described fully in Section 0.1.1, however we will
initially assume a single lumbar scan when describing the processing pipeline.

##  0.1 VB Detection by Vector Field Regression

The first stage in the SpineNetV2 pipeline is to detect vertebral bodies (VBs)
in the raw scans. This is done using a method called Vector Field Regression
(VFR). Operating in each sagittal slice of the scan, a network detects
gaussian responses over the centroids and corners (vertebral landmarks) of
each visible VB. For each corner detected, a corresponding vector field is
also output. At the point of each corner, this corresponding vector field
should ‘point’ to the centroid of the vertebra to which it belongs. This
allows landmarks from the same vertebra to be grouped together in a simple
manner which is robust to rotations, flips and variable vertebra size.

At this point, each detection consists of a centroid and four corners grouped
together to make a quadrilateral in an individual sagittal slice. Slicewise 2D
quadrilaterals corresponding to the same vertebra must then be grouped
together across slices forming 3D volumes. This is done are by measuring the
intersection-over-union (IOU) of quadrilaterals in neighbouring slices. If the
IOU is large, the quadrilaterals are grouped together into a single vertebral
instance. A diagram of this process going from raw scans to 3D volumes is
shown in Figure 1.

Figure 1: An overview of the VFR (Vector Field Regression) method. A U-Net has
13 output channels; 5 detection channels which output gaussian heatmaps over
the four corners and centroid of each visible vertebral body and 8 grouping
channels; the $x$ and $y$ components of vector fields corresponding to each
corner.

Specifically, inference proceeds as follows:

  1. 1.

A raw MR scan of dimension $S\PLH H\PLH W$ is split into $S$ sagittal slices.
These slices are in turn split into square patches of size 25 cm2 with an
overlap of 30% between neighbouring patches on each side. This patch-splitting
allows the network to deal with scans of varying dimensions.

  2. 2.

The patches are resampled by cubic interpolation to resolution of $224\PLH
224$. The resamples patches are fed as input to a ResNet50-encoded U-Net. This
results in a $13\PLH 224\PLH 224$ output per patch; 5 detection channels and 8
grouping channels (see Figure 1).

  3. 3.

By resampling & concatenating these patch-level outputs using the median
output for overlapping regions, an slice-level output tensor is reconstructed
of size $13\PLH H\PLH W$.

  4. 4.

The detection channels are then decoded to find the location of detected VB
centroids and corners. This is done by thresholding each channel & finding all
connected component in the resulting binary map. The exact point-of-detection
for each landmark is the point of maximal response in the detection channel
for the corresponding connected component.

  5. 5.

Now that each VB landmark in the slice has been detected, they must be grouped
into quadrilaterals corresponding to individual VBs. This is done by measuring
the value of the corresponding vector field at the point-of-detection for each
corner landmark. Looping through each detected VB centroid, the corner
landmark of each type which points closest is assigned as the corresponding
corner for that centroid, forming a quadrilateral. If no arrow is within a
distance from the centroid of 50% of the arrow length, that centroid is
discarded as a spurious detection.

  6. 6.

Finally the resulting VB polygons are grouped across neighbouring slices if
they have an IOU over 0.25. If more than one polygon in a neighbouring slice
overlap, then the one with the greatest IOU is chosen.

###  0.1.1 Splitting Large Scans Into Patches

Larger non-square scans (such as whole spine scans) are split into patches
before VFR is applied. This is done by splitting the scan into a grid of
overlapping squares with edge length $50$cm (as determined by the pixel
spacing parameter in the DICOM header) and an overlap of 40% between
neighbouring patches. The output from the detection and grouping channels are
then used to find landmarks in each patch. These landmarks are then
transformed back into the frame of the original scan. At this point, the
algorithm proceeds as before by grrouping landmarks into slicewise polygons
and then across saggital slices. Figure 2 shows the process of patch-splitting
for a whole spine scan.

Figure 2: Splitting a large scan into patches before applying VFR.

##  0.2 Convolutional Labelling of Vertebral Levels

The VFR method described above allows us to detect vertebral bodies in a
sagittally sliced MR scan. The next step of the pipeline is to determine which
vertebral levels the detections correspond to (e.g. S1, L5, L4 etc). This is
made more challenging by the fact that SpineNetV2 is not constrained to a
single field of view. As such, ‘counting up’ methods, i.e. those that rely on
an anchor vertebra being visible (such S1/C2 at the bottom/top of the
vertebral column) are unsuitable. Furthermore, such methods are not robust to
missed detections or variations in the number of vertebra in the vertebral
column (e.g. in cases where a transitional vertebra is present).

There are two pieces of information to consider when labelling a vertebra -
its appearance (e.g. intensity pattern, shape, size etc.) and its context (the
vertebra’s position relative to other detections in the scan). For example, S1
usually has a very distinctive shape which allows it to be labelled from
appearance alone. On the other hand, L5 looks very similar to the other lumbar
vertebrae, however can be easily identified from it’s context - it is the next
vertebra up from S1. Our method attempts to use both these sources of
information when labelling a vertebra.

Firstly, a 3D volume around each detected VB is fed as input to an appearance
network. This outputs a 23-element (from S1 to C3) probability vector
predicting the level of the VB from appearance alone. To include information
about the spatial configuration of the detections, this is then used to
contruct a probability-height map, $P$. At the height of each detection, $P$
has value equal to the output from the appearance network for that detection.
Using this as input, a convolutional context network refines the height-
probability map, taking into account appearance predictions from spatially
nearby VBs to update the probability vectors for each detection. The result is
a refined height-probability map $P^{\prime}$. This process is shown in Figure
3. Finally, $P^{\prime}$ needs to be decoded into discrete level predictions.
Naïvely, this could be done by taking the maximum probability level at the
height of each detection in $P^{\prime}$. However, this would allow for
nonsensical outputs, such as the same level for two detections. Ideally, we
also want to build in soft contraints such that successive detections are
labelled as successive labels. For example, we would expect S1 to be the
detection below L5. However, we also wish to remain robust to missed
detection.

To build these constraints into our approach, we take inspiration from
language modelling. Using a beam search, we can find the most probable valid
sequence of levels for the detections. Penalties are added to the sequences
probability score in the case of transitional vertebrae or numerical
variations to reflect the unlikeliness, yet possibility, of such events. The
exact parameterisation of this search is specified in the next paragraph.

Figure 3: The language modelling-inspired approach to vertebra labelling
outlined in this section. We begin by constructing a probability-height map
$P$ based on appearance information alone. Using a context-aware CNN this map
is refined ($P^{\prime}$). Finally $P^{\prime}$ is decoded into a valid
sequence of levels using a penalised beam search.

The exact specification of this labelling stage is as follows:

  1. 1.

Given VB detections from the previous stage, a 2-D bounding cube is fit
tightly around the union each of detection’s slicewise series of polygons. By
expanding this cube by 100% on each side in the axial and coronal directions
and by 50% on each side in the sagittal direction a volume around each
detection and its nearby anatomical features is created.

  2. 2.

All extracted volume are resized by cubic interpolation to $224\times
224\times 16$ (where 16 is the number of sagittal slices). These are then fed
to an appearance network which outputs a 23-element probability vector,
attempting to classify the VB as a level from S1-C3.

  3. 3.

The output probability vectors are re-calibrated [guo_calibration_2017] using
a softmax temperature $T=10$. These recalibrated vectors are then used to
construct an initial probability-height map, $P\in\mathbb{R}^{H\times 23}$.
For a detection polygon with output appearance probability vector
$\mathbf{p}_{a}\in\mathbb{R}^{23}$ and spanning from height $h_{1}$ to $h_{2}$
in the scan, $P$ is equal to this value between these heights, i.e.
$P(h)=\mathbf{p}_{a}\>\forall\>h_{1}\leq h\leq h_{2}$.

  4. 4.

$P$ is then given as input to a image-to-image context CNN. This outputs
$P^{\prime}\in\mathbb{R}^{H\times 23}$, a refined probability-height map which
considers the neighbours of each detection to update its probability vector.
This refinement step can be seen in Figure 3.

  5. 5.

The next step is to decode $P^{\prime}$ into a series of discrete level
predictions. Firstly, the probability vector at the centroid of each
detection, $\mathbf{p}_{i}$ is extracted. The joint probability of a sequence
of levels can then be calculated as the product of elements of these vectors.
For example, if detected VB $i$ has a probability of being S1 of $a$ and VB
$j$ has probability of being L5 of $b$, then the joint probability of the
level sequence is $a\times b$. To reduce search space, we use a beam search to
find the most likely sequence, starting from the bottom detection and only
storing the 100 most probable level sequences at each step.

  6. 6.

To impose constraints on sequences of level predictions, penalties can be
added to the sequence scores. For example, sequences which predict the same
level twice are given a probability of 0. Furthermore, if there is a missed
detection (i.e. L4 directly above S1), the joint probability score is
multiplied by a penalty score to reflect the rarity of this event. To allow
for cases where S2 or C2 are detected, double detections of S1 or C3 are
allowed without a score penalty and then relabelled in post-processing.

At this point, each vertebra in the scan from S1 to C3 should have been
detected and assigned a level label. These detections can be used in a range
of applications. For example in whole spine scans, SpineNet can be used to
measure spinal curvature (for scoliosis measurement) or extract regions of
anatomical interest such as the spinal chord or vertebral bodies for lesion
segmentation (e.g. spinal metastases or ankylosing spondylitis).

One area of particular interest is radiological grading of degenerative
changes. The following section describe how this is done in the context of the
SpineNet.

##  0.3 Extracting Intervertebral Volumes

To perform radiological grading in T2 lumbar scans, first volumes surrounding
each intervertebral disc must be extracted. This is done using the VB
detections from the previous stage. The mid-point between the centroids of two
consecutive VB detections is calculated. This defines the centre of each
extracted IVV. From this the volume is rotated such that the lower endplate of
the upper vertebra is horizontal. The width of the extracted endplate volume
is then defined as double that of the larger VB detection. The height of the
IVV is then chosen such that the aspect ratio of the extracted patch is 2:1.
This is done across all slices in which the vertebrae is detected.

##  0.4 Radiological Grading of Intervertebral Volumes

Once the IVVs are extracted, they are resampled to have a resolution of
$112\times 224\times 9$ (height, width and sagittal slices respectively) These
volumes are then fed to a radiological grading network which outputs scores
for the following radiological gradings:

  1. 1.

Pfirrman Grading (5 classes)

  2. 2.

Disc Narrowing (4 classes)

  3. 3.

Central Canal Stenosis (4 classes)

  4. 4.

Upper & Lower Endplate Defects (Binary)

  5. 5.

Upper & Lower Marrow Changes (Binary)

  6. 6.

Left & Right Foraminal Stenosis (Binary)

  7. 7.

Spondylolisthesis (Binary)

  8. 8.

Disc Herniation (Binary)

The radiological grading network is a standard multi-class classification 3D
CNN. We experiment with a range of different architectures, described in
Section LABEL:sec:results. Each of these grading schemes is described in
further detail in Section LABEL:sec:datasetsGenodisc. Note that the same model
is used for all gradings schemes and vertebral levels.