Realistic Injection Simulations of a Cyclotron Spiral Inflector using OPAL

We present an upgrade to the particle-in-cell ion beam simulation code OPAL that enables us to run highly realistic simulations of the spiral inflector system of a compact cyclotron. This upgrade includes a new geometry class and field solver that can handle the complicated boundary conditions posed by the electrode system in the central region of the cyclotron both in terms of particle termination, and calculation of self-fields. Results are benchmarked against the analytical solution of a coasting beam. As a practical example, the spiral inflector and the first revolution in a 1 MeV/amu test cyclotron, located at Best Cyclotron Systems, Inc., are modeled and compared to the simulation results. We find that OPAL can now handle arbitrary boundary geometries with relative ease. Comparison of simulated injection efficiencies, and beam shape compare well with measured efficiencies and a preliminary measurement of the beam distribution after injection.


I. INTRODUCTION
OPAL [1] is a particle-in-cell (PIC) code developed for the simulation of particle accelerators. It is highly parallel and comes in two distinct flavors: OPAL-CYCL (specific to cyclotrons and rings) and OPAL-T (general purpose). For the presented application, we focus on OPAL-CYCL which has been used very successfully to simulate existing high intensity cyclotrons like the PSI Injector II [2], and PSI Ring Cyclotron [2,3] as well as to design new cyclotrons like CYCIAE [4], DAEδALUS [5,6], and IsoDAR [7,8].
However, one piece has been missing so far: the axial injection using a spiral inflector or an electrostatic mirror. Both are electrostatic devices that bend the beam from the axial direction into the mid-plane of the cyclotron where it is subsequently accelerated. A schematic view of this type of injection for a spiral inflector is shown in Figure 1.
In order to enable OPAL-CYCL to track particles through a complicated set of electrodes like a spiral inflector, the following additions have been made: • The Smooth Aggregation Algebraic Multi Grid (SAAMG) solver [9] has been extended to include arbitrarily shaped boundaries.
• A geometry class has been implemented that provides boundary conditions to the field solver and handles particle termination should their trajectories intersect with the electrodes.
FIG. 1. Schematic of a spiral inflector with particle trajectories from an OPAL simulation. The beam enters axially (from the top) through an aperture (grey) and is bent into the mid-plane by a combination of the electrostatic field generated by the spiral electrodes (green and blue, voltage typically symmetric at +V spiral and −V spiral ) and the cyclotron's main magnetic field. Then it is accelerated by the two Dees (copper, Dummy-Dees not shown). Color online.
• A number of improvements have been made to the internal coordinate transformations and the handling of beam rotations in OPAL-CYCL in order to accommodate the injection off-mid-plane.
These additions will be discussed in Section II and references therein. At the end of this section, a benchmark against the analytical solution of a coasting beam in a grounded beam pipe with variable axial offset will be presented.
A specific example of the usefulness of realistic spiral arXiv:1612.09018v1 [physics.acc-ph] 29 Dec 2016 inflector simulations is the ongoing R&D effort for the DAEδALUS [5,6], and IsoDAR [7,8] experiments (described briefly in Section III). In both cases, a very high intensity beam (≈ 10 − 30 mA average) needs to be injected, which is higher than current state-of-the-art cyclotrons have demonstrated. Part of the R&D for these projects was an experiment in collaboration with Best Cyclotron Systems, Inc. (BCS) in Vancouver, Canada to test a flat-field ECR ion source, transport through the Low Energy Beam Transport System (LEBT) and finally injection into a small test cyclotron through a spiral inflector. The results of this campaign are described in much detail in [10] and the most important points will be reiterated in Section III, before OPAL simulation results are benchmarked against the experimental results.

II. THE PARTICLE-IN-CELL CODE OPAL
For this discussion we briefly introduce OPAL-CYCL [2], one of the four flavours of OPAL.

A. Governing equation
The collision between particles can be neglected in the cyclotron under consideration, because the typical bunch densities are low. The general equations of motion of charged particles in electromagnetic fields can be expressed in time domain by where m 0 , q, γ are rest mass, charge and the relativistic factor. With p = m 0 cγβ we denote the momentum of a particle, c is the speed of light, and β = (β x , β y , β z ) is the normalized velocity vector. The time (t) and position (x) dependent electric and magnetic vector fields are written in abbreviated form as E and B. If p is normalized by m 0 c, Eq. (1) can be rewritten component wise in Cartesian coordinates as The evolution of the beam's distribution function can be expressed by a collisionless Vlasov equation: Here we have assumed that M particles of the same species are within the beam. In this particular case, E and B include both external applied fields and space charge fields, all other fields are neglected.

B. Bunch rotations
In OPAL-CYCL, the coordinate system is different from OPAL-T: The x and y coordinates are the horizontal coordinates (the cyclotron mid-plane) and z is the vertical coordinate. Internally, both Cartesian (x, y, z) and cylindrical (r, Θ, z) coordinate systems are used. For simplicity, in the past, the injection of bunches in OPAL-CYCL had to happen on the cyclotron mid-plane with the option to include an offset in z direction within the particle distribution itself. The global coordinates of the beam centroid were thus restricted to r and Θ. In order to accommodate the injection of a bunch far away from the mid-plane and a mean momentum aligned with the z-axis, the handling of rotations in OPAL-CYCL was updated from 2D rotations in the mid-plane to arbitrary rotations/translations in three-dimensional (3D) space. Quaternions [11] were chosen to avoid gimbal-lock [12] and a set of rotation functions was implemented. These are now used throughout OPAL-CYCL. Outside of the new SPIRAL mode, quaternions are also used to align the mean momentum of the bunch with the y axis before solving for self-fields in order to simplify the inclusion of relativistic effects in the calculation of self-fields. In the SPIRAL mode, no relativistic effects are taken into account due to the low injection energy (typically < 100 keV).

C. Fieldsolver
The space charge fields can be obtained by a quasistatic approximation. In this approach, the relative motion of the particles is non-relativistic in the beam rest frame, so the self-induced magnetic field is practically absent and the electric field can be computed by solving Poisson's equation where φ and ρ are the electrostatic potential and the spatial charge density in the beam rest frame. The electric field can then be calculated by and back transformed to yield both the electric and the magnetic fields, in the lab frame, required in Eq. (3) by means of a Lorentz transformation. As mentioned in the previous section, this step is omitted in the SPIRAL mode of OPAL-CYCL by setting γ = 1.
A parallel 3D Poisson solver for the electrostatic potential computation in a round beam pipe geometry with open-end boundary conditions is presented in [9].
In this section, we present a method for solving the space charge Poisson problem (6) in much more complex geometries of particle accelerators. The problem discretization takes into account the complex geometries of beam pipe elements and the space charge forces are computed within the geometry. This assures that the space charge components are taken properly into account in any type of geometry.
where Ω represents the computational domain, and ∂Ω the boundary surface of the geometry. Our approach solving the Poisson problem in complex geometries includes the unstructured and structured meshes. Unstructured meshes are used to define the complex geometry because they are easy to conform a block to a complicated shape, and structured meshes are used for the finite-difference discretization method.

Spatial discretization
To solve the Poisson problem (6) on Ω, we use a cell-centered discretization for Laplacian. A grid point is called interior if all its neighbors are in Ω, or nearboundary point otherwise. The discrete Laplacian with mesh spacing h on interior points is defined as the regular seven-point finite difference approximation of −∆.
where e i is the ith coordinate vector [13]. At near-boundary point, x in Figure 2 where not all its neighbors are in the domain, we use an approximation for Laplacian and the value ofφ at near-boundary point is obtained by one of the extrapolation method described in [9]. The approximation is based on linear extrapolation of near-boundary points (0 < s ≤ 1): The value ofφ at x L is defined through one of the extrapolation method in Eq. (9).

Implementation
The query class, an interface to search for the points inside of the irregular domain is implemented. Once the points inside the irregular domain are detected, their intersection values in six different directions are stored in containers. The coordinates values are mapped into its intersection values to be used as a fast look-up table. The distances between the near-boundary point and its intersection values are used for the linear extrapolation.
The finite difference approximation of the Poisson problem requires solving a system of linear equations to compute the electrostatic potential. The resulting linear system is solved using the preconditioned conjugate gradient algorithm complemented by an algebraic multigrid preconditioner using the Trilinos framework [14].
Trilinos is a collection of software packages that support parallel linear algebra computations on distributed memory architectures, in particular the solution of linear systems of equations. Epetra provides the data structures that are needed in the linear algebra libraries. Amesos, AztecOO, and Belos are packages providing direct and iterative solvers. ML is the multi-level package, that constructs and applies the smoothed aggregationbased multigrid preconditioners.

D. External fields
With respect to the external magnetic field two possible situations can be considered: in the first situation, the real field map is available on the median plane of the existing cyclotron machine using measurement equipment. In most cases concerning cyclotrons, the vertical field, B z , is measured on the median plane (z = 0) only. Since the magnetic field outside the median plane is required to compute trajectories with z = 0, the field needs to be expanded in the Z direction. According to the approach given by Gordon and Taivassalo [16], by using a magnetic potential and measured B z on the median plane at the point (r, θ, z) in cylindrical polar coordinates, the 3rd order field can be written as The four stages of OPAL geometry preparation: (a) Initial CAD model of the electrodes in the system. (b) An inverted solid is created by subtraction of all elements in (a) from a "master volume" (in this case a cylindrical outer chamber). (c) The inverted geometry is saved in stp format and a mesh is generated using GMSH [15]. (d) During the initialization phase, OPAL creates a voxel mesh to speed up the tests that have to be performed every time-step (cf. text).
All the partial differential coefficients are computed on the median plane data by interpolation, using Lagrange's 5-point formula. In the other situation, the 3D field for the region of interest is calculated numerically by building a 3D model using commercial software during the design phase of a new cyclotron. In this case the calculated field will be more accurate, especially at large distances from the median plane, i.e. a full 3D field map can be calculated. For all calculations in this paper, we use the second method. Fields (where applicable) are generated as 3D field maps with VectorFields OPERA [17]. In case of RF fields, OPAL-CYCL varies the field with a cosine function: where t is the time of flight, ω rf the RF frequency, and φ S the starting phase of the particle. Finally, in this paper, both the external fields and selffields are used to track particles during each time step using a 4 th order Runge-Kutta (RK) integrator, in which the fields are evaluated four times in each time step. Space charge fields are assumed to be constant during one time step, because their variation is typically much slower than that of external fields.

E. Geometry
For the simulation of precise beam dynamics, an exact modelling of the accelerator geometry is essential. Usually a CAD model of the accelerator or part of it is already available (Figure 3(a)). From these models we need to create the vacuum chamber or beam-pipe for specifying the boundary, i.e. the simulation domain (Ω in Eq. 6). In case of the spiral inflector, we have to add a cylinder to limit the vacuum chamber (figure Figure 3(b)). This modified CAD model can be used to create a triangle mesh T modeling the vacuum chamber of the inflector (Figure 3(c)).
Facilitating a meshed vacuum chamber or beam-pipe, OPAL is able to model arbitrary accelerator geometries and provides methods for 1. testing whether a particle will collide with the inner surface of the geometry (boundary, ∂Ω) in the next time step 2. computing the distance d = |x 0 − I| from a given point x 0 to the boundary intersection point I with ∂Ω (c.f. Figure 4) 3. testing whether a given point x 0 is inside the geometry. Only points inside the geometry are in the computational domain Ω.
The geometry can consist of multiple parts. Each part must be modeled as a 3D closed volume. The used methods are based on well known methods in computer graphics, especially ray tracing [18].

Initializing the geometry
For testing whether a particle will collide with the boundary in then next time step, we can run a line segment/triangle intersection test for all triangles in the mesh. Even to be able to model simple structures, triangle meshes with thousands of elements are required. Applying a brute force algorithm, we have to run this test for all particles per time-step, rendering the naive approach as not feasible due to performance reasons.
In computer graphics this problem is efficiently solved by using voxel meshes. A voxel is a volume pixel representing a value on a regular 3D grid. Voxel meshes are used to render and model 3D objects.
To reduce the number of required line segment -triangle intersection tests, a voxel mesh V covering the triangle mesh T is created during initialization of OPAL. In this step, all triangles are assigned to their intersecting voxels. Whereby a triangle usually intersects with more than one voxel.
For the line segment/triangle intersection tests we can now reduce the required tests to the triangles assigned to the voxels intersecting the line segment. The particle boundary collision test can be improved further by comparing the particle momentum and the inward pointing normal n of the triangles.
In the following, we use the following definitions: T : represents the set of triangles in the triangulated mesh.
V : represents the set of voxels modelling the voxelized triangle mesh.
L: a closed line segment bounded by points x 0 and x 1 (c.f. figure Figure 4).
R: a ray defined by the starting point x 0 passing through x 1 .
T v ⊂ T : represents the subset of triangles t ∈ T which have intersections with v ∈ V .
V L ⊂ V : represents the subset of voxels v ∈ V which intersections with the line segment L.
I t,L : represents an intersection point of a line segment L with a triangle t ∈ T .
T L : represents the set of tuples (t L , I t,L ) with t L ∈ T intersects with L.

Basic ray/line-segment boundary intersection test
In the first step, we have to compute V L , which is the set of voxels in the voxel mesh V which have intersections with the given line-segment or ray L. With we can compute a small subset of triangles which might have intersections with L. We have to run the ray/line-segment triangle test only for all triangles in T V,L .

Particle boundary collision test
To test whether a particle will collide with the boundary we have to test whether the line-segment given by the particle position x 0 at time step s and the expected particle position x 1 at time s+1 intersect with the boundary. The closed line segment L given by x 0 and x 1 is used as input parameter for the boundary collision test.

boundaryCollisionTest(L)
compute T L with basicIntersectionTest() In the collision test we use a slightly modified intersection test. If the particle moves away from a given triangle t, we do not have to run the line-segment triangle intersection test for t. Since we know the inward pointing normals n for all triangles, hence we can compute the dot product of the normal n, and the vector defined by the current particle position x 1 and the position in the next time step x 1 . If both vectors point in the same direction, so there cannot be a collision of the particle with the boundary.

Compute distance from point to boundary
To compute the distance from a point P inside the geometry to the boundary in a given direction v, the same algorithm as for the particle boundary collision test can be used with the ray L = (P, v).

Inside test
Test whether a given point P is inside the geometry.

isInsideTest(P)
select R outside geometry

F. Simple Test Case
In order to test the proper functionality of the SAAMG fieldsolver in combination with an external geometry file, we compared the calculated fields and potentials of a FLATTOP distribution (uniformly populated cylinder) inside a geometry generated according to the previous subsection with the analytical solution of a uniformly charged cylindrical beam inside a conducting pipe. The geometry file we used contained a 1 m long beam pipe with 0.1 m radius. For a concentric beam, a regular Fast Fourier Transform (FFT) field solver is sufficient. However, if the beam is off-centered by an amount ξ (see Figure 5) and especially, when it is close to the conducting walls of the beam pipe, the electric field calculated by the FFT solver does no longer reproduce reality. This is even more true for complicated geometries like a spiral inflector. To compare the simulated results with the analytical solution presented below, all simulations were run in such a way that the bunch frequency f b was adjusted such that for a given bunch length l b and beam velocity v b the given beam current I corresponded to the equivalent DC beam current (i.e. subsequent bunches are head-totail): The bunch radius r b was chosen to be 0.01 m and the length l b to be 0.6 m so that in the center of the bunch to very good approximation the conditions of an infinitely long beam hold. For such an infinitely long beam, E and Φ are independent of z, E z (x, y) = 0, and E x,y (x, y) and Φ(x, y) can be calculated from Poisson's equation using the method of image charges. With (where I is the beam current and v b the beam velocity), the resulting expressions for inside and outside (super-FIG. 6. χ 2 red of the calculated Φ compared to the simulated values. It can be seen that the results are slightly worse for a beam close to the beam pipe. Also, for number of mesh cells in x-direction > 256, only marginal improvement can be seen. These simulations were performed with a total number of particles NP = 8000 · MX where MX denotes the number of mesh cells in x-direction (the independent variable in the plot). script "in" and "out") of the beam envelope are then: where A wide parameter space was mapped in terms of mesh size, number of particles, beam length, and position of the beam inside the beam pipe. As can be expected, the comparison between theory and simulation gets better with higher resolution (i.e. higher number of mesh cells), and larger number of particles. The reduced χsquare was chosen to compare the simulated results to the theoretical prediction and a plot of χ 2 red for variation FIG. 7. Φ and Ex along the x-axis for different offsets ξ. Dots are values calculated by OPAL, while solid lines are the analytical solution. Excellent agreement can be seen with a reduced 0.01 < χ 2 red < 0.02 for the potential and 0.03 < χ 2 red < 1.3 for the field. of mesh size and particle number is shown in Figure 6. It can be seen that the agreement is better for a centered beam and so it is especially important to choose a high enough resolution when the beam is close to the beam pipe (or other electrodes in the system). For this particular case, it was found that a total number of mesh cells of 256 in x-direction (≈ 25 across the beam diameter) together with ≈ 2 · 10 6 particles gave excellent agreement even when the beam was touching the pipe, with only slight or no further improvement at larger numbers. As another representative example, the OPAL results of a 0.6 m long beam in a 10 cm diameter beam pipe, using 2048000 particles and a mesh of dimensions 256 x 128 x 512, are plotted together with the analytical solution from equations 15 -20 for a beam with varying offset ξ in x-direction from the center of the beam pipe in Figure 7.
In summary, the SAAMG solver performed as expected when tested with the simple test-case of a quasi-infinite uniform beam in a conducting pipe. In the next section, the solver will be applied to a real world problem and results will be compared to measurements.

III. BENCH-MARKING AGAINST EXPERIMENTS
An important step in bench-marking new simulation software is the comparison with experiments. During the summers of 2013 and 2014, a measurement campaign was held at Best Cyclotron Systems Inc. (BCS), to test the production of a high intensity H + 2 beam in an off-resonance ECR ion source and its injection into a compact (test) cyclotron through a spiral inflector. These tests were performed within the ongoing R&D effort for the DAEδALUS and IsoDAR experiments (cf. next section) and provided a good opportunity to compare results of injected beam measurements with OPAL simulations using the new spiral inflector capability.

A. DAEδALUS and IsoDAR
The Decay At-rest Experiment for δ CP studies At a Laboratory for Underground Science (DAEδALUS) [5,6] is a proposed experiment to measure CP violation in the neutrino sector. A schematic view of one DAEδALUS complex is shown in Figure 8. H + 2 is produced in an ion source, transported to the DAEδALUS Injector Cyclotrons (DIC), and accelerated to 60 MeV/amu. The reason for using H + 2 instead of protons is to overcome space charge limitations of the high required beam intensity of 10 emA of protons on target. H + 2 gives 2 protons for each unit of charge transported, thus mitigating the risk. The ions are subsequently extracted from the cyclotron and injected into the DAEδALUS Superconducting Ring Cyclotron (DSRC) where they are accelerated to 800 MeV/amu. During the highly efficient stripping extraction, the 5 emA of H + 2 become 10 emA of protons which impinge on the neutrino production target (carbon) producing a neutrino beam virtually devoid ofν e . In a large detector, one can then look forν e appearance through neutrino oscillations. As is depicted in Figure 8, the injector stage of DAEδALUS can be used for another experiment: The Isotope Decay At Rest experiment Iso-DAR [7,8]. In IsoDAR, the 60 MeV/amu H + 2 will impinge on a beryllium target creating a high neutron flux. The neutrons are captured on 7 Li surrounding the target. The resulting 8 Li beta-decays producing a very pure, isotropic ν e beam which can be used forν e disappearance experiments. IsoDAR is a definitive search for so-called "sterile neutrinos", proposed new fundamental particles that could explain anomalies seen in previous neutrino oscillation experiments.
At the moment, OPAL-CYCL is used for the simulation of three very important parts of the DAEδALUS and Iso-DAR systems: 1. The spiral inflector 2. The DAEδALUS Injector Cyclotron (DIC), which is identical to the IsoDAR cyclotron.

The DAEδALUS Superconducting Ring Cyclotron (DSRC) for final acceleration.
For the topic of bench-marking, we will restrict ourselves to item 1., the injection through the spiral inflector.

B. The Teststand
As mentioned before, the results of the injection tests are reported in detail in [10]. Here, we will summarize the items pertinent to a comparison to OPAL, specifically, how we obtain the particle distribution at the end of the LEBT (entrance of the cyclotron), used as initial beam in the subsequent injection simulations with the SAAMG solver.
The test stand was comprised of the following parts: 1. Versatile Ion Source (VIS) [20]. An off-resonance Electron Cyclotron Resonance (ECR) ion source.

Cyclotron with spiral inflector.
During the experiment, it was possible to transport up to 8 mA of H + 2 as a DC beam along the LEBT to the cyclotron and transfer 95% of it through the spiral inflector onto a paddle probe. The 4-rms normalized emittances stayed below 1.25 π-mm-mrad. Capture efficiency into the RF "bucket" was 1-2% because of reduced dee voltage (V dee ) due to an under-performing RF amplifier (cf. discussion in [10]).

C. Initial Conditions
The quality of any simulation result depends on the initial conditions. In the case of the OPAL SAAMG simulations of the injection through the spiral inflector, the initial particle distribution consisted of 66021 particles obtained by carefully simulating the ion source extraction (using KOBRA-INP [21]) and the subsequent LEBT (using WARP [22]) and comparing the simulation results to the measurements, with good agreement as reported in [10]. During the WARP simulations of the LEBT, the "xy-slice-mode" was used in which the selffields are calculated only for the transverse direction (assuming only very slow changes in beam diameter along the z-axis compared to the length of each simulation step) and neglecting longitudinal self-fields (which is a sensible approach for DC beams). Space charge compensation played a big role in order to obtain good agreement and was taken into account using a semi-analytical formula [23]. The final particle distribution that was obtained for the set of parameters recorded during the measurements is shown in Figure 9. It should be noted that the bunch was generated from the xy-slice at a position 13 cm away from the cyclotron mid-plane and coaxial with the cyclotron center by randomly backward and forward-projecting particles according to their respective momenta. It can be seen that this beam enters the spiral inflector converging, which has been found experimentally to give the best injection efficiency. The important parameters of the injected beam are listed in Table I.

D. Results
The beam described in the previous section was then transported through the spiral inflector using the standard FFT and the new SAAMG field solver described in Section II C, and putting in place the geometry seen in Figure 3. Inside the cyclotron, 45 • after the exit of the spiral inflector, a radial probe was placed which had 5 fingers of ≈ 5 mm vertical extent, and ≈ 1 mm radial extent each. On these fingers, the electrical beam current was measured. The probe was slowly moved from a position blocking the beam completely, to just outside of the radial extent of first turn, thereby giving the beam current distribution shown in the top plot of Figure 11. In the same plot, the results from OPAL sim- ulations using the same parameters as recorded during the measurement (see Table I) are plotted. Good qualitative agreement can be seen for both the FFT and the SAAMG solver. Due to the tail towards low radius, two Gaussians are used to fit each peak in the left column.
The full widths at half maximum (FWHM) of the dominant Gaussians are listed in Table II. There seems to be a slight shift towards higher vertical position that is better reproduced in the SAAMG solver, but this is well within the systematic errors of the measurement and how well the initial parameters like magnetic field, spiral inflector voltage and beam distribution were known, hence the conclusion is that both FFT and SAAMG reproduce the measured radial-vertical beam distribution equally well for a 6 mA beam. This shows that the SAAMG solver is working as expected.

Higher beam currents
For initial beam currents up to 36 mA, the results of the FFT and SAAMG solvers start to show stronger (but still fairly subtle) discrepancies which can be attributed to the image charges on the electrodes only included with the SAAMG solver. An example is shown in Figure 10 where the expected spreading of the peak positions is accompanied by a noticeable overall shift towards smaller radii in case of the SAAMG solver.

A. Summary and Discussion
For the first time a comprehensive and precise beam dynamics simulation model, from the ion source, throughout the LEBT, into the central region of a cyclotron was presented. The central region includes the spiral inflector and the first accelerating gap. From the exit of the LEBT, the open source code OPAL was used for the beam transport through the spiral inflector and the first turn of the test cyclotron. Key ingredients of the model are the the flexible handling of the complex  . Measurement and simulation of a 5-finger probe (for description cf. to text) in the first turn of the BCS test cyclotron central region, 45 • after the exit of the spiral inflector. Each finger has a main peak and a tail towards lower radius. These are the particles that are not sufficiently accelerated. To guide the eye and to compare with the simulations, the data for each finger is fitted with two Gaussians. geometry, and the field solvers for space charge calculation. In comparison with first measurements, both the FFT and the SAAMG solver perform well, with hints that image charge effects become more important at higher currents, where use of the SAAMG solver allows including the complicated boundary conditions posed by the electrode system. These ingredients -geometry and field solvers (FFT and SAAMG) are now included with OPAL. Validation of the model included simple test cases and comparison to measurements from a dedicated cyclotron test stand, injecting a DC beam of 7 mA of H + 2 . Both yielded good agreement. The level of detail available in this model now allows us to obtain a detailed understanding, and predict the complicated beam dynamics in the high current compact IsoDAR cyclotron.

B. Outlook
Given this benchmarked model, a full start to end simulation of the IsoDAR cyclotron is ongoing, using first the detailed geometry for acceleration up to 1 MeV/amu, and then a simplified model with FFT only for the subsequent acceleration up to 60 MeV/amu. Recently, a proposal was put forward to test direct injection into the compact IsoDAR cyclotron using a Radio Frequency Quadrupole (RFQ) [24]. For the design of this device, OPAL-CYCL and the new SPIRAL mode will play an essential role.