Realistic simulations of a cyclotron spiral inflector within a particle-in-cell framework

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. 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 midplane of the cyclotron where it is subsequently accelerated.A schematic view of this type of injection for a spiral inflector is shown in Fig. 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: (i) The Smooth Aggregation Algebraic Multi Grid (SAAMG) solver [9] has been extended to include arbitrarily shaped boundaries.(ii) 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.(iii) A number of improvements have been made to the internal coordinate transformations and the 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 midplane 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.
handling of beam rotations in OPAL-CYCL in order to accommodate the injection off-midplane.These additions will be discussed in Sec.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 inflector simulations is the ongoing R&D effort for the DAEδALUS and IsoDAR experiments (described briefly in Sec.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 Sec.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 flavors of OPAL.OPAL uses the PIC method to solve the collisionless Vlasov equation (cf.Sec.II A) in the presence of external and self-fields (cf.Sec.II B and Sec.II C).Particles are tracked using a 4th order Runge-Kutta (RK) integrator, in which the external fields are evaluated four times in each time step.Self-fields are assumed to be constant during one time step, because their variation is typically much slower than that of external fields.Boundary conditions for the field solver and particle termination are provided by the geometry class (cf.Sec.II D).

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 the time domain by with p ¼ m 0 cγβ the mechanical momentum of the particle and m 0 , q; γ; c, the rest mass, charge, relativistic factor, and speed of light, respectively β ¼ ðβ x ; β y ; β z Þ is the normalized velocity vector.The time and position dependent electric and magnetic vector fields Eðx; tÞ and Bðx; tÞ are written as E and B for brevity.
For M particles and the canonical variables: x j the position and P j ¼ p j − qA j the momentum of the jth particle, the evolution of the beam's distribution function fðx j ; P j ; tÞ∶ ðℜ 3M × ℜ 3M × ℜÞ → ℜ can be expressed by a collisionless Vlasov equation where A denotes the vector potential.In this particular case, the vector fields E and B, include both external (applied) fields and space-charge (self) fields, all other fields are neglected.

B. External fields
With respect to the external magnetic fields there are two options: (1) OPAL can calculate the off-median plane (z ≠ 0) components from a field map obtained for the cyclotron midplane, (2) a 3D fieldmap can be calculated externally and loaded into OPAL.For external electric fields only option two is available.
In option one, the vertical magnetic field component B z alone is known in the median plane(z ¼ 0) only.The reason for this mode is that often the cyclotron field is obtained by measurement and this simplifies the process.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.Using a magnetic potential and (measured) B z on the median plane at the point ðr; θ; zÞ in cylindrical coordinates, the 3rd order field can be written as with B z ≡ B z ðr; θ; 0Þ and All the partial differential coefficients are computed from the median plane data by interpolation using Lagrange's 5-point formula.
In option two, a full 3D field map for the region of interest is calculated numerically by 3D CAD modeling and finite elements analysis (FEA) using commercial software.In this case the calculated field will be more accurate, especially at large distances from the median plane, but not necessarily correspond as well to the field of an existing machine which might include manufacturing errors.For all considerations in this paper, we use the second method.Fields (where applicable) were generated as 3D field maps with VectorFields OPERA [11].
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.

C. Space-charge fields
The space-charge fields can be obtained by a quasistatic approximation.In this approach, the relative motion of the particles is nonrelativistic 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 in the computational domain Ω with (arbitrary) boundary surface ∂Ω.Here, ϕ 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 (E sc , B sc ) in the lab frame, required in Eq. (1), by means of a Lorentz transformation.As beams through a spiral inflector are typically nonrelativistic, this step is omitted in the SPIRAL mode of OPAL-CYCL by setting γ ¼ 1.
A parallel 3D Poisson solver for the electrostatic potential computation was presented in [9] for the special case of Ω and ∂Ω being a round beam pipe with open-end boundary conditions.In this section, we present a method for solving Eq. ( 2) in more complex geometries by taking into account arbitrary boundary shapes ∂Ω during the discretization of the problem and the computation of the space-charge forces within Ω.

Spatial discretization
To solve the Poisson problem in Eq. (2) in Ω, we use a cell-centered discretization for the Laplacian.A grid point is called interior point if all its neighbors are in Ω, or near-boundary point otherwise.The discrete Laplacian with mesh spacing h ¼ h i ¼ ðh x ; h y ; h z Þ for interior points is defined as the regular seven-point finite difference approximation of −∇ 2 , where êi is the ith coordinate vector and h i denotes the mesh spacing in the ith dimension [12].At a near-boundary point (e.g., x in FIG. 2) where not all its neighbors are in the domain, we use an approximation for the Laplacian and the value of ϕ at the near-boundary point is obtained by one of the extrapolation methods described in [9].With s L=R;i denoting the six possible cases for a boundary near x ([Left, Right] × [x, y, z]) and 0 < s L=R;i ≤ h i , we have where the value of ϕ at x L=R;i is defined through one of the following extrapolation methods:

Implementation
The query class, an interface to search for the points inside of the irregular domain was implemented.Once the points inside the irregular domain are detected, their intersection values in six different directions are stored in containers.The coordinate values are mapped onto their 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 [13].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 multilevel package, that constructs and applies the smoothed aggregation-based multigrid preconditioners.

D. Geometry
For the precise simulation of beam dynamics, exact modeling of the accelerator geometry is essential.Usually a CAD model of the accelerator (or part of it) is already available [see Fig. 3(a)].From these models we need to create a mesh, to specify the boundary, ∂Ω in Eq. ( 2).In case of the spiral inflector, we have to add a cylinder as an outer boundary [see Fig. 3(b)].This modified CAD model can then be used to create a triangle mesh T modeling the inner surface of the system [see Fig. 3(c)].
Provided with a closed, meshed boundary, 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 ¼ jX 0 − Ij from a given point X 0 ∈ Ω to the boundary intersection point I ∈ ∂Ω (see FIG. 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 [15].

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 in the case of very simple structures, triangle meshes with thousands of elements are required.Applying a brute force algorithm, we have to run this test each timestep for all particles, 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 the 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 test 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 modeling the voxelized triangle mesh.L: represents a closed line segment bounded by points X 0 and X 1 (see FIG. 4).R: represents 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 have 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 Þ where 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 T v ¼ ft ∈ Tjt intersects with vg 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
In the collision test we use a slightly modified intersection 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 intersects 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.However, if the particle moves away from a given triangle t, we do not have to run the line-segment triangle intersection test for t.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 n • ðX 1 − X 0 Þ > 0, 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, given a directional vector d, the same algorithm as for the particle boundary collision test can be used with the ray L ¼ ðP; dÞ.

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

E. 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 midplane 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 [16] were chosen to avoid gimbal-lock [17] 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 through Lorentz backtransformation into the lab frame.In the SPIRAL mode, no relativistic effects are taken into account due to the low injection energy (typically < 100 keV).

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 Fig. 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 OPAL results with the analytical solution presented below, all simulations were run in a way where 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-to-tail) 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 ϕ and E inside and outside (superscript "in" and "out") of the beam envelope are and 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 of mesh size and particle number is shown in Fig. 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 × 128 × 512, are plotted together with the analytical solution from Eqs. (3) for a beam with varying offset ξ in x-direction from the center of the beam pipe in FIG. 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. BENCHMARKING AGAINST EXPERIMENTS
An important step in benchmarking 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 offresonance 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 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).
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 Fig. 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 due to the required high 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.The injector stage (ion source and DIC) of DAEδALUS can be used for another experiment: The Isotope Decay At Rest experiment IsoDAR [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 source which can be used for νe disappearance experiments.IsoDAR is a definitive search for socalled "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 IsoDAR systems: (1) The spiral inflector (2) The DAEδALUS Injector Cyclotron (DIC), which is identical to the IsoDAR cyclotron.(3) The DAEδALUS Superconducting Ring Cyclotron (DSRC) for final acceleration.For the topic FIG. 7. ϕ and E x 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 0.01 < χ 2 red < 0.02 for the potential and 0.03 < χ 2 red < 1.3 for the field. of benchmarking, 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) [19].An off-resonance electron cyclotron resonance (ECR) ion source.( 2 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 underperforming 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 [20]) and the subsequent LEBT (using WARP [21]) and comparing the simulation results to measurements, with good agreement as reported in [10].During the WARP simulations of the LEBT, the "xy-slice-mode" was used in which the self-fields are calculated only for the transverse direction (assuming only very slow changes in beam diameter along the zaxis 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 semianalytical formula [22].The final particle distribution that was obtained for the set of parameters recorded during the measurements is shown in Fig. 9.It should be noted that the bunch was generated from the xy-slice at a position 13 cm away from the cyclotron midplane 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 solvers described in Sec.II C, and putting in place the geometry seen in Fig. 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.10.Measurement and simulation of a 5-finger probe (for description cf.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 toward 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.
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 Fig. 10.In the same plot, the results from OPAL simulations 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 Fig. 11 where the expected spreading of the peak positions is accompanied by a noticeable overall shift toward smaller radii in case of the SAAMG solver.

IV. CONCLUSION A. Summary and discussion
A comprehensive spiral inflector model within the particle-in-cell framework OPAL has been presented that allows injection from the axial direction through the combined 3D electric and magnetic fields of the spiral inflector (or mirror) and the cyclotron main field, respectively.Key features of the model are the flexible handling of the complex geometry (for boundary conditions and particle termination) 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 better understanding of, 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) [23].For the design of this device, OPAL-CYCL and the new SPIRAL mode will play an essential role.

ACKNOWLEDGMENTS
FIG. 2. Representative example of a boundary left of x inx-direction.x L;x ¼ x − h x • êx , x B;x ¼ x − s L;x • êx , x and x R;x ¼ x þ h x • êxare the point outside of Ω, boundary point, nearboundary point, and interior point, respectively.Similar for the y and z-directions.

FIG. 3 .
FIG.3.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 chamber).(c) The inverted geometry is saved in stp format and a mesh is generated using GMSH[14].(d) During initialization, OPAL creates a voxel mesh to speed up the tests that have to be performed every time-step (cf.text).

FIG. 4 .
FIG. 4. The boundary ∂Ω is discretized by a set of triangles T. Shown is the line segment/triangle intersection test with the line segment ðX 0 ; X 1 Þ, the triangle t ∈ T and the intersection at I ∈ ∂Ω.

FIG. 8 .
FIG. 8. Schematic of the DAEδALUS facility.H þ 2 is produced in the ion source, transported to the DAEδALUS Injector Cyclotron (DIC) and accelerated to 60 MeV=amu.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, 5 emA of H þ 2 becomes 10 emA of protons which impinge on the neutrino production target.

FIG. 9 .
FIG.9.Initial distribution for injection through the spiral inflector, obtained by carefully simulating the LEBT.The beam is slightly hollow and has a clearly visible halo.Both are caused by overfocused protons, as discussed in[10].The length of the bunch corresponds to one full rf period at 49.2 MHz and injection energy of 62.7 keV, centered at the synchronous phase, thus, to first order, representing a DC beam.
This work was supported by the U.S. National Science Foundation under Grant No. 1505858 and the corresponding author was partly supported by the Bose Foundation.The research at PSI leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under Grant No. 290605 (PSI-FELLOW/COFUND).Furthermore, the authors would like to express their gratitude to Best Cyclotron Systems, Inc. in Vancouver, for hosting the 1 MeV cyclotron injection tests, and the INFN-LNS ion source group in Catania, for the loan of the VIS.

TABLE I .
Beam and cyclotron parameters for inflection and acceleration studies.

TABLE II .
Full width at half maximum (FWHM) for the measured distribution and the OPAL results using the FFT and the SAAMG solver.