34 Linear Transfer Maps

OPTION, ENABLELINEARTRANSFERMAPS=TRUE enables numerical linear maps during the design-reference OrbitThreader pass. This describes the implementation on the 546-add-linear-transfer-matrix-to-each-elemen development branch. No analytic element map is used by the calculator: every matrix is reconstructed from external-field tracking, including fields supplied by supported field maps. Collective fields are excluded; RF cavities and traveling-wave structures are currently rejected. The requested LINE interval, or one RING turn, is threaded independently of the production MAXSTEPS budget.

34.1 Coordinates and moving frame

The reported coordinates are

\[ X=(x,x',y,y',\zeta,\delta)^T,\qquad x'=\frac{u_x}{u_s},\quad y'=\frac{u_y}{u_s},\quad \zeta=-\beta_0c(t-t_0),\quad \delta=\frac{|\mathbf u|}{u_0}-1. \]

Here \(\mathbf u=\mathbf p/(mc)=\vec\beta\gamma\) is OPALX’s dimensionless mechanical momentum, \(u_0=|\mathbf u_0|\), and \(\beta_0=u_0/\sqrt{1+u_0^2}\). Positions and \(\zeta\) are in metres; slopes and \(\delta\) are dimensionless. Ordinary entrance/exit planes are normal to the corresponding reference tangent \(\mathbf e_s=\mathbf u_0/u_0\). At a RING return, the starting section and axes are reused even when the returned momentum is not parallel to its normal. Each ray is recorded at its own crossing time, not at a shared laboratory time.

For private map rays, the last nominal step brackets the exit-plane crossing. Boris retains tracked bisection because changing trial durations produced measurable map differences in the stable DBA benchmark. For RK4 and DOP853, a safeguarded secant search proposes trial arrival times from the signed plane distances; a trial that fails to halve the bracket forces bisection on the next iteration. Each trial is fully integrated from the original bracket start, including field-support subdivision: coordinates are never interpolated. The search retains the absolute time-bracket tolerance \(10^{-12}|\mathrm{DT}|\) (seconds), with early termination for an exact plane hit or adjacent representable times. Nonfinite residuals or failure to converge are errors. This changes trial times and final roundoff relative to fixed bisection, without changing integration methods, support tolerances, or the crossing direction and return-time gate.

The transverse frame follows minimum-rotation (Bishop) transport. Between two reference tangents, rotate all three axes about \(\mathbf a\propto\mathbf e_s\times\mathbf e_s'\) through

\[ \theta=\operatorname{atan2} \left(|\mathbf e_s\times\mathbf e_s'|,\mathbf e_s\cdot\mathbf e_s'\right). \]

There is no added roll about the tangent. Straight sections therefore retain their transverse axes without requiring a nonzero curvature.

34.2 Finite differences and overlaps

For each coordinate \(j\), launch two rays with \(X_j=\pm\epsilon_j\) and all other coordinates zero. LINEARTRANSFERMAPSTEPS sets the six starting amplitudes; the defaults are \(10^{-3}\) in each coordinate’s units. Construct the six columns

\[ D_{\mathrm{in},j}=\frac{X^+_{\mathrm{in},j}-X^-_{\mathrm{in},j}}{2},\qquad D_{\mathrm{out},j}=\frac{X^+_{\mathrm{out},j}-X^-_{\mathrm{out},j}}{2},\qquad M D_{\mathrm{in}}=D_{\mathrm{out}}. \]

A pivoted inversion of the input differences gives the matrix. Centered differences have an \(O(\epsilon_j^2)\) truncation error in a differentiable map. Reducing DT alone cannot remove this error floor.

34.2.2 Segment ownership

The reference path is partitioned wherever its nominal body-owner set changes, independently of field-support boundaries. Each segment has a unique identifier and is attached to every nominal owner. Ordinarily an element owns its entrance-to-exit map; overlapping nominal bodies or a ring-seam split can produce several segment maps. activeElements records the union of support-selected elements encountered by the reference in that interval, not the owner set. Private rays always evaluate the fields at their own positions:

\[ \mathbf E(\mathbf r,t)=\sum_{i\in A(\mathbf r)}\mathbf E_i(\mathbf r,t),\qquad \mathbf B(\mathbf r,t)=\sum_{i\in A(\mathbf r)}\mathbf B_i(\mathbf r,t). \]

The reference segment’s owners do not restrict this ray-dependent field set. A nominal drift map includes neighbouring fringe fields when present. Shared segments are composed once, and intervals without a nominal owner are included even when they contain a field tail:

\[ M_{\mathrm{total}}=M_N\cdots M_2M_1. \]

The combined matrix and diagnostics are printed to stdout. For a periodic reference pass this is a same-section return map, not proof of a closed orbit. Element attachments can be traversed with OpalBeamline::getLinearTransferMapsInReferenceOrder(); that traversal deduplicates shared segments but does not contain unowned intervals. The complete product is available from OrbitThreader’s combined-map getter.

34.3 Geometry, fringe support and ring return

Three quantities must remain separate:

Design circumference: \(C_{\mathrm{design}}=\sum_i L_i\), counting each nominal occurrence once. Bend lengths are design arc lengths (including the chord-to-arc conversion for an RBEND); drifts fill the nominal gaps.
Field support: the spatial domain in which a source contributes to \(\mathbf E\) and \(\mathbf B\). This can extend into a drift, another magnet, or across the ring seam, without adding geometric length.
Tracked return length: the distance actually travelled to the starting transverse plane. It is a closed-orbit circumference only after position and momentum closure have been established.

Use explicit 6D entrance poses for new lattices; ELEMEDGE is deprecated and is not part of this design. RING is a beam-sequence declaration, not an element, and has no CIRCUMFERENCE input attribute. Its computed length does not check that the independently supplied 3D poses form a contiguous, closed geometry.

The ring reference search detects a return to

\[ (\mathbf r-\mathbf r_{\mathrm{start}})\cdot\mathbf n_{\mathrm{start}}=0, \qquad \mathbf n_{\mathrm{start}}=\mathbf u_{\mathrm{start}}/|\mathbf u_{\mathrm{start}}|, \]

in the launch crossing direction, after travelling more than \(C_{\mathrm{design}}/2\). The crossing is refined by reintegration and bisection. Failure to return within \(2C_{\mathrm{design}}\) (or four nominal flight times) raises an error. This bounded search targets a simple one-circuit ring; it is not a general multi-loop return detector or a closed-orbit solver. The final map plane and axes are the starting plane and axes, rather than a new plane normal to a possibly mismatched returned momentum.

Stdout reports design circumference, reference return length, position closure residual in metres, and relative momentum mismatch. The distance-keyed IndexMap wraps by the measured reference return length, not the design circumference. Periodic reuse still assumes a sufficiently repeatable reference orbit; a non-closed launch is not made periodic by this bookkeeping. Production TURNS stops on directed crossings of the fixed launch plane, not on nominal design length. When no explicit return-count or energy stop is active, ring progress is reported relative to the design circumference.

Changing a numerical field cutoff must not move element poses or nominal map boundaries. At fixed geometry the design circumference must remain invariant, while tracked orbits and matrices should converge as the cutoff is enlarged. Do not add an equivalent thin-edge kick on top of a distributed representation of the same fringe effect. The present change does not alter the native SBEND Enge profile or its existing FINT coefficient convention.

The regression coverage includes an analytic circular orbit whose travelled circumference deliberately differs from the supplied design-length guard, a nominal drift traversed by a neighbouring field, and native Enge support cutoffs of two, three, four and five full gaps at fixed body geometry. The latter compares all 36 matrix entries: the three- and four-gap results agree with the five-gap result within \(10^{-8}\) in the test’s coordinate units, while the two-gap cutoff produces a resolved change above \(10^{-6}\). Nominal lengths and body endpoints remain unchanged. These are case-specific assertions, not universal cutoff recommendations. The relevant C++ tests pass, including the OrbitThreader suite on two ranks. The six-cell DBA smoke run separately records its non-closed reference return.

34.3.1 Documentation comparison: MAD-X, Bmad and MaryLie

This is a documentation comparison, not an executed numerical code-to-code benchmark. The conclusions below concern distinct models, not interchangeable predictions for an arbitrary finite fringe field.

Code and primary source	Geometry and fringe treatment	Relevant cross-check
MAD-X element definitions, SBEND and DIPEDGE	SBEND length is the design arc; DIPEDGE is a zero-length edge correction. Fringe focusing does not add circumference.	Hard-edge and convention-matched thin-edge maps; not a distributed-overlap oracle.
Bmad manual, 4 August 2026, §§5.13, 5.16, 5.18	Element length is separate from field-map extent. `field_overlaps` associates a source extending past its ends with other elements.	Closest architectural analogue; use a supported integration method with the same distributed field. `bmad_standard` tracking ignores field maps.
MaryLie 3.0 manual, December 2003, §§6.7, 6.31	Separate entrance/body/exit maps; `arc` supplies geometry bookkeeping for imported maps. The manual warns about hard-edge fringe approximations and recommends GENMAP for realistic treatment.	Independent map-formalism and edge-model comparison; these sections do not establish automatic distributed-field overlap handling.

The CERN CAS treatment likewise factors a finite fringe transport into body/drift transport and an edge correction referred to a fixed boundary; a fringe correction is not an extra lattice interval. See the CAS numerical tracking lecture discussion of fringe fields. The code-specific definitions above remain the authority for an actual numerical comparison.

Before comparing full \(6\times6\) matrices, match entrance/exit planes, reference rigidity, magnetic field integral, edge angles, and phase-space coordinates. OPALX’s \(\boldsymbol\beta\gamma\) momentum and reported slopes must be converted consistently to another code’s canonical momentum and longitudinal conventions. In particular, MAD-X defines \(g=2\,\mathrm{HGAP}\) and

\[ \mathrm{FINT}=\int_{-\infty}^{\infty} \frac{B_y(s)[B_0-B_y(s)]}{gB_0^2}\,ds, \]

whereas MaryLie frng takes the full pole-to-pole gap directly. Equal parameter names or numbers do not guarantee equal fringe fields. The intended numerical validation sequence is a shared hard-edge DBA baseline, a convention-matched thin-edge comparison, and a shared distributed-field OPALX/Bmad overlap and ring-seam comparison. Those external-code runs remain future work.

34.4 Field-support boundary subdivision

Both the map-enabled reference pass and its private rays use the selected numerical integrator. LINEARTRANSFERMAPINTEGRATOR="BORIS" (default, alias "LF2") uses drift–Boris-kick–drift, evaluating fields at the drift midpoint. "RK4" and "DOP853" select the fixed-step Runge–Kutta methods described below. The selection is independent of Richardson levels. It does not change the production-particle integrator or select RK methods for map-disabled or secondary-species threading. RING reference threading always resolves support boundaries, including with Boris and with map calculation disabled. Its IndexMap records accepted subintervals; adjacent intervals with identical support sets are coalesced. The support set is obtained from each element’s isInside() query, using that element’s placement and, for a sector bend, its curved containment chart. Support can extend beyond the mechanical body because of fringe fields.

A Boris trial step is accepted if support sets at the start, midpoint and end agree. Otherwise its duration is halved; after advancing the first half, the second half is recomputed from the resulting state. This continues until a boundary straddle satisfies

\[ c|h|\leq\max\left(10^{-12}c|\Delta t|, 64\epsilon_{\mathrm{mach}}\max(1\,\mathrm m,|\mathbf r|)\right). \]

Here \(h\) is the trial duration and \(\Delta t\) the requested ray step. The floating-point term prevents endless subdivision when position updates become unresolvable. Trial steps are also capped at \(L_{\min}/(4c)\), where \(L_{\min}\) is the shortest positive longitudinal support extent or nominal body arc length. The body-length cap resolves map ownership when support is wider than the body. Thus a regular longitudinal crossing cannot jump over a whole thin field region without a support sample. The existing input time-step safety check remains active.

Smooth regions retain the selected integrator’s order. The boundary operation is numerical subdivision, not an exact element map, and it handles each ray’s crossings independently. Grazing crossings, transverse holes smaller than the sampling scale, and discontinuities hidden inside a field table still require adequate nominal resolution. For RK methods, all field-evaluation stages, including the independently integrated diagnostic midpoint, must have the same support set as the start and end; otherwise the trial is subdivided using the same tolerance. This also applies to map-enabled RK pre-roll. This algorithm does not estimate interpolation error or guarantee detection of arbitrary geometry at an arbitrarily large DT.

During one ray advance, an already evaluated support set can be reused at an identical recursive start or when an accepted endpoint becomes the next half-step’s start. This assumes containment is a fixed function of position during that advance. The reuse retains all original support predicates, field-evaluation order, accepted steps and tolerances; no field values or integrated trajectory prefixes are cached. Rejected trial endpoints are discarded, and no membership cache survives into a separate advance call or another ray. This host-only optimization adds no MPI communication, OpenMP scheduling, or device transfers.

34.5 Map-only Runge–Kutta integration

RK4 and DOP853 integrate the relativistic Lorentz equations and signed path length:

\[ \begin{aligned} \dot{\mathbf r}&=c\mathbf u/\gamma, & \dot{\mathbf u}&=\frac{qc}{\mathcal E_0} \left(\mathbf E+\frac{c}{\gamma}\mathbf u\times\mathbf B\right),\\ \dot s&=c|\mathbf u|/\gamma, & \gamma&=\sqrt{1+|\mathbf u|^2}. \end{aligned} \]

Here \(\mathbf u=\boldsymbol\beta\gamma\), \(\mathbf r\) and \(s\) are in metres, \(t\) in seconds, \(\mathbf E\) in V/m, and \(\mathbf B\) in tesla. \(\mathcal E_0=mc^2\) is the numerical rest energy in eV and \(q\) is the signed charge in elementary-charge units, matching OPALX’s Boris convention.

For \(\mathbf y=(\mathbf r,\mathbf u,s)\) and signed duration \(h\):

\[ k_i=h f\!\left(t_0+c_i h,\mathbf y_0+\sum_{j<i}a_{ij}k_j\right), \qquad \mathbf y_1=\mathbf y_0+\sum_i b_i k_i. \]

Every stage evaluates the summed external fields at its own position, momentum and time. Classical RK4 uses four stages and has smooth global order four. DOP853 uses the twelve-stage, eighth-order integration formula from Hairer and Wanner’s implementation, with coefficients recorded in Algorithms/RungeKuttaTableau.h.

Both selections use fixed nominal DT with support-boundary subdivision. DOP853’s embedded fifth-/third-order estimators, adaptive error controller, stiffness detection and dense output are not implemented here. This makes DT directly comparable across the convergence study; selecting DOP853 does not request an automatically controlled error tolerance.

An independent half-duration solve supplies a same-order diagnostic midpoint; the accepted endpoint still comes from the full-duration solve, not two half steps. Including the midpoint field sample, a smooth trial costs 9 field evaluations for RK4 and 25 for DOP853, versus 1 for Boris. Subdivision and endpoint localization add further evaluations. Weighted stage increments use compensated sums; position, time and path retain the compensation described below. Momentum updates and field evaluations remain double precision.

These explicit RK methods are not symplectic, are not exactly time-reversible, and do not preserve magnetic energy exactly. Their orders require sufficiently smooth fields and sufficiently resolved steps; field-table interpolation, event localization and finite-difference noise can dominate map errors. The host-only implementation adds no MPI reductions, OpenMP scheduling changes, GPU kernels or host/device transfers. Each rank computes its own identical map.

34.6 Position, clock and path bookkeeping

The host reference/map tracker performs each Boris half drift directly in metres, without repeatedly scaling absolute positions by \(c\Delta t\) and back:

\[ \mathbf r_{1/2}=\mathbf r_0+\frac{ch}{2}\frac{\mathbf u_0}{\gamma_0}, \qquad \mathbf r_1=\mathbf r_{1/2}+\frac{ch}{2}\frac{\mathbf u_1}{\gamma_1}, \qquad \gamma_i=\sqrt{1+|\mathbf u_i|^2}. \]

Here \(\mathbf u=\mathbf p/(mc)=\boldsymbol\beta\gamma\) and \(h\) is a signed step duration in seconds. The Boris kick between the half drifts is unchanged. Signed path length is integrated from speed:

\[ \frac{ds}{dt}=c\frac{|\mathbf u|}{\sqrt{1+|\mathbf u|^2}}, \qquad \Delta s=\frac{ch}{2} \left(\frac{|\mathbf u_0|}{\gamma_0}+\frac{|\mathbf u_1|}{\gamma_1}\right). \]

Unlike the endpoint chord, this does not shorten the accumulated path just because the velocity turns. It is exact for constant speed, apart from momentum and floating-point errors, and second-order quadrature under acceleration. It does not increase the integration order of Boris.

Position, time and path additions use Kahan compensation. Each state retains a high part \(a\) and correction \(e\), representing \(a-e\). Both parts survive support subdivision, reference sampling and shadow-ray initialization. Coordinate and time differences use

\[ (a-b)-(e_a-e_b), \]

so small flight-time differences are not discarded when subtracting large absolute clocks. Field evaluators still receive ordinary double-precision positions and absolute time; this is not extended-precision field interpolation.

The requested start/stop path is located by bisection of numerically tracked states, including field-support subdivision, rather than converting a path fraction directly into a time fraction. The time tolerance uses the support localization formula above with the enclosing interval as \(\Delta t\); the endpoint is labelled by the requested path within that tolerance. Element boundary refinement likewise stores the path integrated by the trial ray.

These changes affect host OrbitThreader/reference and private map rays, including reference bookkeeping with map calculation disabled. They do not modify the production-particle pusher, MPI reductions, OpenMP scheduling or GPU kernels. Floating-point trajectories and maps can therefore change relative to the earlier implementation; the DBA comparison must use the same input, grid and analytic oracle. Symplecticity is not guaranteed by compensated accumulation.

34.7 Diagnostics and validation

The calculator reports both

\[ |\det M-1|,\qquad \max_{ij}|(M^TJM-J)_{ij}|,\qquad J=\operatorname{diag}(J_2,J_2,J_2),\quad J_2=\begin{pmatrix}0&1\\-1&0\end{pmatrix}. \]

Unit determinant is necessary but not sufficient for symplecticity in canonical coordinates. Normalizing momentum to \(\beta\gamma\) does not by itself make mechanical momentum or slopes canonical. These reported coordinates therefore make the canonical-\(J\) test a diagnostic, not a universal invariant. Field interpolation, finite differences and state-dependent subdivision also preclude an unconditional symplecticity claim.

The source repository’s sandbox/map-2 cases compare all 36 entries against closed-form drift, quadrupole, FODO and DBA matrices. The DBA convergence script records stdout with --info 2, signed \(R_{16}\) and \(R_{26}\), determinant and canonical-\(J\) residuals, and the difference from the finest numerical map. That last comparison helps distinguish time discretization from the fixed finite-difference floor.

34.8 Fixed-section return-map service (experimental)

The OneTurnMap C++ service is the first building block for the general RING closed-orbit finder. The experimental COF input command now uses it; it does not replace the segmented transfer-map calculator described above.

Its state is \(v=(x,p_x,y,p_y)\) on a fixed local \(z=0\) section. Positions are in metres and momenta are mechanical \(p/(mc)\), not slopes. Given fixed total normalized momentum \(p_0>0\), each launch uses \(p_z=\sqrt{p_0^2-p_x^2-p_y^2}>0\). The same section and momentum constraint apply to every finite-difference perturbation.

The shared ExternalFieldRayTracker advances each ray independently. A negative section-distance excursion exceeding the default 1 nm hysteresis arms the next positive crossing; neither the launch nor the reverse crossing counts. Accepted support-resolved substeps are examined, and the return is located by tracked bisection to a time-bracket width of \(10^{-12}DT\) or representable resolution. Coordinates are never teleported or projected onto the plane. DT must still resolve the section excursions; a self-intersecting orbit may need a more selective section than the first directed return. Path and step budgets bound failed searches. Losses propagate as errors; field-free intervals are allowed.

The derivative uses eight returns,

\[ M_{:,j}=\frac{T(v+h_j e_j)-T(v-h_j e_j)}{2h_j}, \]

with separately specified positive position/momentum steps. It is not obtained by extracting a block of the slope-coordinate 6D matrix. Returned momenta are not renormalized: integration energy drift remains a separate diagnostic. Topology, geometric closure and static-magnetic physics validation belong to the caller (CofCmd for the input command). This host-only service constructs no bunch or collective solver and does not alter existing tracking or parallel execution.

34.8.1 Damped fixed-point solver (C++ interface)

ClosedOrbitSolver accepts a deterministic four-coordinate return map, including OneTurnMap; CofCmd provides its input-file adapter. For \(F(v)=T(v)-v\), central differences form \(M=dT/dv\) and each Newton correction solves

\[ \left[S^{-1}(M-I)S\right]w=-S^{-1}F,\qquad \Delta v=Sw. \]

\(S\) is diagonal with positive, user-supplied characteristic coordinate scales (metres or normalized momentum); its default entries are one. A small 4-by-4 column-pivoted Householder QR implementation recomputes trailing column norms and solves the triangular system without forming an inverse. The ratio of the smallest to largest QR diagonal magnitude is reported as a rank diagnostic, not a condition number. A ratio at or below the default \(10^{-12}\) rejects the Newton system. This threshold is configurable and is not a bound on errors in the numerically differentiated map.

Backtracking tries \(\lambda=1,1/2,\ldots\) and requires

\[ \|S^{-1}F(v+\lambda\Delta v)\|_2 \leq(1-10^{-4}\lambda)\|S^{-1}F(v)\|_2. \]

The default permits 20 halvings and 20 accepted iterations. With damping disabled, only the full step is tried and must still reduce the residual. Lost line-search trials are rejected without changing the last accepted state. A failed initial evaluation or derivative probe aborts with a diagnostic; finite-difference steps are not silently reduced after losses.

Both the residual and the undamped Newton correction must meet separate position and normalized-momentum tolerances (both numerically \(10^{-10}\) by default). A fresh map evaluation must then confirm residual closure. Singular systems are rejected even if the supplied point already has zero residual, since the solver cannot establish a locally isolated solution. Defaults for central-difference steps are \(10^{-6}\) in each coordinate’s units; these are starting values, not validated settings for a particular accelerator.

The result contains status, last accepted coordinates, residual, final Jacobian on success, evaluation count, and accepted-iteration diagnostics. Eigenanalysis is a separate C++ operation described below. Analytic coupled affine maps and nonlinear fixed-point tests validate the solver; the first native cyclotron benchmark is described below. A stable hard-edge DBA-ring benchmark is described next; it does not validate finite-fringe DBA rings or arbitrary ring optics. Timestep, finite-difference and launch-amplitude convergence remain independent requirements. No production pusher, reduction order, MPI collective, GPU kernel, or OpenMP scheduling is changed by this serial algorithm.

34.8.2 Stable non-achromatic hard-edge DBA ring

The sandbox/Regression-Tests/dba-ring-cof benchmark repeats six cells of two 30-degree, 2 m-radius sector bends, two 1 m drifts and a 0.2 m quadrupole. Changing the quadrupole to \(K_1=+2.84\,\mathrm{m}^{-2}\) in OPALX’s electron convention gives stable, non-achromatic optics; circumference is 25.76637061436 m. The independent hard-edge reference is \(M_{\mathrm{cell}}^{6}\). For the on-axis fixed-momentum orbit its transverse slope matrix is converted to mechanical coordinates by \(M_u=T M_{\mathrm{slope}}T^{-1}\), with \(T=\operatorname{diag}(1,u_0,1,u_0)\) and \(u_0=p/(mc)\) from the actual BEAM.

DOP853 tracking at \(\Delta t=5\times10^{-12}\) s and central differences \(h_x=h_y=3\times10^{-5}\) m, \(h_{u_x}=h_{u_y}=u_0(3\times10^{-5})\), gives:

Principal fractional mode	Analytic	Tracked
Horizontal	0.229892430867	0.229892431493
Vertical	0.310463602980	0.310463603175

The maximum slope-matrix entry error, with positions scaled by 1 m, is \(6.32\times10^{-9}\); maximum eigenvalue-modulus error is \(5.13\times10^{-10}\). The unchanged \(10^{-8}\) unit-circle tolerance classifies both modes as stable. The twice-larger timestep also agrees, and a seed displaced by 0.1 mm in both planes converges to the design orbit. Smaller finite-difference steps can amplify tracking noise enough to fail the strict unit-circle test; the saved scan retains those outcomes. Integer tune branches, spectral tunes, off-momentum optics and distributed fringe fields are not validated here.

For conventional elements, COF assigns aperture checks to the nearest finite design centreline, independently of aperture dimensions. This prevents a remote ring arm’s infinite longitudinal slab from causing a false loss. Circular SBEND arcs and straight/RBEND axes include endpoints; ties within floating-point roundoff are all checked. The assignment is stateless for independent shadow rays and return-root reintegration. It assumes a local orbit neighbourhood of a distinct path, not a solid model for intersecting beam pipes. Cyclotron-sector native domain checks and all field-selection/pusher logic are unchanged.

34.8.3 Eigenvalues, stability and fractional modes

LinearMapEigenAnalysis::analyze takes the final 4-by-4 matrix; the caller must first require a converged closed orbit. writeReport emits the eigenvalues, moduli, stability classification, residual and eigenbasis diagnostics, and mode phases. The COF input command calls these interfaces after verified convergence.

The full real nonsymmetric matrix, including transverse coupling, is analyzed with reference LAPACK’s dgeev through LAPACKE. No plane blocks are discarded. The routine works on a copy of \(A=S^{-1}MS\), retaining LAPACK’s error status. For normalized complex right eigenvectors, it verifies

\[ \frac{\|Av-\lambda v\|_2}{\max(\|A\|_F,|\lambda|)}\leq 10^{-10}. \]

For the zero matrix the unscaled residual is used. Nonfinite data, a failed LAPACK call or excessive residual cause an exception. Singular values from dgesvd provide the reciprocal condition of the real-packed eigenvector basis; this flags nearly dependent eigenvectors, including defective complex pairs. It is a scaled-basis diagnostic, not an eigenvalue error bound.

Classification uses a configurable unit-modulus band (default \(10^{-8}\)):

Unstable: at least one eigenvalue has modulus above \(1+10^{-8}\).
Non-unit-circle: some eigenvalue lies outside the band, but none grows. This is not certified as a stable conservative map; purely damped spectra are not incorrectly labelled growing.
Marginal: the spectrum is in the band but has real unit eigenvalues, a near-integer phase, or eigenbasis reciprocal condition below \(10^{-10}\).
Stable: the remaining well-conditioned, unit-circle complex-pair cases. This is a linear spectral classification, not a nonlinear stability proof.

Each positive-imaginary eigenvalue gives a representative phase \(q=\arg(\lambda)/(2\pi)\) and its conjugate branch \(1-q\). The report explicitly leaves the integer tune and oriented branch undetermined. For example, the same pair can represent fractional tune 0.11 or 0.89. Modes are sorted by this representative phase, not labelled radial/vertical or horizontal/vertical; coupled modes can exchange order in a parameter scan. No mode-continuation algorithm is implied. Off-unit-circle pairs do not receive stable tune values.

A unit-circle phase within \(10^{-6}\) turns of an integer is flagged explicitly. Exact real \(+1\) or \(-1\) eigenvalues are marginal and do not create oscillatory mode entries. Numeric phases reported for other marginal cases must not be interpreted as validated stable tunes. All diagnostic tolerances are configurable and require separate timestep/finite-difference studies for a real machine.

Tests cover stable rotations, coupled/scaled similarities, real reciprocal and complex unstable pairs, near-integer modes, repeated/defective spectra, and an analytic closed-orbit-solver-to-eigenanalysis workflow.

34.8.4 Closed-orbit launch handover

A named COF retains the full laboratory launch \((\mathbf r_0,\mathbf p_0,t_0)\), where positions are in metres and \(\mathbf p_0\) is mechanical momentum divided by \(mc\). For the solved section coordinates \((x,p_x,y,p_y)\) and fixed magnitude \(P_0=|\mathbf p_0|\), \(p_z=\sqrt{P_0^2-p_x^2-p_y^2}>0\); the section transform maps this state to the lab.

For TRACK initialisation, let \(Q\) map the orbit-local frame to the laboratory. Its +z axis is parallel to \(\mathbf p_0\), with transverse roll fixed by the shortest rotation from the COF section +z axis. Generated particles are placed as

\[ \mathbf r_i=\mathbf r_0+Q\boldsymbol\xi_i,\qquad \mathbf p_i=Q\boldsymbol\pi_i. \]

Here \(\boldsymbol\pi_i\) is the full local normalized momentum, not a deviation. This rigid change preserves distances, momentum norms and the distribution’s local correlations. It does not modify the integration algorithm or establish space-charge equilibrium. The reference remains the solved state even when the finite sample has a nonzero centroid. Scalar frame construction occurs on the host; particle kernels continue using the existing local frame and device views. All MPI ranks receive the same launch state; no new particle reduction is needed. Energy/species compatibility permits only floating-point roundoff (64 machine epsilons relative), not a physics matching tolerance.

The INITIALORBIT interface transfers this state without transferring the COF integrator or timestep. COF’s fixed-section ray return and TRACK’s Boris trajectory are different numerical calculations. TRACK also defines its return normal by \(\mathbf p_0\), whereas COF retains the selected section normal; these need not coincide for nonzero section-frame transverse launch momenta. Terminal-step localization controls distance from TRACK’s return plane, not transverse closure or agreement with COF. Establish TRACK timestep convergence independently of Newton residual, finite-difference convergence, and eigenvalue stability.

34.8.5 First native cyclotron closed-orbit benchmark

The opt-in OrbitThreaderTest.Cyclotron72ClosedOrbitBenchmark composes eight native CyclotronSector elements, the PSI magnetic map and mirrored trim coil, the shared RK4 ray tracker, fixed-section return map, Newton solver and LAPACK eigenanalysis. This C++ benchmark bypasses input parsing; separate command-level regressions exercise the COF input interface, native proton/electron SBEND rings and MPI status propagation. Supply the local map explicitly and run:

OPALX_COF_CYCLOTRON_MAP=/absolute/path/to/bfield.dat \
OMP_NUM_THREADS=1 OMP_PROC_BIND=false \
./omp-build/unit_tests/Algorithms/TestOrbitThreader \
  --gtest_filter=OrbitThreaderTest.Cyclotron72ClosedOrbitBenchmark

Without the map environment variable this optional machine benchmark is skipped. The static 72 MeV proton case starts from the first legacy ic.dat row, \((r,p_r,y,p_y)=(2.1314,-0.00024,0,0)\) in metres and mechanical \(p/(mc)\). The section is global \(Z=0\), with positive-\(Z\) return. The trim profile uses \(R_{min}=4.35\) m, \(R_{max}=4.47\) m, \(B_{max}=0.0014\) T and slope \(600\) m\(^{-1}\); its radial tails are retained. No RF or collective field is present.

Seven timestep levels, \(\Delta t=6/(50.65\times10^6\,720\,2^k)\) seconds, \(k=0,\ldots,6\), converge to numerical closed orbits with the unchanged \(10^{-10}\) componentwise closure tolerances. The finest radius is approximately 2.131503545 m and \(p_r/(mc)=-0.0002391852\). The final timestep halving changes the radius by about 18 nm; this discretization difference is larger than the nonlinear closure residual and must not be confused with it.

The radial Jacobian is sensitive to both timestep and central-difference step. At the finest timestep its positive phase is approximately 0.1574, but the pair modulus still differs from unity by order \(10^{-4}\)–\(10^{-3}\) for small derivative steps. The vertical complementary phase is approximately 0.89338. These are diagnostic phases, not a certified stable spectrum: the default \(10^{-8}\) unit-circle tolerance is unchanged. The coarse-to-fine radial modulus alternates between artificial growth and damping. Resolving the differentiated tracking error is required before interpreting a stability label physically; the field map’s piecewise interpolation is a candidate to investigate, not a confirmed explanation. This benchmark does not validate the separate segmented transfer-map seam treatment or an ISIS/DBA-ring closed orbit.

A displaced launch (1 mm radial, 0.5 mm vertical, with momentum offsets) recovers the same orbit within \(10^{-9}\) in each coordinate’s units. An independent 0.1 mm amplitude, 100-nominal-turn spectral run gives 1.155 and 0.890, also the previous old-OPAL spectral values. Do not directly subtract these from physical per-return map phases: the legacy estimator assigns its entire sampled span exactly 100 turns. If \(t_s\) is the last-minus-first sample time and \(T_{rev}\) is the measured COF period, compare the legacy spectral value against

\[ \nu_{legacy,expected}=\nu_{return}\frac{t_s}{100T_{rev}}. \]

Here \(T_{rev}\simeq118.6675\) ns and the factor is approximately 0.9975584, including the omitted final sampling interval. The predicted legacy-axis phases are approximately 1.15455 and 0.89120, within one frequency bin (0.0025) of the measured peaks. This independent normalization explains the comparison offset without fitting the tunes, changing the legacy estimator or relaxing the eigenvalue stability threshold. No new old-OPAL run is implied.

34.8.6 Higher-energy checks with trim coils

Set OPALX_COF_CASE=400 or OPALX_COF_CASE=590 in the same benchmark invocation (the default remains 72). The 400 case uses the actual 399.995 MeV ic.dat row. The table ends at 529.993 MeV, so the 590 case first performs fixed-energy Newton continuation from that row in increments no greater than 5 MeV. This is a 590 MeV coasting orbit, not a closed orbit during acceleration.

Both cases use the same seven timestep levels, derivative-step sweep, trim coils, displaced-launch check and 0.1 mm spectral excitation described above. The finest solutions are:

Kinetic energy [MeV]	Radius [m]	\(p_r/(mc)\)	Radius change on final timestep halving
399.995	4.0203122462	0.1966177730	2.0 nm
590	4.4014800681	0.2101917477	16.5 nm

Both displaced starts recover the same orbit within \(10^{-9}\) in coordinate units. Fresh closure residuals remain below \(10^{-10}\) and finest-step relative energy drift is below \(3\times10^{-14}\). These are numerical checks, not bounds on magnetic-map or physical-model errors.

New old-OPAL RK4 runs use the OPALX closed-orbit coordinates, the same map, explicitly full-azimuth trim coil, 0.1 mm excitation, 2880 steps per nominal turn and 100 nominal turns. The legacy Lomb estimator applied to the old-OPAL signals gives the following comparison on the legacy frequency axis, with spacing 0.0025:

Energy [MeV]	Mode	Map phase, converted	OPALX spectral	Old OPAL spectral
399.995	Radial	1.488155	1.4875	1.4875
399.995	Vertical	0.872986	0.8725	0.8725
590	Radial	1.326722	1.3275	1.3275
590	Vertical	1.173548	1.1725	1.1725

All four map-to-spectral differences are smaller than half a frequency bin. Eigenvector weights identify the planes in this uncoupled median-plane case; the spectral result supplies the integer/conjugate branch. Sorted eigenvalue order alone is not a plane label. The corresponding selected physical per-return branches are about 1.489114/0.873549 and 1.328672/1.175272, respectively; eigenanalysis alone still cannot select these branches.

At derivative steps \((10^{-5}\,\mathrm{m},10^{-6},10^{-5}\,\mathrm{m},10^{-6})\), the finest radial moduli are approximately 0.999660 and 0.998556. They remain sensitive to differentiation and timestep: the tune agreement does not certify conservative linear stability. No unit-circle tolerance has been relaxed. sandbox/cyclotron/compare_cof_energies.py reproduces the matched old-OPAL run and comparison from a successful benchmark log; inputs are stored under sandbox/cyclotron/opal/cof-tunes/ and results under sandbox/cyclotron/cof-results/.