Ring-mode benchmarks¶
Ring mode's element maps and optics against independent references: MAD-X
5.09 (TWISS and its exact-model PTC_TWISS), two reference codes, and
analytic results. Reference code A tracks with element maps whose bend
multipoles are expanded in Cartesian coordinates; reference code B with
exact elements whose bend multipoles are curvilinear; reference code C, used
for the stripping foil and for painting, tracks thick elements and scatters
particles in a foil one screened-Rutherford collision at a time. Every
reference is generated by a script that never imports HELIX and is stored
with a manifest (the SHA-256 of the data — of the inputs too for some — and
the options) in tests/ring/fixtures/refs/; the tests in tests/ring/ check
against the stored files, so they run without the reference codes.
The test ring is the simplified Fermilab Booster (24 cells, 474.2 m,
combined-function bends with a sextupole component, pole faces of half the
bend angle, 800 MeV protons), a public deck under a BSD licence
(tests/ring/fixtures/).
Element maps against reference code A¶
One set of 64 particles (|x|, |y| ≤ 3 mm, |δ| ≤ 2×10⁻³, |ζ| ≤ 0.3 m) through 24 single elements of reference code A: drifts, quadrupoles, sextupole, octupole, a thin multipole up to order 4 (normal and skew), combined-function bends in each of its bend models and integrators, with a field error, with linear and full pole faces at both bend signs, a straight bend, a cavity, frame rotations and shifts, a wedge and a linear edge.
| Agreement | Cases |
|---|---|
| 10⁻¹³ relative (10⁻¹⁴ m in ζ, a sum of path-length differences) | 21 |
| Set by rounding: the exact curved-bend formula, which rounds in both codes (each about 1.5×10⁻¹³ m from a 40-digit evaluation; 1.4×10⁻¹³ m apart); the reference's straight bend (10⁻¹¹ relative) and δ after a cavity (5×10⁻¹⁶), which HELIX evaluates without cancellation | 3 |
The plain-Python and compiled kernels give identical results bit for bit.
Symplecticity¶
The Jacobian of every node's map in canonical coordinates, by central differences in 50-digit arithmetic, satisfies |JᵀSJ − S| ≤ 2×10⁻¹⁶ for all 31 maps tested (the largest values come from sines and cosines stored as float64 parameters); a deliberately non-symplectic map (a kick that depends on px) shows 3.6×10⁻⁴ — the check sees it.
Simple Booster against MAD-X¶
RF voltage zero (4-D). Chromaticity per unit δ (MAD-X's per-PT values times β0).
| Quantity | HELIX | MAD-X TWISS |
MAD-X PTC_TWISS (model=2, method=6, nst=10, exact=true, fringe off) |
|---|---|---|---|
| Qx | 6.709996465 | 6.709996463 | 6.709996476 |
| Qy | 6.809997589 | 6.809997588 | 6.809997597 |
| Q′x | −8.6531 | −8.6505 | −8.6531 |
| Q′y | −8.0892 | −8.0892 | −8.0892 |
| Q′x, bends without K2 | −4.3918 | −4.3891 | −4.3918 |
| Q′y, bends without K2 | −21.8051 | −21.8051 | −21.8051 |
| α_c | 0.03367748970 | 0.03367748976 | |
| γt | 5.449167 | 5.449167 | |
| β_x, D_x at the start | 33.736455 m, 3.186395 m | 33.736455 m, 3.186395 m |
β, α, the phase advance and D at all 72 quadrupoles and cell ends agree with
MAD-X's table to 6.3×10⁻⁹ relative in β, 7×10⁻⁹ in α, 1.7×10⁻⁹ in the phase
advance (in units of the tune) and 8×10⁻¹⁰ relative in D (the tests check
10⁻⁶). (MAD-X's twiss, chrom gives −8.891 / −7.956 — its
finite difference over off-momentum linear maps — and is not the reference.)
Pole faces. The same ring with each pole-face model, and reference code A's own optics for the same model (its default bend model, 98 kicks per bend):
| Pole faces | Q′x | Q′y | Reference code A |
|---|---|---|---|
exact (MAD-X TWISS and PTC_TWISS) |
−8.6531 | −8.0892 | — |
| full | −8.2440 | −8.0270 | −8.2440, −8.0270 |
| linear | −7.0144 | −7.8338 | −7.0144, −7.8338 |
With RF (200 kV, h = 84): the synchrotron tune from the 6-D one-turn map, 0.023656, against the small-amplitude formula, 0.023650.
Thin-lens FODO ring¶
Twelve cells of thin quadrupoles (f = 4 m) and 3 m drifts: tunes (with the integer part), β and α at the start equal the analytic values to 10⁻¹², and the chromaticity equals the thin-lens sum −(1/4π) Σ β k to 10⁻⁶.
Against HELIX's own linear matrices¶
Every element of the Booster (its RF cavity at its own voltage), and vertical and negative bends, tilted, skew and shifted magnets, thin lenses and steerers, for a proton lattice and an H⁻ one: the ring nodes' Jacobian at the reference (50-digit arithmetic) equals HELIX's linac element matrix to 10⁻¹⁰. The whole turn agrees to 2×10⁻⁸, the floor of a float64 finite difference. Two differences are by design, and what both describe is compared:
- an off-centre magnet also makes dispersion through its chromatic focusing, which the linac matrix leaves to the orbit — its transverse block is compared, and its dispersion terms are checked against the closed form;
- a thin ring cavity has no transverse map (as MAD-X's); the linac gap adds its transverse RF defocusing (10⁻⁶ /m at the Booster's 9 kV) — the longitudinal rows are compared.
A pole face at zero angle with a gap acts through its fringe field (a weak defocusing across the bend plane, MAD-X's −h tan(−ψ)) in ring mode and in the linac's edge alike: the two element matrices agree to 10⁻¹⁰ there too. Every edge model reproduces MAD-X's tunes for the Booster with sector faces, a 10 cm gap and FINT 0.5 (Qy 6.6707484).
Against reference code B¶
The simple Booster as reference code B runs it, with the curvilinear model: CODATA 2022 constants (the masses formed from their values in kilograms), the reference's integration (sixth order, six steps per magnet), its RF frequency and voltage.
1000 particles, two turns (|x| up to 15 mm, |δ| up to 4.1×10⁻³), the largest difference over all particles after the second turn:
| x | y | c t | px | py | pt |
|---|---|---|---|---|---|
| 4.6×10⁻¹² m | 8.6×10⁻¹⁵ m | 8.6×10⁻¹³ m | 1.5×10⁻¹³ | 2.2×10⁻¹⁵ | 1.4×10⁻¹⁶ |
This is the reference's own floating-point rounding: against a 40-digit evaluation of the same maps, the reference is off by about 3×10⁻¹² m in x after one turn and HELIX by about 10⁻¹⁶ m. Part of it comes from the form of the reference's sector bend and curvilinear kicks, which subtract terms of the order of the bending radius (HELIX evaluates the same expressions rearranged); the rest is the rounding of its compiled code (fused multiply-adds, its math library). The difference has no structure in the particles' coordinates: a cubic polynomial in the six initial coordinates explains no more of it than it would of noise.
Single elements: the quadrupoles, the multipole and the linear edge with a gap agree to the last bits (3×10⁻¹⁷ m in c t); the exact drift to 1.1×10⁻¹⁵ m in c t (both codes round there, each within 8×10⁻¹⁶ m of the 40-digit value). The reference's sector bend rounds: 2.5×10⁻¹⁴ m in x with a sextupole at negative curvature, 3×10⁻¹⁰ m with octupole and decapole components, where HELIX is within 2×10⁻¹⁸ m of the 40-digit result. Running the reference's published settings (sextupoles at fourth order, five steps) instead of the parity settings moves its particles by 4×10⁻¹² m.
Optics. The reference's tunes and chromaticity from its own tracking, against the same model's optics in HELIX:
| Qx (fraction) | Qy (fraction) | Q′x | Q′y | |
|---|---|---|---|---|
| reference code B | 0.70999676 | 0.80999779 | −7.9983 | −7.9987 |
| HELIX, curvilinear model | 0.70999675 | 0.80999779 | −7.9983 | −7.9987 |
Why the curvilinear chromaticity is not MAD-X's. Two differences, each measured: linear pole faces (Cartesian bends with linear faces: −7.0144 / −7.8338 against −8.6531 / −8.0892 with exact ones), and the curvilinear quadrupole field inside the bends, which differs from the Cartesian field, to the first order beyond it (cubic terms of the Hamiltonian), by a normal sextupole −h·k1. Cartesian bends with linear faces and that sextupole added give −7.9983 / −7.9987, the curvilinear chromaticity (the curvilinear model gives −7.99827 / −7.99872 with six integration steps per bend and with sixty).
Backends¶
The compiled kernels on one CPU thread and on all of them agree bit for bit
over three Booster turns of 256 particles, and the plain Python (ring mode
without numba) agrees with them bit for bit over one turn. A particle that
leaves a map's domain (the square root of a negative number in a drift, say)
is lost as invalid at the same node, with the same coordinates, on every CPU
backend. The CUDA kernel, run in numba's CUDA simulator on a machine without
a GPU, agrees with the CPU to 10⁻¹⁵ relative; on a GPU,
scripts/verify_cuda_ring.py makes the comparison (GPU runs
— it has not yet run on one).
Storage at fixed energy¶
The turn loop of ring runs, single particles over many turns. The tunes of a tracked particle are found by frequency analysis (NAFF; Laskar 1993).
| Check | Result |
|---|---|
| Invariants: a small-amplitude particle in a thin-lens FODO ring, 10⁶ turns | its Courant–Snyder invariants change by 3.6×10⁻¹² and 3.3×10⁻¹² (rounding; no secular drift) |
| A FODO ring from reference code B's examples (drift 0.25 m, quadrupole of 1 m at k = +1 m⁻², drift 0.5 m, quadrupole of 1 m at k = −1 m⁻², drift 0.25 m; 100 periods) | NAFF tune 0.188242184, the cell's tune to 2×10⁻¹² (the reference rounds it to 0.1883); α_x −1.5905003549, β_x 2.8216194080 m against the reference's −1.59050035, 2.82161941; the thick-lens cell analytically to 10⁻¹² |
| Simple Booster, 1000 turns, one probe, RF on | NAFF tunes = the eigentunes of the 6-D one-turn map: Qx to 5×10⁻¹⁰, Qy to 2×10⁻⁹, Qs to 2×10⁻¹⁰; the 6-D Qx differs from the 4-D one by 8×10⁻⁶ — the synchrotron motion coupled in through the horizontal dispersion at the cavities (Qy by 3×10⁻¹¹) |
| The same, RF off | NAFF tunes = the 4-D optics tunes: Qx to 1×10⁻⁹, Qy to 2×10⁻⁹ |
| NAFF on pure tones | the frequency to rounding (10⁻¹⁶) from three periods in the record on, with a constant beside the tone or not, and for a real signal a period from 0 or 0.5 |
| Stationary RF bucket | a particle started at ζ = 0 with 0.97 of the bucket height oscillates, one with 1.03 never turns back; above transition the bucket is at a sine lag of 180° and its synchrotron tune agrees with the small-amplitude formula to 2×10⁻³ |
| The Beam tab's distribution in ring coordinates | with the ring's periodic Twiss, its betatron moments match the ring (BMAG 1.000 ± 0.005, 20 000 particles); with β_x 40 % off, BMAG 1.057 = (1.4 + 1/1.4)/2 |
| The matched beam | cut at 3 rms in each normal mode and then set: its normal-mode emittances are the ones asked for to 10⁻¹², with and without RF (uncut-and-unset, a 3-rms cut leaves 0.949 of them) |
| ζ wrapped into the RF period | particles outside the bucket, slipping over several periods in 40 turns, stay in [−C/2h, C/2h); the motion is that of the unwrapped run to 10⁻⁸ |
| Checkpoint and restart | the restarted run ends bit for bit as the uninterrupted one — also when every particle was lost at the checkpoint's turn; one thread and all of them give the same bits; a restart with other probes or another statistics interval is refused |
| CUDA, numba's simulator | the device-resident run gives the CPU's particles to 10⁻¹⁵ (the simulator runs the kernels on the CPU; the tune probes are tracked on the CPU on either device); launches split over particles and turns give the same bits as one launch |
Speed (the simple Booster, no space charge, a 14-core Mac, every thread; time per turn of the tracking):
| 58 000 particles | 5 000 particles | |
|---|---|---|
| reference code B, its published settings (6th order with 6 steps in the bends and quadrupoles, 4th order with 5 steps in the sextupoles) | 2.16 s | 0.323 s |
| HELIX, curvilinear model, 6th order with 6 steps in every magnet (the sextupoles too) | 2.46 s (1.14 × the reference) | 0.207 s (0.64 × the reference) |
| HELIX, Cartesian model, expanded bend maps | 2.10 s | 0.176 s |
| HELIX, Cartesian model, exact bends (default: fourth order, about 100 kicks per bend) | 3.56 s | 0.300 s |
HELIX's time is proportional to the number of particles (60 µs per particle and turn in the default model, 14 threads); the reference's per particle rises from 37 to 65 µs between the two bunches (each element is a parallel region of its own). With 5 000 particles, 1000 turns of a whole run of the default model (statistics every turn, a tune probe: 0.320 s per turn) take 5.3 minutes, the reference 5.4 — the target is 10 minutes. With 58 000 particles they take 59 minutes (the reference 36): out of reach for both codes on this machine.
The ramp¶
The programs and the turn loop of the ramp: first the reference particle alone, then bunches through the simple Booster. The test program is a generic resonant cycle (15 Hz, 800 MeV to 6 GeV) with a voltage table rising from 0.4 to 1.2 MV.
| Check | Result |
|---|---|
| Programs | the resonant cycle meets its start momentum at t_start and its top at t_min + 1/(2f), flat before t_start; tables linear between their points with their ends held, in kinetic energy, momentum or total energy alike; programs that cannot run are refused with the reason |
| The phase rule | the reference's gain over a turn is the gain the turn needs, φs = asin(ΔE/V); above transition 180° − φs; a gain beyond the voltage capped at 90° and counted, or refused |
| Paraphasing and voltage signs | below and above transition, groups chosen by number or by name pattern: the cavities' gains add up to the gain the turn needs with the net voltage; so they do on an H⁻ deck (every voltage negative) and with one cavity's voltage reversed; groups whose voltages cancel, a cavity in two groups and cavity numbers beyond the ring's are refused, and so are cavities at another harmonic |
| Running totals | over 2000 turns the turn times add up to the reference's clock and the cavities' gains to its energy (10⁻¹³); the reference stays on the program (10⁻⁶) |
| Reference change | exact: the particle's energy, momentum and time are unchanged by the re-expression; with no gain it changes no bit |
| The map over ramped turns | over 20 turns its determinant is (p_start/p_end)³ to 2×10⁻⁶ (central differences, whose floor here is about 3×10⁻⁷) and that of its vertical block p_start/p_end to 10⁻¹⁰: the product of the cavities' momentum ratios |
| The beam follows the reference | 300 turns: the run's record of the reference is that of the reference alone bit for bit; the bunch's mean momentum deviation stays below 10⁻⁴ (6×10⁻⁶ measured) |
| Adiabatic damping | over the same 300 turns β0γ0 grows by 3.0 %: the normalized vertical emittance is conserved (3.3×10⁻⁴ measured; the geometric one shrinks as 1/β0γ0), and so is the longitudinal emittance in eV·s (3.7×10⁻⁴), both checked to 2×10⁻³ |
| Synchrotron tune | the first turn's linear RF focusing (each cavity at V_k cos φ_k) gives the accelerating bucket's small-amplitude Qs to 1.2×10⁻³ (22 cavities at the dispersion they sit at, against a smooth formula); a tracked probe's NAFF Qs lies within 2 % of the range between the bucket's at the start and at the end of the run |
| The accelerating bucket against reference code A | the reference's test of a bucket whose reference gains energy (1000 m, γt = 4, h = 40, 100 kV, 50 keV a turn, φs = 30°): Qs 0.00327284 and 0.00176524 against its 0.0032729 and 0.00176525 below and above transition; the bucket's ends −8.33333 m and 4.76991 m, mirrored above transition, as the reference's (to its 10⁻⁴ m) |
| The separatrix | 0.2 MV and 0.1 MeV a turn (φs = 30°): particles started at 0.97 of the first turn's separatrix stay in the bucket of every turn for 2.5 synchrotron periods; at 1.03 every one leaves it |
| Transition | started 150 turns below it: the phase jumps to 180° − φs on turn 150, the bunch is shortest within 50 turns of it (4.5 % shorter than at either end), no particle is lost; a ramp started above transition reports no crossing, a decelerating one the turn it crosses downwards |
| A flat ramp | the reference-change machinery with no gain tracks as the fixed-energy ring: the centroid to 10⁻¹², the second moments to 10⁻⁹ relative |
| Restart, threads, CUDA | a restarted ramp ends bit for bit as the uninterrupted one, a checkpoint is not continued with another ramp; one thread and four give the same bits; the CUDA kernel in numba's simulator gives the CPU's particles to 10⁻¹⁴ |
Against reference code B. The same program on the simple Booster as reference code B runs it (the ramp rule written again for it in a generator that never imports HELIX; the curvilinear model at the reference's settings, sixth order with six steps in every magnet, CODATA 2022), in two segments, each from 200 protons, the same in both codes: from the start of the ramp, 600 turns; and across transition, 400 turns with the jump at turn 200.
| Start of the ramp | Across transition | |
|---|---|---|
| reference γ, every turn | 4×10⁻¹⁵ relative | 8×10⁻¹⁵ |
| reference time, every turn | 1.1×10⁻¹³ relative | |
| bunch, x after 1 turn | 2.5×10⁻¹² m | |
| bunch, x at the end | 5.9×10⁻¹¹ m (600 turns) | 1.0×10⁻¹⁰ m, c t 6.3×10⁻¹⁰ m (400 turns) |
The bunch difference after a turn is the reference's own rounding, as at fixed energy; it grows with the turns as rounding differences do in a nonlinear map.
Space charge¶
The 2.5-D model of Space charge, against analytic fields, kicks and tune shifts.
| Check | Result |
|---|---|
| Field of a Gaussian beam | the solver on the exact charge of each cell against the Bassetti–Erskine field, inside 3 rms sizes, for a round beam and one three times wider than high: 3.0 and 3.4 × 10⁻³ of the largest field on 64 × 64 nodes over 8 rms sizes, 7.0 and 8.0 × 10⁻⁴ on 128 × 128 (the error falls by 4.2–4.5 each time the cell is halved: second order); the potential of the round beam, −(E1(r²/2σ²) + ln r²)/2, to 10⁻³ of its span |
| No self-force, momentum | a particle alone feels no kick (10⁻¹² of the kick its charge gives at one cell); the kicks of a bunch add up to zero transverse momentum (10⁻¹²), with the TSC and the CIC shapes; the field's Green's functions are exactly antisymmetric |
| One kick of a long Gaussian bunch | a cold round Gaussian bunch, one kick, against the analytic 2.5-D kick (J. Qiang, Phys. Rev. Accel. Beams 28, 114602 (2025), eqs. 31–32) with 400 000 particles: the transverse kick to 0.45 % rms of its largest value (1.5 % at most), the longitudinal one to 2.6 % rms; reference code B's own test of this case allows 2 % and 20 %. Run on the same case (400 000 particles), reference code B's longitudinal kicks fit a potential zero at 1 m (fitted ratio 1.007; with 6 m, 0.85) — HELIX's default reference radius |
| Round pipe | images of node charges against image charges at b²/w̄: field and potential to 4 × 10⁻¹⁶; on the wall the potential is zero to 10⁻¹²; charge beyond the wall is left out of the images and each such kick counted — in the frozen model too, whose fields are solved once per measurement (36 of 36 kicks over 3 turns; counted per field solve before: 12, or 34 with the fields solved at every kick); particles beyond the wall (60 to 100 mm from the axis of a 50 mm pipe) feel the beam's own field — the same bits as with no wall, within 2 % of a line charge's; particles 0.2 mm to 20 µm inside a 20 mm wall (0.99 to 0.999 of the radius), with the grid's nodes on both sides of the wall: a centred beam's potential against Φ_free(r) − Φ_free(b) to 2 × 10⁻⁴ of the centre's, an off-centre beam's field against its free field plus its image line charge at b²/x̄ to 0.1 % (with the images interpolated from the nodes, masked beyond the wall, as before: the potential 37 to 790 times too large, the field 2 to 10 % low) |
| Off the grid | the grid charge's multipole series (order 16) against the direct sums of its charges, at 2 to 6 times the charges' extent: to 2(r_max/r)¹⁷ of the field and the potential; the kernel's pipe series equals the host's to 10⁻¹²; beams 3 and 6 times wider than high (200 000 particles, particle-in-cell and frozen), particles just past the grid in the narrow plane and far out against Bassetti–Erskine, also for a beam 30 times wider than high (cells 30 times wider than high): within 3 % (with the series used everywhere off the grid they were wrong by 10² to 10⁸) |
| Rectangle, ellipse | a line charge off centre in a rectangle: the fitted harmonic (32 terms) against the exact series solution, 6 × 10⁻¹¹ of the free potential and 1 × 10⁻⁸ of the free field (the series' finite differences); a Gaussian beam's potential cancelled on 400 wall points to 2 × 10⁻³ of its span (rectangle and ellipse); an ellipse of equal axes is the round pipe to 10⁻⁹; the potential unit line charges at the centre and out to half the aperture leave on the wall (the check before a run, on sixteen times the fit's wall points with a rectangle's corners and an ellipse's axis ends among them — what a walk of 200 000 points around the wall finds, to 0.1 % — refused above 0.01; the fit's error is harmonic, so it bounds the potential's error inside): an ellipse of 28.5 × 21 mm with 24 harmonics 2.8 × 10⁻⁶, a rectangle of that size 6.9 × 10⁻⁴, a rectangle of 60 × 25 mm 5.5 × 10⁻³, a rectangle of 40 × 10 mm with 40 harmonics 5.6 × 10⁻³ (48: 3.5 × 10⁻³); refused: that rectangle with 24 and 32 harmonics (0.026, 0.011 — a check half a step from the corners passed 32 at 8.8 × 10⁻³), rectangles of 80 × 10 and 100 × 10 mm with any number up to 64 (0.054 and 0.095 at best) |
| Unequal weights | a beam of two populations of different weights: the frozen model's Gaussian has the charge's centroid and rms sizes to 10⁻⁹ (counted as particles its rms was 2.4 times the charge's); a particle-in-cell grid at 4 rms sizes holds the charge of a few heavy particles far out — the field beyond them within 0.3 % of the beam's (sized by the particles' count, a third of the charge fell off the grid and that field was 0.71 of the beam's); a beam at the start (5 mA) with injection (2 mA): the frozen Gaussian's rms is the circulating charge's on every turn it measures (the particles' rms is 9 % smaller); with equal weights the moments are the plain sums, bit for bit |
| Longitudinal kick | a uniform round beam (radius a) in a round pipe (radius b): the kick against −K_l f′(ζ) [½(1 − r²/a²) + ln(b/a)], 1.0014 (fitted ratio over the core); the geometry-factor model with g = 1 + 2 ln(b/a), 1.0018; none with the kick off; none for a coasting beam |
| Tune shift (L. J. Laslett, BNL-7534, 1963) | a round beam in a ring of 40 identical cells focusing both planes alike (β_x = β_y everywhere: ∮ β/(σ_x(σ_x + σ_y)) ds = C/(2ε) exactly), coasting, dQ ≈ −0.01, probes tracked 256 turns: a KV beam (100 000 particles, 32 × 32, probe at 0.5 rms sizes) 0.998 and 0.999 of the analytic shift; a Gaussian beam (100 000 particles, 64 × 64, probe at 0.2 rms sizes, the analytic shift at that amplitude) 0.976 and 0.994; the frozen model (20 000 particles, probe at 0.1 rms sizes) 0.998 and 0.991 |
| Threads, restart, CUDA | one thread and four give the same bits; a restarted run (particle-in-cell and frozen) ends bit for bit as the uninterrupted one, also through a ramp, and reports its counters; the field solve, the series and the images give the same bits over 300 repeats from arrays placed differently in memory (numpy's complex multiply, used before, did not on a 16 × 16 grid), three runs on a 16 × 16 grid with a round pipe repeat bit for bit, a fitted pipe's operator is the same for one and three BLAS threads and over 200 builds among allocations of all sizes; the frozen model's fields solved again give the same bits, its cache keeping the first kicks' fields and solving the others each time; in numba's CUDA simulator the deposit and the moments are the CPU's integers and sums for the same particles, with equal and unequal weights, the kicks within 10⁻¹⁴ |
| One particle | one particle runs as without space charge — at fixed energy and through a ramp across the edge of the RF period (ζ wrapped once a turn: a wrap after every stretch between kicks would be off by up to 3.6 µm in x in 40 turns) |
| Kicks | the kicks sit at element boundaries, the boundary nearest each of equally spaced positions: the same as the plain rule for 1 to 5000 asked, never more than the element starts, lengths adding up to the circumference; a zero-length element at the end of the ring (s = C) gets no kick of its own; the targets k C/n start at the start of the ring: a ring of 40 identical cells gets 40 and 80 kicks of equal length (targets halfway between two starts, used before, dropped one of 80 by rounding); a request denser than the starts gets every start without its targets being made (10⁶ asked on a ring of 80 starts); an element of negative length is refused |
Against reference codes on the simple Booster¶
The simple Booster at 800 MeV with its RF on (h = 84, 200 kV), the same 50 000 initial particles in every code (drawn from one seed with one recipe: Gaussian, normalized rms emittances 2.5 µm, the lattice's periodic Twiss), 1000 turns; emittances normalized rms, the betatron ones without the part correlated with the momentum. Reference code B: its 2.5-D model (64 × 64 cells over 1.1 times the beam's extent, quadratic shape, 128 line-density bins over the bunch's extent, its longitudinal kick, a kick at every slice of every element: 666 a turn); HELIX at the same element settings (curvilinear model, sixth order, six steps), 192 kicks asked (179 placed), 64 × 64 cells over the beam's extent, the potential zero at 1 m (reference code B's open boundary). Reference code A: its averaged 2.5-D model (the whole bunch's transverse distribution times its line density) with its longitudinal kick on, 96 kicks, fixed grids over ±5 rms sizes, 128 bins over the RF period, one thread (it reached fewer turns); HELIX's Cartesian model at its settings (96 kicks, 4 rms sizes).
| Check | Result |
|---|---|
| The example's bunch (reference code B's simple-Booster example: 5 × 10¹⁰, 0.1 eV s — 6π times the rms longitudinal emittance, the example's convention; σ_δ 1.1 × 10⁻³ against the bucket's half height of 2.2 × 10⁻³, so part of it lies outside the bucket and spreads over the RF period; analytic tune shift −0.137 / −0.150) | the first 10 turns agree to 2.4 × 10⁻⁴ (longitudinal) and 1.8 × 10⁻⁴ (betatron); over 1000 turns reference code B's longitudinal emittance grows by 18 % after turn 10, HELIX's not (−0.1 %), and its betatron emittances by 0.9 % and 1.1 %, HELIX's by 0.4 % and 0.4 %. Reference code A agrees with HELIX: over the 150 turns it reached, its longitudinal emittance stays flat (+0.2 % after turn 10, HELIX's Cartesian run −0.1 %), the emittances within 5 × 10⁻³ of each other, the betatron ones within 1.3 × 10⁻³ |
| A bunch matched inside the bucket (1.5 × 10¹⁰, 0.01 eV s (6π times the rms), β_t 1160.25 m — the small-amplitude match of the 6-D optics; σ_z 0.30 m, σ_δ 3.6 × 10⁻⁴, the bucket 6 σ_δ high; analytic tune shift −0.159 / −0.157; its longitudinal space charge at the bunch's centre about as strong as the RF's focusing, so that the bunch length swings between 0.30 and 0.41 m) | reference code A and HELIX's Cartesian run agree through the swings: over the 150 turns reference code A reached the longitudinal emittances to 2.2 × 10⁻³, the betatron ones to 9 × 10⁻⁴, the bunch lengths to 4.2 × 10⁻³ and the energy spreads to 2.9 × 10⁻³ (HELIX's curvilinear run gives the same longitudinal emittance). Against reference code B, two seeds each: the betatron emittances agree to 7 × 10⁻⁴ over the first 100 turns and grow over 1000 turns by 0.17–0.25 % (x) and 0.14–0.34 % (y) in HELIX, 0.40–0.44 % and 0.50–0.52 % in reference code B — the codes about as far apart as the seeds; the longitudinal emittance grows by 2.6 % in HELIX (2.0 % in the first 5 turns, as the bunch takes up its space charge, 0.6 % after) and by 91–96 % in reference code B, from the first turns on — its line density again: 5.3 % over 300 turns with 32 bins instead of 128 (48 %) |
| Convergence (HELIX, the bunch inside the bucket, 200 turns) | the longitudinal emittance changes by 3.5 % — 3.8, 3.5 and 3.5 % with 25 000, 50 000 and 100 000 particles, 3.5 % on 128 × 128 cells, 3.2 % with 128 line-density bins, 3.5 % with 96 kicks: the bunch taking up its space charge (its rms emittance oscillates: 2.6 % at turn 1000), as in reference code A; without space charge it holds to 0.1 % (reference code B, 650 turns). The betatron emittances grow by a few hundredths of a percent, less with more particles (x: 0.090, 0.061, 0.052 %) — the heating by the sampling noise —, the same on 128 × 128 cells, more in y with 96 kicks (0.12 % against 0.05 %) |
| Why reference code B's longitudinal emittance grows | its line density — binned over the beam's extent at every kick, differenced forward with no filter, read with no interpolation — is noisy: without its longitudinal kick its longitudinal emittance stays flat (as without space charge: 0.00543 eV s rms); with 32 bins instead of 128 it grows by 0.25 % in 300 turns instead of 2.4 %; HELIX given as noisy a line density (512 bins, no filter) grows by 1.6 % (with 64 bins and a one-bin filter: −0.2 %). HELIX's runs in these comparisons fix their bins: 256 against reference code B, 128 against reference code A. With the longitudinal kick off in both codes the longitudinal emittances agree to 6 × 10⁻⁵ over 300 turns and the betatron ones to 3 × 10⁻³ |
Injection¶
| Check | Result |
|---|---|
| H⁻ stripping | an 800 MeV H⁻ leaves the foil as a 799.1296 MeV proton (the same velocity: its momentum m_p/m_H⁻ = 1 − 1.0880 × 10⁻³ of the ion's, δ = −7.067 × 10⁻⁴ against an 800 MeV proton), exactly; the injected centroid placed exactly |
| Bookkeeping | injected = circulating + lost + waste, in macroparticles and real particles, exactly; each macroparticle I T / (e B n) real particles; a beam beside the foil all waste (missed), the two-step stripping's waste within 3 σ of its fractions |
| A Gaussian source's cut | at least 1 σ — a smaller one is refused (a pair of normal numbers is kept with probability 1 − exp(−c²/2): at 10⁻⁴ σ the draw never ended); at 1 and 3 σ the numbers drawn are those drawn before, bit for bit |
| The foil's statistics | 10⁶ crossings of 6 mg/cm² carbon at 800 MeV: the rms angle kick is Highland's θ0 and the mean loss the foil element's to 0.5 % |
| Repeated crossings of a foil element | a matched beam crossing a foil element every turn: after N crossings its vertical action has grown by β_y θ0(N x)²/2, Highland at the total thickness (each crossing takes the growth of θ0² from one thickness to the next), to 5 % over 20 crossings of 6 mg/cm² — 1.37 times the sum of 20 one-crossing kicks |
| A bump of exact dipoles | four thick kickers (20 mrad, 0.3 m) compiled as straight-frame dipoles around a foil in the simple Booster: the bump's height at the foil is the arcs' exact one, 2 L θ/(1 + p_z) + (g − L) θ/p_z, to 10⁻¹⁴ m, in the optics and in the run's closed-orbit record turn by turn as the waveform takes it down |
| Painting | the bump collapsing as 1 − √(t/T) under a beam injected at its full height: the particles' actions are those the closed orbit's distance gives each turn (uniform in action), to 1 % |
| Restart, threads, devices | a run restarted in the middle of the injection ends bit for bit as the run that was not interrupted; the same bits for 1 and 3 threads; the CUDA path (simulator) as the CPU to 10⁻¹⁴ — its foil step finds the particles inside the foil on the GPU and moves only those to the host and back, and the injected ones to the GPU (the foil's record the CPU's exactly, with circulating particles crossing it turn after turn) |
| The bump dipole against reference code A | a thick kicker as a bump dipole — horizontal and vertical, with and without fringe, with face angles, both signs — tracks 48 particles crossing it at up to 70 mrad as reference code A's straight-frame dipole with dipole-only faces, to 3 × 10⁻¹⁵ m (the reference's position update (p_z′ − p_z)/k0 rounds at about 10⁻¹⁶/k0) |
| The foil against reference code C | a pencil beam of 800 MeV protons through carbon of 380 µg/cm² to 20 mg/cm²: the core of reference code C's multiple scattering (single screened-Rutherford scatters; a Gaussian fitted to the central 98 % of the angles, the definition Highland's formula is fitted to) is HELIX's θ0 (Lynch & Dahl's form) to 1.1 % at 6 and 20 mg/cm² (0.8–1.1 % and 0.3 % narrower) and to 6 % at 380 and 600 µg/cm² (6.0 % and 4.7 % narrower), where the formula is extrapolated two orders below the thicknesses it was fitted to (x/X0 from 10⁻³). At thin foils the distribution is not a Gaussian: its central 68 % is 17–20 % narrower than θ0 and its single-scattering tail — not modelled in HELIX — makes its full rms 1.42–1.59 θ0. Reference code C's full model takes Rossi's width (no logarithmic term: 1.75 θ0 at 600 µg/cm²) and loses its mean ionization loss with no straggling: HELIX's Bethe mean loss is 2.7 % below it, the default tabulated rate 18.8 % below (a floor below about 2 GeV) |
| Multi-turn injection into the CERN PS Booster against reference code A (a developers' check: the lattice is not distributed with HELIX) | ring 1 at 160 MeV (flattened by MAD-X), the measured injected distribution seven particles a turn for 142 turns, the four-magnet injection chicane collapsing over 500 turns, no RF; both codes thick, at the same bend integration, the chicane magnets as straight-frame dipoles with dipole-only faces: the closed orbit at the injection point agrees to 9 × 10⁻¹² m over the bump from 46 mm to 0; with the chicane at full strength the tunes agree to 6 × 10⁻¹⁰ (Q_y 4.4706 against 4.45 without: the chicane magnets' edge focusing); turn by turn the circulating beam's centroid within 0.3 % (x) and 1.1 % (y) of its rms size and the rms sizes within 0.35 % and 1 % |
Losses, space charge and painting over many turns¶
Whole stretches of a cycle, each code from the same particles.
| Check | Result |
|---|---|
| Losses through transition against reference code B | the generic program of the ramp check from 300 turns before transition for 600 turns, a rectangular aperture (19 mm × 4 mm half widths) at the start of every cell, where the dispersion is 3.19 m; 1000 protons in a bunch not matched to the bucket, whose momentum spread swings to 2.7 × 10⁻³ within 100 turns (reference code B tracked cell by cell, so that each loss has its cell): the same 584 particles are lost, each on the same turn at the same cell (483 before transition, 101 after; 567 beyond the aperture's horizontal half width, from turn 10 on, and 17 beyond its vertical one, all in the first turn — 7 of them outside it from the start); the 416 survivors end within 1.6 × 10⁻¹⁰ m in x, 3.7 × 10⁻¹¹ m in y and 1.2 × 10⁻¹⁰ m in ζ of the reference's, the reference particle's γ within 10⁻¹⁴ |
| Space charge through transition against reference code B | the generic program of the ramp check from 300 turns before transition for 600 turns; 5000 protons standing for 1.5 × 10¹⁰ in a bunch 0.1 m long (rms) with about a third of the momentum spread its bucket matches, so that it rotates to 3.4 cm at turn 146 while its momentum spread swells; the reference code's 2.5-D space charge with its transverse kick alone (its longitudinal kick heats a bunch by itself, above). Without space charge the codes' particles end within 1.3 × 10⁻⁵ m of each other; with it the horizontal betatron emittance grows by more than 40 % before transition. With the line density over 4096 bins of the RF period (1.4 mm), HELIX agrees turn by turn to 2.7 % in the horizontal betatron emittance (+44 % to transition against +47 %), 0.64 % in the vertical, 0.15 % in the bunch length, 0.30 % in the energy spread and 1.8 % in the longitudinal emittance, which both codes see end 16 times as large. The growth depends on how finely the line density is resolved — with 256 bins (2.2 cm) HELIX's peak line density is 20 % low at the shortest bunch and the growth +13 % — and on its noise: at 5000 particles and fixed bins, four samples of the bunch grow by 39.2–42.9 % (mean 41.2 %) with 2048 bins and by 44.7–49.9 % (48.5 %) with 4096, while the line density the particles see differs by less than 0.3 % between the two; with 20 000 particles the reference gives +42 % and HELIX +41 % with 4096 bins, the two within 1 % turn by turn. HELIX's default line density ("auto": its bins follow the bunch, its smoothing the number of particles — The line density) gives +36.2 % at 5000 particles (four samples, 34.4–37.9 %) and +36.1 % at 20 000: the same with four times the particles, where the reference's growth falls from +47 % to +42 % |
| Painting through a foil against reference code C | a generic FODO ring (12 cells, 201.6 m, 800 MeV, no RF: a coasting beam; the tunes of both codes agree to 10⁻¹⁰), its quadrupoles without reference code C's hard-edge fringe (HELIX's have none) and HELIX's drifts expanded as the reference's; a vertical bump of four thin kickers (6 mm at the foil) falling as 1 − √(t/T) over 300 turns, held through each turn at its value at the turn's start; 20 protons a turn injected at 6 mm above the final orbit — the same 6000 in both codes, in the same order and at the same energies — through a foil 2–20 mm in y. With a foil that only counts its hits: the same 5996 particles circulate (four miss the foil in both); the foil is hit 604 089 times against 604 109 (at most 7 hits apart in a turn: particles at its edge counted a turn earlier or later), 100.75 times per particle on average; the emittances agree turn by turn to 3 × 10⁻⁶ (y) and 1.7 × 10⁻⁵ (x), each particle's final vertical action to 6 × 10⁻⁷ at the median; the particles end 5 µm apart at the median in y (2 µm in x) |
| Repeated crossings of a foil | the same with 600 µg/cm² of carbon: each circulating particle crosses it about 100 times. HELIX gives a particle's k-th crossing the growth of Highland's θ0² from k − 1 to k thicknesses, so that its crossings add up to Highland at their total thickness (θ0² is not additive over layers: kicks of the one-crossing width would be 1.30 times too narrow in angle after the median particle's 99 crossings). Reference code C scatters by single collisions: the painted vertical action agrees to within 1.4 % at the median over five of its runs (its foil's random numbers are not seeded), and its foil widens the median betatron action in x — the plane the painting leaves alone — 0.84–0.90 times as much as HELIX's. A single-scatter estimate (screened Rutherford, carbon, 800 MeV) gives 0.88: Highland at the total thickness (Lynch & Dahl's form; x/X0 = 1.4 × 10⁻³, inside the range it was fitted over, to about 11 %) is 6 % wider than the core of 100 crossings' single scatters, and with the crossings at different betatron phases that core grows 5 % less than in one plane. At the 99th percentile the reference's foil share is 1.8–2.5 times HELIX's: its Rutherford tail, which HELIX does not model. (HELIX's foil also slows the particles — a mean δ of −8.6 × 10⁻⁵ after 100 crossings — and the reference's scattering-only foil does not: the actions are taken with the dispersion out.) |
References¶
- J. Laskar, "Frequency analysis for multi-dimensional systems. Global dynamics and diffusion", Physica D 67, 257 (1993) — NAFF.
- H. Yoshida, "Construction of higher order symplectic integrators", Phys. Lett. A 150, 262 (1990).
- M. Bassetti and G. A. Erskine, "Closed expression for the electrical field of a two-dimensional Gaussian charge", CERN-ISR-TH/80-06 (1980).
- J. Qiang, S. Lidia, R. D. Ryne and C. Limborg-Deprey, "Three-dimensional quasistatic model for high brightness beam dynamics simulation", Phys. Rev. ST Accel. Beams 9, 044204 (2006); erratum 10, 129901 (2007) — the integrated Green's function.
- J. Qiang, Phys. Rev. Accel. Beams 28, 114602 (2025), doi:10.1103/npx1-nclt, eqs. 31–32 — the analytic 2.5-D Gaussian kick.
- L. J. Laslett, "On intensity limitations imposed by transverse space-charge effects in circular particle accelerators", BNL-7534, p. 324 (1963).
- F. W. Jones, "A method for incorporating image forces in multiparticle tracking with space charge", Proc. EPAC 2000, p. 1381.
- R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles (1988).
- V. L. Highland, "Some practical remarks on multiple scattering", Nucl. Instrum. Methods 129, 497 (1975); G. R. Lynch and O. I. Dahl, "Approximations to multiple Coulomb scattering", Nucl. Instrum. Methods B 58, 6 (1991) — the foil's θ0.
- M. S. Gulley et al., "Measurement of H⁻, H⁰, and H⁺ yields produced by foil stripping of 800-MeV H⁻ ions", Phys. Rev. A 53, 3201 (1996).