How Fast Can High-Fidelity Gear Rattle Be Simulated, and Where Does Trajectory Simulation Give Way to Structure?

1. Introduction

Gear rattle is the audible signature of an unloaded or lightly loaded gear pair whose teeth repeatedly cross the backlash gap and impact on drive and coast flanks. It is a canonical piecewise-smooth mechanical system: between impacts the mesh contributes no restoring force, while inside contact the teeth engage through a time-varying, load-dependent mesh stiffness. The resulting equations of motion switch between distinct smooth vector fields at state-dependent events, and the closed-loop response is strongly nonlinear, exhibiting multi-flank impacts, subharmonics, and chaos. This behaviour has been documented and modelled for more than three decades: the periodically forced piecewise-linear oscillator provides the reduced archetype [1]; impact-map and gearbox studies traced the origin of chaotic rattle [2]; the gear-pair model with periodic stiffness and backlash established the reference nonlinear dynamics [3]; motor-driven configurations added the actuation loop [4]; and both refined models [5] and combined experimental–theoretical work [6] have continued to sharpen the picture of backlash-driven vibro-impact. The physics is, in this sense, well understood. What is not settled is a pair of engineering questions that the title of this paper poses directly.

How fast can high-fidelity gear rattle be simulated? By high-fidelity we mean a model that retains the features practitioners actually care about — the three-state backlash (drive contact, free flight, coast contact), a time-varying mesh stiffness excitation, tooth friction, and a flexible reduction of the gear bodies rather than a lumped two-mass caricature. Simulating such a model faster than real time is a prerequisite for hardware-in-the-loop powertrain testing [7] and for the real-time multibody workflows now being catalogued systematically [8]. The naive route — a stiff adaptive integrator chattering across every flank crossing — is far too slow, and it is also needlessly inaccurate, because the exact solution of each smooth segment is available in closed form.

Where does trajectory simulation give way to structure? Speed is not without consequence. Rattle is chaotic, so a computed trajectory is trustworthy only for a finite predictability horizon set by the largest Lyapunov exponent [9,10]. Past that horizon a faster point-trajectory buys essentially nothing: the meaningful object of computation ceases to be the trajectory and becomes the invariant structure — periodic-orbit skeletons, their stability, and the statistics of the attractor. One of our aims is to locate this crossover quantitatively and to show that beyond it real-time effort is better spent on structure than on integrating one more decorrelating path.

We are deliberately explicit about what is and is not new here. The enabling tools are all prior art. Segment-exact propagation by exponential integrators is a mature body of work [11,12], including breakpoint-stepping variants for switched systems [13]; the general theory and numerics of nonsmooth integration are likewise established [14,15]. Linearising the flow across an event uses the saltation matrix [15,16], and the performance ceiling of any such kernel is read off the Roofline model [17]. The predictability argument rests on finite-time Lyapunov analysis and its recent extensions [9,10]. We introduce none of these. Our contribution is instead to assemble them into a coherent, measured account of a single hard problem, and specifically:

(i)

high-fidelity super-real-time integration of a nonsmooth gear-rattle model, achieved by propagating each backlash state with a segment-exact exponential step and fusing many mesh cycles into a single spectral map — the constant-contact fusion advances one mesh cycle in a single matrix–vector product (∼61 ns) and jumps ${10}^6$ cycles in 634 ns, while the sparse-regime event-driven path reaches $16\times$ real time with a brute-force cross-check at $7\times {10}^{-15}$;

(ii)

the observation that single-trajectory depth and ensemble breadth hit the same wall — both saturate against memory bandwidth rather than arithmetic, so a 56-core rattle atlas returns $16.6\times$ speedup (about $30\%$ efficiency, not $56\times$), the same bandwidth ceiling that also caps the deepest single-path fusion;

(iii)

the identification of a Lyapunov wall — with a largest exponent of ∼0.55 per rattle cycle, exponential compute buys only a linear gain in predictable horizon, so beyond the wall real time must shift from trajectories to structure (Floquet skeletons, basin and attractor statistics);

(iv)

honest negative results — where fusion does not pay, where added modal fidelity is not worth its cost, and where the flexible model admits no stable period-1 orbit at all ($0/80$ operating points), so the accessible range is entirely chaotic.

These findings are summarised by the flagship layered-ceiling map of Figure 1, which collects the achievable real-time factor as a function of model fidelity and integration strategy and marks the bandwidth and Lyapunov walls that bound the whole design space. As a headline calibration of the two questions above: at fixed accuracy the exact-fused kernel reaches an error of $1.2\times {10}^{-14}$ at $5\times {10}^{-4}$ ms/solve, against a high-order adaptive reference (Vern7 at tolerance ${10}^{-10}$) that returns error $1.5\times {10}^{-9}$ at $68.9$ ms — roughly ${10}^5$ times faster and${10}^5$ times more accurate, because exactness on each smooth segment removes both the cost and the local error of adaptive stepping.

The remainder of the paper is organised as follows. Section 2 states the nonsmooth gear-rattle model and its flexible reduction. Section 3 develops the segment-exact and spectral-fusion kernels and the event-driven path (Algorithms A1–A3), reports the tier ladder (Table 1) and the event-path progression (Table 2), and closes with the Roofline analysis of the shared depth/breadth memory wall (Figure 4). Section 4 turns to structure past the Lyapunov wall — the predictability crossover (Figure 7, Figure 8), the chaos atlas and bifurcation sweep (Figure 9, Figure 10), the Floquet skeleton (Figure 11), the NVH and phase-space signatures (Figure 12, Figure 13), and the degree-of-freedom convergence (Figure 14). Section 5 collects the negative results (Table 5), Section 6 the related work, and Section 7 concludes. Derivations, the optimisation history, validation, and the computing environment are deferred to the appendices.

2. A High-Fidelity Nonsmooth Gear Model

We study a spur gear pair driven near a mesh-order resonance, modelled at a fidelity that resolves the deformation-compatible tooth compliance, the flexibility of the connecting shafts, tooth friction, and the backlash-induced loss and re-establishment of contact. Figure 2 lays out the five concrete test cases used throughout the paper side by side — from a 2-DOF torsional pair to the flexible-shaft, three-state-backlash configuration — together with the shared mesh-stiffness engine that feeds all of them; the tiers differ only in which degrees of freedom are released. The full derivation, parameter set, and the exact segment propagator are deferred to A; here we state only the structure that the rest of the paper exploits, namely that the model is piecewise-linear time-invariant between well-defined switching events.

Deformation-compatible time-varying mesh stiffness.

The mesh compliance is obtained by an improved potential-energy method [18,19], in which each engaged tooth pair contributes bending, shear, axial and Hertzian-contact energy, summed over the instantaneous contact lines to give the time-varying mesh stiffness (TVMS) $k_m\left(\varphi \right)$ as a function of mesh phase $\varphi =\Omega t$. Profile modification and misalignment are admitted through the deformation-compatibility condition: under a misalignment measured by a single effective angular deviation ${\epsilon }_{\mathrm{eff}}$, part of the flank may unload, so that the tangent mesh stiffness $k_m(\varphi ;{\epsilon }_{\mathrm{eff}},T)$ and the load centroid $z_c(\varphi ;{\epsilon }_{\mathrm{eff}},T)$ depend jointly on ${\epsilon }_{\mathrm{eff}}$ and the transmitted torque T. Both the loaded static transmission error and this coupling between backlash, transmission error and TVMS follow the pathway documented in [20]. In the dynamic model $k_m$, $z_c$ and the static transmission error $e\left(\varphi \right)$ enter as prescribed, tabulated functions of mesh phase; within a mesh sub-interval they are treated as constant, which is the first source of piecewise structure. Two verified properties of this excitation pathway matter later. First, the misalignment loop is self-consistent: closing the slow feedback (dynamic response $\to {\epsilon }_{\mathrm{eff}}\to$ stiffness table) at ${\epsilon }_0=1$ mrad and 50 N m converges to a mean mesh stiffness $41\%$ below nominal, matching the $42\%$ partial-contact drop that the static table predicts independently. Second, parametric stiffness variation alone cannot produce rattle: under pure $k_m\left(\varphi \right)$ excitation the separation margin scales with the response itself and the teeth never leave contact — it is the displacement excitation $e\left(\varphi \right)$ near a mesh resonance that opens the gap, so the rattle studied here is STE-driven, not parametrically driven.

Flexible-shaft modal degrees of freedom.

The pinion and gear are carried on flexible shafts, reduced by modal truncation to N retained shaft/disc modes per member. The generalized coordinate is $q={[{\eta }^{\top },\delta ]}^{\top }$, collecting the modal amplitudes $\eta$ together with the mesh coordinate; the dynamic-transmission-error contribution of every mode enters the mesh only through a rank-one coupling vector $w_c={\tilde{w}}_0+z_c{\tilde{w}}_1$, so that misalignment (through $z_c$) rotates a single coupling direction rather than restructuring the operator. Truncated modes are handled by residual (mode-acceleration / attachment) augmentation to restore the quasi-static flexibility they carry [21,22], consistent with finite-element geared-rotor practice [23] and component-mode reduction [24]. The retained system has dimension $dim=2N+1$.

Tooth friction.

Sliding friction on the flanks is included as a Coulomb law $s=\mu F\Lambda \left(\varphi \right)$, where $F=k_m\delta +c_m\dot{\delta }$ is the mesh load, $\mu$ the friction coefficient and $\Lambda \left(\varphi \right)$ the sign/lever function that reverses at the pitch point. Because $\Lambda$ is fixed on each sub-interval, the friction generalized force is affine in the state there; it renders the segment operator non-symmetric (non-conservative) but does not break its linearity.

Three contact states and the backlash half-gap.

Backlash of half-gap b makes the mesh restoring law a three-state piecewise-linear map [1,3,4]: a free state ($\left|\delta \right|<b$, no tooth load), a drive-side contact state ${\mathrm{contact}}^{+}$ ($\delta \ge b$) and a back-side contact state ${\mathrm{contact}}^{-}$ ($\delta \le -b$). The switching surfaces $\delta =\pm b$ separate three affine vector fields, giving a nonsmooth (Filippov-type) mechanical system [14,15]. Only in the contact states do the mesh stiffness and damping act; in the free state the mesh transmits no load and the two members interact only through the shafts. Crossing a switching surface is an impact-like event that must be located, not stepped over.

Augmented piecewise-linear state.

Collecting the modal and mesh coordinates with an appended constant “carry” coordinate that transports the external torque and the static-transmission-error excitation, we obtain an augmented state $x\in {\mathbb{R}}^{dim+1}$ whose evolution on each segment $\sigma$ (a fixed contact state and a fixed mesh sub-interval) is linear time-invariant,

$$\dot{x}=A_{\sigma }x,A_{\sigma }=A_{\sigma }\left(\mathrm{state},k_m,c_m,z_c,\mu ,e\right),$$

(1)

with the affine excitation absorbed into the carry coordinate so that $A_{\sigma }$ is a genuine matrix and $x\left(t\right)=exp\left(A_{\sigma }(t-t_0)\right)x\left(t_0\right)$ within the segment. A trajectory is thus a concatenation of exact LTI flights joined at two kinds of breakpoints: mesh sub-interval boundaries, where $A_{\sigma }$ changes because $(k_m,z_c,e)$ change; and contact-state switches at $\delta =\pm b$, where both the active vector field and (through the damping discontinuity) the transition map change. The exact per-segment propagator $exp\left(A_{\sigma }\Delta t\right)$, the switching-surface event localization, and the state-transition (saltation) correction at $\delta =\pm b$ are given in Appendix A and summarized as Algorithm 1; the whole of the present paper is an investigation of how far this exact piecewise-linear structure can be pushed before physics, rather than arithmetic, sets the limit.

3. Exact Propagation and the Real-Time Frontier

The model of Section 2 is piecewise-affine: within each contact state the dynamics are linear, and the only nonlinearity is the switching between states. This structure is not new, and neither is the observation that a linear segment can be advanced by its matrix exponential rather than by a general-purpose stepper; exponential integrators and their piecewise-linear specialisations are mature technology [11,12,13], and closed-form treatments of clearance and piecewise-linear oscillators have been given semi-analytically [25,26]. Our contribution here is not a new integrator but a measurement: applied to a high-fidelity nonsmooth gear-rattle model and pushed to its limits, exact propagation runs faster than real time, and we can say precisely how fast, where the speed comes from, and where it stops. Throughout, “×” denotes the wall-clock speed-up over real time, and all timings are single-core unless stated otherwise.

3.1. Constant-contact Regime: a Mesh Cycle in One Matvec

When the tooth pair stays in contact, the mesh cycle is a single affine segment and its time-T flow is one constant propagator $\Phi =e^{AT}$ acting on the state, precomputed once and reused for every cycle. Advancing a whole mesh period then costs a single matrix–vector product, $\sim 61$ ns for the low-dimensional lumped model. Because contact is smooth over many cycles, the composition ${\Phi }^n$ can be assembled by repeated squaring: dyadic power-jumping advances ${10}^6$ mesh cycles in 634 ns, a $\sim 38000\times$ speed-up over cycle-by-cycle marching, at the cost of storing ${log}_{2}n$ dyadic factors. The Tier ladder of Table 1 reports the small-step (non-fused) real-time factor as the lumped dimension grows: $3.3\times$ at 2-DOF, $2.6\times$ at 6-DOF, $2.0\times$ at 12-DOF, $1.3\times$ at 20-DOF, and $0.53\times$ at the 40-mode flexible case (state dimension 81), the first point that falls below real time. The fused propagator recovers all of this and more: on the Tier-5 flexible configuration with Coulomb friction ($\mu =0.1$), constant-contact simulation moves from $0.4\times$ small-step to $8.0\times$ with $\Phi$-segment propagation and to $487.7\times$ with $\Phi$ fusion, while segment- and fused-propagation agree to $5.7\times {10}^{-12}$ and the transmission-error RMS differs by only $5.8\%$ ($5.89$ vs. $6.25\mu$m) against the reference model. Constant contact is therefore essentially free; the cost of the problem lives entirely in the switching.

3.2. Nonsmooth Regime: Closed-Form Event Location and Segment Fusion

Under rattle the trajectory crosses backlash boundaries, and the wall-clock cost is set by how expensively each crossing is located. Event-driven integration of nonsmooth systems is standard [14,27], as is the use of switch detection to place discontinuities exactly [28], and the impact-map view of gearbox chaos dates to Hongler and Streit [2]. What matters for speed is that inside an affine segment the guard function is a known combination of exponentials, so the next crossing is the smallest root of a scalar transcendental equation rather than the output of a bisection over full state steps. Table 2 traces the resulting optimisation path for strong rattle at $\Omega =1200$ (tooth-passage time $t_{\mathrm{tooth}}=262\mu$s). A naive segmented propagator (V1) needs 371 matvec per cycle ($5\times$); reusing propagators (V2) drops this to 193 ($8\times$); and closed-form spectral event location (V3) reaches 74 matvec plus 124 cheap scalar guard evaluations per cycle, still $8\times$ but now dominated by the physics of the impacts rather than by root-finding overhead. In the sparse-rattle regime ($\Omega =1400$) the same V3 path takes $29.1$ matvec and $14.6\mu$s per cycle for a $16\times$ real-time factor, and cross-checks against a brute-force reference to $7\times {10}^{-15}$. The segment-fusion procedure that composes consecutive same-state segments is given as Algorithm 2 in Appendix C, with the spectral event solver in Algorithm 1. Exact event location is thus what keeps the nonsmooth regime above real time, but its per-cycle cost is bounded below by the event rate, which is a property of the trajectory, not of the algorithm.

3.3. High-Dimensional Flexible Configuration and the Memory-Bandwidth Wall

The most demanding configuration couples flexible modal DOFs, a three-state backlash, and Coulomb friction, with state dimension $2N+1$ for N retained modes. Here two facts dominate. First, the constant-contact propagator is far better conditioned than the free-flight one—$\mathrm{cond}\left(V_1\right)\approx {10}^3$ in contact versus $8.1\times {10}^5$ free—and the spectral eigenbasis reconstructs the state to $1.5\times {10}^{-11}$, so exactness is not the obstacle. The spectral event table costs 18 MB, against 190 MB if the same coverage were stored as dyadic power factors, and the spectral solution tracks brute force at the expected first-order rate ($\sim {10}^{-3}$ at the working step) with an event rate of $2.7$–$2.9$ per cycle. Second, and decisively, the propagator matvec is memory-bound. Measured matvec cost grows as $3.16\mu$s at dimension 31, $3.7\mu$s at 41, and $10.7\mu$s at 81: the jump is the L1 crossing, since an $81\times 81$ double-precision matrix is 52 KB and no longer fits the 32 KB L1 cache, placing the knee near 22 DOF. Model-order reduction is therefore not accuracy hygiene but a throughput device—keeping the working matrix resident in L1—and the relevant reduction machinery (mode acceleration, modal truncation augmentation, residual-attachment and dual Craig–Bampton synthesis, and data-driven reduction for piecewise-linear structures) is established prior art [21,22,24,29]. The optimisation chain of Figure 3 and Table 3 shows the payoff near resonance: the raw $N=40$ model (dimension 81) runs at $0.20\times$ ($1282\mu$s, 120 matvec per solve), below real time; truncating to $N=20$ (dimension 41) with a simple projection reaches $0.57\times$; adding modal truncation augmentation (Algorithm 3) holds accuracy at $0.54\times$; and combining MTA with segment fusion reaches $1.53\times$ ($171\mu$s, 52 matvec)—the same physics now above real time. The MTA residual flexibility amounts to a $0.19\%$ static softening ($c_{\mathrm{res}}k=1.9\times {10}^{-3}$), and the displacement-RMS is preserved to better than $2\%$ across the chain ($0.0167$ mm full-$N40$, $0.0164$ simple-$N20$, $0.0165$ MTA-$N20$). Internal consistency is exact—fusion versus spectral agree to $4.0\times {10}^{-15}$ and spectral versus brute to $5.06\times {10}^{-6}$. Table 4 collects the full-configuration figures. The fidelity–cost trade is lopsided: fidelity is cheap until the working set leaves cache (Figure 5a), after which every added mode is paid for in bandwidth.

3.4. The Unified Wall

The L1 crossing of Section 3.3 is one face of a single limit. Making one trajectory deeper (more modes, longer horizons) and running more trajectories at once (an ensemble sweep) are throttled by the same resource, and the Roofline model [17] together with memory-bound performance prediction for stencil-type kernels [30] names it: at the propagator’s low arithmetic intensity the kernel sits under the bandwidth roof, not the compute roof. Figure 4 shows both faces meeting the same ceiling. On the ensemble side we compute a rattle atlas of 80 operating points (8 loads × 10 speeds) with 8 contact-state tables preloaded: the 56-core parallel run takes $0.66$ s against $11.0$ s serial, a $16.6\times$ speed-up—not the $56\times$ that core count would suggest, but roughly $30\%$ parallel efficiency, because the shared propagator and table traffic saturate memory bandwidth long before the cores saturate. Depth and breadth are throttled by one wall, and neither exponential-in-cost strategy escapes it.

Because hardware performance counters are unavailable under WSL2 (the Microsoft kernel does not virtualize the PMU), we characterize the memory system directly, by microbenchmark, using the production matvec kernel itself (Figure 5). Three measurements pin the wall down. First, a memory mountain: applying the solver’s own column-major saxpy kernel to a rotating set of operators — emulating how a mesh cycle walks through its spectral table — shows a cost that is flat while the rotating working set stays inside the 35 MB socket L3, then jumps $2.7\times$ (dim 81: $1.8\to 4.8\mu$s per matvec) once it spills to DRAM. At matched residency the per-operator cost scales essentially arithmetically with ${\dim }^2$; the hard cliff belongs to the table, not to the single operator. Second, aggregate streaming bandwidth saturates at ≈32 GB/s with only 14 of the 56 hardware threads — barely $3\times$ one core’s $10.2$ GB/s — and then decreases to 22 GB/s at 56 threads under NUMA and oversubscription pressure. Third, thread scaling of the full spectral solver on a shared table saturates at 25–$28\times$ on 56 threads, and the harder 8-table atlas workload at the $16.6\times$ already quoted. The traffic arithmetic closes the argument: a raw-N40 mesh cycle performs 120 matvecs against 52 KB operators, about $6.3$ MB of table traffic per cycle, so the cycle is a bandwidth workload rather than a flop workload. Appendix F maps every working set of the solver onto the machine’s cache hierarchy explicitly (Figure A2); read this way, MTA reduction, segment fusion, and the spectral-table layout are all the same move — relocating this traffic upward in the hierarchy.

3.5. Accuracy with Speed

None of the above trades accuracy for speed; exact propagation delivers both at once. Figure 6 places the fused exact method on the work–precision plane against a high-order adaptive reference. The exact-fused solution reaches error $1.2\times {10}^{-14}$ at $5\times {10}^{-4}$ ms per solve, while Vern7 at tolerance ${10}^{-10}$ reaches error $1.5\times {10}^{-9}$ at $68.9$ ms per solve—about ${10}^5\times$ faster and about ${10}^5\times$ more accurate simultaneously, because exactness removes local truncation error rather than balancing it against step count. This is the strongest statement the method supports on a single trajectory. It is also the point at which the question stops being “how fast” and becomes “how far,” which Section 4 takes up: when accuracy is essentially machine-exact and speed is bandwidth-bounded, the remaining limit is set not by the solver but by the physics of the trajectory itself.

4. Where Trajectories Give Way to Structure

The preceding sections answered “how fast”: exact piecewise propagation and its fused variants place a high-fidelity nonsmooth gear model deep into super-real-time territory. This section answers “where” – where the value of ever-faster trajectory computation runs out, and what real-time compute should produce instead. The boundary is not set by our solver, our hardware, or our memory bandwidth; it is set by the physics of the system itself.

4.1. The Lyapunov Wall

Rattle in the backlash regime is chaotic, and chaos imposes a finite horizon of predictability that no amount of computation can push back. For the strong-rattle 6-DOF configuration we estimate a largest Lyapunov exponent of $\lambda \approx 0.55$ per rattle cycle: an initial state separation of $4.57\times {10}^{-4}$ at three cycles grows to $2.09\times {10}^{-2}$ by ten cycles, consistent with exponential divergence at that rate. The consequence is the central asymmetry of chaotic simulation. Because trajectory error grows as $e^{\lambda t}$, the horizon $T_{\mathrm{pred}}$ over which a computed trajectory tracks the true one to a fixed tolerance scales only as $T_{\mathrm{pred}}\sim {\lambda }^{-1}ln({\epsilon }_{\mathrm{tol}}/{\epsilon }_0)$. Reducing the numerical error floor ${\epsilon }_0$ – or, equivalently, spending exponentially more compute to integrate more accurately or further – buys only an additive, logarithmic extension of the predictable horizon (Figure 7). This is a hard physical wall: exponential computational effort yields at best linear gain in trustworthy horizon, and beyond a few tens of cycles the pointwise trajectory is not merely inaccurate but, for the underlying continuous system, unknowable in any practical sense. The same asymmetry motivates finite-time Lyapunov diagnostics of predictability [9], reliability studies of computed chaotic solutions [31], and recent work on how far the horizon of a chaotic simulation can legitimately be trusted [10]. Our super-real-time integrator does not evade this wall – it reaches it faster.

The methodological conclusion follows directly. Once real-time computation has outrun the Lyapunov horizon, additional throughput spent on a single long trajectory is wasted, because that trajectory has already lost its individual meaning. The productive use of the surplus is to compute invariant structure – the statistical and geometric skeleton that is stable under the very sensitivity that destroys individual orbits. Past the wall, real time should stop chasing trajectories and start mapping structure. The remainder of this section demonstrates three such structural products, all computed within a real-time budget on 56 CPU cores.

This crossover can be made quantitative on the fully resolved flexible model (Figure 8, $N=40$). A perturbed twin’s normalised state separation $\delta \left(n\right)/R$ grows to the attractor scale by a predictability horizon $n_h\approx 17$ mesh cycles (Benettin exponent $\lambda \approx 0.9$cyc${}^{-1}$, order-one and consistent with the 6-DOF skeleton). Over the same trajectory lengths, the relative standard error of a time-averaged structural observable (the mean contact force) falls as $n^{-1/2}$, from order unity toward $0.087$ by 400 cycles. The two curves cross near $n_h$: to the left the pointwise orbit is the better-determined object and integration is worthwhile; to the right the orbit has decorrelated while the statistical structure has sharpened, so the real-time budget is better spent on structure. The horizon itself is a property of the resolved model: coarser modal truncations over-estimate the exponent (Section 4.5), so $N=40$ gives the honest, and longest, horizon.

4.2. The Chaos Atlas

The first structural product is a global map of dynamical regime over the operating envelope. We sweep an $8\times 10$ grid of 80 operating points (eight load levels × ten speeds) and, at each point, estimate the largest Lyapunov exponent from the exact-propagation trajectory to classify the long-term regime (Figure 9). A structural property of the spectral tables makes the speed axis free: the block exponential is evaluable at any elapsed time, so one table per load serves every nominal speed — the $\Omega$ sweep costs no re-tabulation at all, whereas a dyadic power table is welded to the single time step it was squared for. The result is unambiguous: zero of the eighty points settle onto a stable period-1 attractor ($0/80$ stable). Within the resolved envelope the flexible-mesh rattle model is all-chaotic; there is no quiet corner where a single trajectory would suffice. This is itself an argument for the structural view – the atlas is the deliverable, and it is meaningful precisely because it is an ensemble statistic rather than a set of individual orbits, each of which is unpredictable past the wall of §Section 4.1. The all-chaotic finding is consistent with the classical picture of backlash-driven gear dynamics as generically chaotic in the impacting regime [4,32].

Computing the atlas exercises ensemble breadth rather than trajectory depth, and it is instructive that both hit the same ceiling. Running the 80 points across 56 cores with the eight per-load spectral tables preloaded takes $0.66$ s, against $11.0$ s serial: a speedup of $16.6\times$, not the nominal $56\times$. The parallel efficiency is roughly $30\%$, and the shortfall is not scheduling overhead but memory bandwidth – the independent integrations contend for the same shared path to the preloaded tables, exactly the bottleneck that governs single-trajectory depth. The ensemble is embarrassingly parallel in principle, and the many-independent-ODE literature on GPUs and CPUs [33,34,35] treats such sweeps as the canonical throughput case; but our workload is table-bound rather than arithmetic-bound, and its behaviour is set by the memory-bound Roofline regime [36] rather than by peak flop counts. Depth and breadth are two faces of one wall.

A one-parameter cut through the atlas sharpens the picture. Sweeping the nominal speed ratio $\Omega /{\Omega }_0$ from $0.55$ to $1.70$ and stroboscopically sampling the mesh deflection once per cycle yields the Poincaré bifurcation diagram of Figure 10(a); the largest Lyapunov exponent over the same sweep is in Figure 10(b). The response is chaotic across almost the whole range (${\lambda }_{max}\gtrsim 2$), interrupted only by a narrow near-periodic window at $\Omega /{\Omega }_0\approx 1.23$ where ${\lambda }_{max}$ dips to $\approx 0.1$ and the Poincaré set collapses to a thin band. We probed this window explicitly for multistability: 32 stratified initial conditions (perturbation amplitudes spanning four decades) at the window centre and at its two shoulders, plus a control point at $\Omega /{\Omega }_0=1$ — 128 trajectories in all — every one of which converged to a chaotic attractor ($\lambda \ge 0.11$ even at the window centre). The dip is weak chaos, not a periodic island: no coexisting periodic attractor, and hence no resolvable basin structure, exists even in the most quiescent corner of the sweep. This closes the multistability question for this configuration and reinforces the verdict of the two-dimensional atlas: the structural, rather than the trajectory, view is the correct one everywhere in the resolved envelope.

4.3. Unstable Periodic Orbits and Their Bifurcation Type

The second structural product is the organizing skeleton of the chaos itself: the unstable periodic orbits (UPOs) embedded in the attractor and their local stability type. Where individual trajectories are unknowable, the UPOs and their invariant manifolds are fixed, computable objects that shape the statistics of the flow. We locate them by Newton shooting and continue them in operating parameter, following standard multiple-shooting / Newton–Krylov continuation of periodic orbits [37,38] adapted to the piecewise-smooth setting [39,40]. The essential ingredient for a nonsmooth model is correct linearization across the impact and backlash-crossing events: the monodromy matrix is assembled from the smooth-segment fundamental solutions bridged by saltation (jump) matrices at each event [15,16], and the Floquet multipliers are read off the resulting one-period map [41,42]. Full method detail — including the neutral-direction regularization that the gear pair’s rotational symmetry forces on the Newton solve — is given in Appendix E; here we report the structural verdict.

For the strong-rattle 6-DOF case, Newton shooting converges on a period-1 UPO – the skeleton about which the chaotic band is organized. Its monodromy spectrum carries a neutral pair at unity together with a dominant complex-conjugate multiplier $-0.373\pm 0.897\mathrm{i}$ of modulus $\left|\lambda \right|=0.972$. This pair sits just inside the unit circle, at a distance of $0.028$ from it (Figure 11): the period-1 skeleton is weakly stable and poised immediately below a Neimark–Sacker (torus) bifurcation, the codimension-one route by which a complex pair crosses the unit circle. Such near-torus and grazing-induced Neimark–Sacker structure is precisely what has been catalogued for impacting oscillators [42,43]. In the sparse-rattle regime the corresponding period-2 orbit is strongly unstable, with a leading multiplier ${\left|\lambda \right|}_{max}=1.92$ outside the unit circle. The two cases together sketch how the periodic skeleton loses stability across the envelope, and they do so with objects that are reproducible to numerical precision – unlike the trajectories they organize.

One measured caveat deserves emphasis, because it is easy to get wrong in any nonsmooth system. The saltation-corrected monodromy linearizes the flow along the fixed event sequence of the orbit; an independent ground-truth test — kick the converged orbit by ${10}^{-5}$ and iterate the true map — shows the strong-rattle perturbation exploding to $6.3\times {10}^{-2}$ within a single return and growing by $1.51\times$ per map thereafter, even though $\left|\lambda \right|=0.972<1$ promises weak stability. The discrepancy is not a saltation error: the kick changes which events occur, and this grazing-induced reshuffling of the event sequence [40,44,45] lies outside any linearization about the unperturbed sequence. Floquet multipliers in a nonsmooth system are therefore trustworthy only as far as the event sequence is robust — one more reason why, past the wall, the Lyapunov exponent rather than the multiplier is the honest measure of sensitivity (Appendix E quantifies both cases).

4.4. Engineering NVH Signatures

The third structural product is the one an NVH engineer actually consumes: statistical and spectral signatures of the vibration, which are stable ensemble properties even though the underlying time series is chaotic. From the dynamic transmission error (DTE) we obtain an RMS of $\approx 13\mu \mathrm{m}$ against a peak-to-peak of $\approx 60\mu \mathrm{m}$; the crest-like ratio $\mathrm{pp}/\mathrm{RMS}\approx 4.6$ quantifies the strongly non-Gaussian, impact-dominated character of rattle (Figure 12). Spectrally, the mesh-fundamental harmonic sits at $h_1=5\mu \mathrm{m}$, matching the imposed static-transmission-error excitation and confirming that the model reproduces the correct forcing line rather than a numerical artefact. The peak contact force scales from $6.0$ to $9.6\mathrm{kN}$ as load rises from $5\%$ to $30\%$, a monotone load signature usable for durability screening. These are the quantities that survive the Lyapunov wall: aggregated over the chaotic ensemble they are the correct real-time deliverable once pointwise prediction has ceased to be meaningful, and they connect the fast nonsmooth solver of the previous sections to the engineering questions – rattle severity, spectral content, load-dependent contact loading – that motivate the simulation in the first place.

The same motion admits the two classical phase-space and spectral diagnostics of nonlinear dynamics (Figure 13). The phase portrait $(\delta ,\dot{\delta })$ with its once-per-cycle Poincaré section (Figure 13a) shows an orbit that fills a bounded region and crosses the backlash band $\pm b$ repeatedly, with a Poincaré set that is a fractal-like scatter rather than a finite point cloud — the geometric fingerprint of the chaos. The power spectra of the DTE and contact force (Figure 13b) place the mesh fundamental and its harmonics on a broadband, subharmonic-rich floor, the spectral fingerprint of the same non-periodic response. Both are ensemble-stable even though the trajectory that generates them is not.

4.5. Degrees of Freedom: Response and Chaos Convergence

A recurring question for a reduced flexible model is how many modes are enough. The question has a subtle answer here, because two different response measures converge at different rates. Expanding the model from $N=15$ to 40 retained modes (state dimension $31\to 81$), the impact-severity metric ${\delta }_{\mathrm{RMS}}$ converges quickly ($14.8$, $14.5$, $13.7\mu$m; Figure 14a), so a coarse truncation already captures the amplitude. The largest Lyapunov exponent does not: coarse truncations over-estimate the sensitivity, returning ${\lambda }_{max}\approx 2.3$ per mesh cycle at $N=15$ and 20, and only the resolved $N=40$ model settles to ${\lambda }_{max}\approx 0.7$ (Figure 14b), consistent with the 6-DOF skeleton’s $\approx 0.55$ per rattle cycle. Enough degrees of freedom are therefore needed not for amplitude fidelity, which is cheap, but to avoid a spurious over-prediction of chaotic sensitivity — and it is the resolved exponent that sets the honest predictability horizon of Figure 8. This is the “few-to-many degrees of freedom” axis made quantitative: adding modes buys not more amplitude but a truer measure of unpredictability.

5. Honest Limits and Negative Results

The preceding sections reported where the method wins. This section reports, with equal care, five places where a back-of-the-envelope estimate was wrong by a factor we could not have derived on paper. Table 5 collects them. In every case the wall was measured, not predicted, and in every case it is instructive: the gap between the naive estimate and the silicon is itself the physics — or the hardware — telling us something.

The memory wall: $0.20\times$ where we expected $\sim 1\times$.

Extrapolating the small-step tier ladder (Table 1), a naive count of floating-point work put the flexible $N=40$ near-resonance case at roughly real time; it ran at $0.20\times$. The estimate failed because it counted flops, not bytes: at dimension 81 the operator spills L1 and the kernel becomes bandwidth-bound, as quantified in Section 3.3 and Figure 4. This is not a defect to be optimised away by better vectorisation — it is a Roofline position [17] — and what recovered real time was reformulation, not tuning: MTA reduction plus segment fusion took the same case to $1.53\times$ (Figure 3).

Single precision buys nothing: the event-dive bound.

Halving the floating-point width is the reflex fix for a bandwidth-bound kernel, and we expected a corresponding speedup. We measured none worth reporting. The reason is that in the event-rich regime the wall is not matvec throughput but event localisation: the spectral event path (Alg. 1) spends 74 matvec/cycle alongside 124 scalar backlash-crossing evaluations, and the crossing detection — root bracketing on a piecewise-smooth switching surface [14,28] — must stay in double precision, or the $O\left(\Delta t\right)$ crossing accuracy (spectral-vs-brute $\sim {10}^{-3}$) degrades below the tolerance that makes the exact segment propagation meaningful. The kernel is dive-bound, not width-bound; f32 only shrinks the bytes we are not waiting on.

Modal truncation augmentation corrects $0.19\%$.

We carried MTA [21,22] expecting the residual of the high, truncated modes to matter at near-resonance. It does not: the residual flexibility measures $c_{\mathrm{res}}k=1.9\times {10}^{-3}$, a $0.19\%$ static softening. The end-to-end evidence agrees — the impact-severity metric $\delta$-RMS is $0.0167$ mm at full $N=40$, $0.0164$ at plain $N=20$, and $0.0165$ at MTA $N=20$, all within $2\%$. The negative result is the useful one: for this rattle model the high modes are dynamically inert, the backlash nonlinearity dominates the response, and MTA earns its place as insurance and as a clean justification for truncation, not as a correction that changes the answer.

Ensemble: $16.6\times$ on 56 cores, not $56\times$.

The heatmap of Figure 9 (80 operating points, 8 loads × 10 speeds, 8 mesh tables preloaded) ran in $0.66$ s parallel versus $11.0$ s serial: $16.6\times$, about $30\%$ parallel efficiency on 56 cores. The linear-speedup estimate ignored that each trajectory is the same bandwidth-bound kernel as the depth problem. Breadth and depth do not hit two different walls; they hit the same memory-bandwidth wall from two sides. This is the expected outcome for many independent low-dimensional integrations once each instance’s working set spills cache [33,35,36], and it is why the flagship ceiling map (Figure 1) draws one wall, not two.

Dynamic misalignment does not break real time.

One might expect that letting the mesh stiffness track a drifting misalignment forces a re-tabulation and forfeits real time. Measured, it does not: a five-slice $\epsilon$-grid driven by a slow misalignment drift runs at $\approx 1.37\times$, indistinguishable from the fixed-slice value, because a single slice stays cache-resident between the (only 8 per 2000) switches. The price is $N_{\epsilon }\times$ table memory, not throughput — the same memory-versus-time trade in a new guise. Section 6.8 gives the measurement in full.

Zero of eighty stable: there is no Floquet object here.

Across the 80-point flexible ensemble, the count of stable period-1 orbits was $0/80$: the flexible rattle is chaotic throughout the operating range. We had expected islands of periodic motion on which the monodromy/saltation machinery of our Floquet analysis (Figure 11) [16,41] could be applied. Their absence is not a failure of the solver but a statement about the object: with grazing-induced and vibro-impact chaos filling the range [32,44,46], classical Floquet analysis has, generically, no periodic orbit to linearise about. The one period-1 skeleton we could isolate (strong-rattle 6-DOF UPO) sits at $\left|\lambda \right|=0.972$, distance $0.028$ from a Neimark–Sacker bifurcation — neutral, on the verge of losing even that. Where no periodic object exists, the right invariant is not a multiplier but an exponent: the finite-time and largest Lyapunov exponents [9,10,31] become the language of stability, and the largest exponent of $\approx 0.55$ per rattle cycle is the physical wall of Figure 7. This is the paper’s central pivot, stated here as a negative result: past this point real time must move from resolving a trajectory to characterising structure, because the trajectory itself is no longer the answerable question.

6. Discussion and Related Work

The results of Sections 4–5 sit at the intersection of four mature literatures — piecewise-exact propagation, nonsmooth event handling, performance modelling, and predictability theory — none of which we extend. This section states precisely what is borrowed, what is new, and, most importantly, where the speed we report stops being an engineering quantity and becomes a physical one.

6.1. Piecewise-Exact Propagation and Semi-Analytical Solutions

The backbone of our fast path — integrating a linear segment by its matrix exponential and re-solving only at the mesh/backlash switches — is not new. Closed-form propagation of forced piecewise-linear oscillators dates to Shaw and Holmes [1], and the modern semi-analytical programme [25,26] constructs periodic orbits of clearance systems by stitching such analytic segments across their switching manifolds. Breakpoint-stepping exponential integration is standard practice in stiff network simulation [13], and rests on the classical exponential-integrator theory [11,12,47]. Our contribution is not the segment solution but its amortization: for a constant-contact regime the mesh cycle collapses to a single matrix–vector product ($\approx 61$ ns), and dyadic power-jumping composes ${10}^6$ cycles in 634 ns — an effective $\sim 38000\times$ over cycle-by-cycle stepping (Figure 1, Alg. 2). This is a memoization of prior art, not a new integrator; against a high-order adaptive reference it is simultaneously $\sim {10}^5\times$ faster and $\sim {10}^5\times$ more accurate ($1.2\times {10}^{-14}$ at $5\times {10}^{-4}$ ms/solve versus Vern7@${10}^{-10}$ at $1.5\times {10}^{-9}$/$68.9$ ms; Figure 6), precisely because it removes discretization error inside each segment rather than reducing it.

6.2. Nonsmooth Numerics and Event Handling

Where the trajectory actually crosses backlash faces we adopt event-driven, switch-detecting formulations rather than time-stepping through the discontinuity, following the established nonsmooth-dynamics toolbox [14,27,28]. The spectral event locator of Alg. 1 is an implementation refinement in this tradition: on the strong-rattle orbit ($\Omega =1200$, $t_{\mathrm{tooth}}=262\mu$s) it drives the per-cycle cost from 371 matvec (V1, $5\times$ real time) to 74 matvec plus 124 scalar root evaluations (V3, $8\times$), and on the sparse orbit ($\Omega =1400$) reaches $29.1$ matvec / $14.6\mu$s per cycle ($16\times$) while cross-checking against brute-force integration at $7\times {10}^{-15}$ (Table 2). The linearization of the resulting hybrid Poincaré map uses saltation matrices [16] and piecewise-smooth Floquet theory [41]; the broader nonsmooth-bifurcation setting is that of Leine and Nijmeijer [15] and di Bernardo et al. [42]. We claim no new theory here — only that the analytic segment structure makes each of these standard operations cheap.

6.3. Performance Modelling: The Roofline as an Interpretive Frame

We read all of our timings through the Roofline model [17], exactly as it is used to diagnose memory-bound stencil and finite-difference kernels [30,36]. This framing is what lets us state the paper’s central empirical claim without hand-waving: single-trajectory depth and ensemble breadth strike the same memory-bandwidth wall, with the depth-side cache cliffs and the breadth-side bandwidth saturation measured directly in Section 3.4 (Figure 4, Figure 5). Large independent-ODE ensembles on GPUs [33,34,35] would raise the absolute throughput but not the arithmetic intensity of a matvec-bound segment, so the same ceiling reappears one level up.

6.4. Model Reduction and the Flexible Tier

The flexible tiers use standard component-mode and truncation-augmentation machinery: mode-acceleration and modal truncation augmentation [21], dual Craig–Bampton synthesis [24], and residual-attachment enrichment [22], in the spirit of finite-element geared-rotor analysis [23]. Data-driven reduction of piecewise-linear structures via DMD [29] is a complementary route we do not take. Our only quantitative addition is to price reduction against the same memory wall: MTA retains super-real-time behaviour ($N20$ MTA+segment-fusion reaches $1.53\times$, $171\mu$s, 52 matvec) at a residual-flexibility softening of $c_{\mathrm{res}}k=1.9\times {10}^{-3}$ (0.19%), with $\delta$-RMS agreeing across full-$N40$/simple-$N20$/MTA-$N20$ to within 2% ($0.0167$/$0.0164$/$0.0165$ mm; Figure 3, Alg. 3). Reduction here buys cache residency, not accuracy.

6.5. Continuation, Grazing, and Structural Real Time

Once trajectory error is dominated by sensitivity rather than discretization, the useful objects are invariant structures, and these too are prior art: numerical continuation [38,39], multiple-shooting Newton–Krylov location of periodic orbits [37], and the grazing-bifurcation theory of vibro-impact systems [40,44,45,46]. Our Floquet analysis of the 6-DOF strong-rattle period-1 skeleton (Newton shooting) returns a neutral pair at 1 and a dominant complex multiplier $-0.373-0.897\mathrm{i}$, $\left|\lambda \right|=0.972$, at distance $0.028$ from the unit circle — a near Neimark–Sacker configuration — while the sparse period-2 orbit is unstable (${\left|\lambda \right|}_{max}=1.92$) (Figure 11). Computing this skeleton is cheap precisely because the segment propagator and its saltation-corrected monodromy are analytic; the point is methodological, that at high $\Omega$ the real-time budget is better spent continuing structure than propagating any one trajectory.

6.6. Real-Time Multibody and HIL Context

Framed as a real-time-capable multibody problem [8,48], our tiers occupy the same design space as HIL powertrain gearshift models [7], with the fidelity-versus-real-time trade made explicit by the tier ladder of Figure 2 and Table 1, and the physical output is the genuine NVH signature of Section 4.4 (Figure 12). This grounds the model in the classical gear-rattle and gearbox-chaos literature [2,3,4,5,6,20,32] with a physically calibrated TVMS [18,19].

6.7. Engineering Ceiling Versus Physical Ceiling

The above lets us separate two kinds of limit that are routinely conflated. Everything in Section 6.1, Section 6.2, Section 6.3, Section 6.4, Section 6.5 and Section 6.6 is an engineering ceiling: it is set by discretization, cache size, arithmetic intensity, and bandwidth, and it is therefore optimizable — the $0.20\times \to 1.53\times$ chain (Figure 3) is exactly the recovery of a super-real-time budget by moving work below the memory wall. The Lyapunov ceiling is different in kind (Figure 7). The flexible rattle is all-chaotic across the operating range studied — $0/80$ atlas points carry a stable period-1 orbit — and the largest exponent is $\approx 0.55$ per rattle cycle (state separation $4.57\times {10}^{-4}$ at 3 cycles growing to $2.09\times {10}^{-2}$ at 10 cycles). No integrator removes this: it is a property of the flow, not of the numerics, and the finite-time predictability literature [9,10,31] states the consequence exactly — exponential compute buys only a linear extension of the predictable horizon. This is the boundary the title names. Below it, faster propagation returns proportionally more usable trajectory; at it, per-trajectory real time is physically defeated and the real-time target must shift from the trajectory to the invariant structure (periodic-orbit skeleton, multipliers, basins), which is what Section 6.5 computes. We do not claim to beat the Lyapunov wall; we claim to locate it and to show that structure remains real-time-tractable on its far side.

6.8. Limitations

Several limitations bound these claims, and we report them as such rather than smoothing them over. The segment-fusion and spectral-table machinery assumes the operating regime is well described by a small set of linear pieces; where regime dwell times shorten the event rate rises to $2.7$–$2.9$ per cycle (Table 4) and the amortization advantage erodes toward the raw matvec cost. The spectral tables trade memory for time — 18 MB in the most complex flexible+backlash+friction configuration, versus 190 MB for a dyadic layout — so the depth speedups are contingent on this table remaining cache- and memory-resident, precisely the wall of Section 6.3. Accuracy of the reduced flexible tiers is validated only in an $O\left(\mathrm{d}t\right)$ sense (spectral-vs-brute $\sim {10}^{-3}$ for the full backlash configuration, $5.06\times {10}^{-6}$ near resonance; reconstruction $1.5\times {10}^{-11}$; conditioning cond$\left(V_1\right)$ ranging $8.1\times {10}^5$ free to $\sim {10}^3$ in contact), and the friction model is a fixed Coulomb $\mu =0.1$ giving a TE-RMS within 5.8% of the small-step reference ($5.89$ vs $6.25\mu$m); richer tribology, temperature, and lubrication are out of scope. A further modelling choice bounds the most complex tier. The deformation-compatible stiffness engine resolves the full misalignment dependence $k_m(\varphi ,{\epsilon }_{\mathrm{eff}},T)$ — a 1 mrad misalignment unloads part of the flank and drops the tangent stiffness by $42\%$, and the 12-DOF Tier 3 runs this as a closed loop (Figure 2) — but the flexible Tier 5 spectral tables reported here are built at the nominal ($\epsilon =0$) slice, so misalignment-varying stiffness and flexible-modal backlash are not yet simultaneously active in a single run. This is a wiring gap, not a barrier, and importantly not a throughput one — which we verified directly. A fixed misalignment is a free re-tabulation, since the segment count and propagator structure are unchanged and the runtime is therefore identical. For the dynamic case we pre-tabulated a five-slice $\epsilon$-grid ($\epsilon \in [-1,+1]$ mrad, at the light $T=50$ N m load where misalignment drops the mesh stiffness by up to $42\%$) and drove it with a slowly drifting ${\epsilon }_{\mathrm{eff}}$, selecting the nearest slice each cycle. Because ${\epsilon }_{\mathrm{eff}}$ drifts over many mesh cycles, the run switched tables only 8 times over 2000 cycles (one slice stays cache-resident for $\approx 250$ cycles between switches), and the slice-switching kernel ran at $\approx 1.35\times$ real time — statistically indistinguishable from the fixed-slice $\approx 1.37\times$, the difference lying within run-to-run scatter — while the response changed measurably ($\delta$-RMS $+4.3\%$). The cost of the dynamic case is thus $N_{\epsilon }\times$ table memory ($5\times 4.8$ MB for the reduced tier here), which pushes on the same Roofline wall of Section 6.3 once the slice set spills L3, rather than lost single-trajectory real time. The ensemble efficiency of $\sim 30\%$ is bandwidth-bound on a single 56-core node and we make no distributed- or GPU-scaling claims. Finally, the Floquet result is a local statement about one period-1 skeleton at one operating point; a full unfolding of the near-Neimark–Sacker structure across the atlas, and validation against experiment, remain open. The negative results (Table 5) — no stable period-1 orbit in range, sub-real-time performance at dim $\approx 81$, and sub-linear ensemble scaling — are as much a part of the contribution as the speedups.

7. Conclusion

We set out to answer two questions about high-fidelity gear rattle: how fast can it be simulated, and where does simulating single trajectories stop being the right thing to do. The two answers turn out to be tightly coupled, and neither rests on a new integrator or a new framework.

How fast. For the nonsmooth gear model in its constant-contact regime, exact piecewise propagation collapses one mesh cycle to a single matrix–vector product, and dyadic power-jumping advances the state across a million cycles in the time a naive solver spends on a handful. This is the source of the super-real-time factors reported along the tier ladder (Table 1) and the layered ceiling of Figure 1. But the ceiling is not open sky. Once the impact-driven event path dominates, or the modal dimension grows past the point where the propagator matrix leaves cache, the throughput saturates: the same memory-bandwidth wall that limits a single deep trajectory also limits the breadth of a many-core ensemble, where the parameter heatmap runs at roughly a third of ideal core-scaling rather than near-linear (Figure 4). Depth and breadth are two views of one Roofline-type constraint [17,30,36]. Super-real-time is achievable and useful, but it is bounded, and we have tried to report the bound rather than the peak.

Where trajectories give way to structure. The second wall is physical, not computational. With a largest Lyapunov exponent of order one per rattle cycle, an exponential increase in compute buys only a linear increase in the horizon over which a computed trajectory tracks the true one [9,10,31]. Figure 7 makes this concrete: past a few rattle cycles the individual orbit is no longer the meaningful object. The right real-time deliverable there is structure that survives the loss of the orbit – invariant sets, Floquet stability of the underlying periodic skeleton (Figure 11), and the Lyapunov atlas over operating conditions (Figure 9), which in the range studied is uniformly chaotic. The speed of the propagator is what makes computing these structural objects in real time feasible; the Lyapunov wall is what makes computing them the right goal.

Within this framing our contributions are modest and specific. First, we assemble known ingredients – exact propagation of piecewise-linear dynamics, event localization, modal truncation augmentation – into a high-fidelity gear rattle simulation that runs super-real-time (Table 3, Figure 3). Second, we show that single-trajectory depth and ensemble breadth are limited by the same memory-bandwidth wall. Third, we argue and illustrate that at the Lyapunov wall the real-time target must move from trajectories to structure. Fourth, we report the negative results honestly: the configurations where the fusion advantage evaporates, the flexible cases that do not clear real time, and the absence of stable period-one motion across the swept range (Table 5). None of the underlying methods – exact piecewise propagation, saltation and Floquet analysis, shooting and continuation, the Roofline model, Lyapunov-based predictability – is ours; the increment is in composing them for this problem and in mapping where each stops helping.

Two directions follow naturally. Nonsmooth continuation of periodic orbits, including grazing branches, would replace the sampled Lyapunov atlas with a cleaner bifurcation atlas of the rattle skeleton [38,39,40]; and porting the memory-bound ensemble to a GPU many-ODE solver [33,34,35] would push breadth further along – though, if our reading of the wall is right, toward the same bandwidth limit rather than through it.

Appendix A.1. Augmented Per-Segment LTI State

Within a single contact state $s\in \{\mathrm{free},\mathrm{drive},\mathrm{coast}\}$ the reduced dynamics are linear time-invariant with periodic excitation. We fold the static-transmission-error (STE) forcing into the state so that each segment is autonomous. For the flexible N-modal configuration the physical state $(\mathbf{q},\dot{\mathbf{q}})\in {\mathbb{R}}^{2N}$ is augmented by the accumulated relative rotation (the mesh phase), giving the augmented dimension $dim=2N+1$. Writing $\mathbf{z}={(\mathbf{q},\dot{\mathbf{q}},\varphi )}^{\top }$, the segment obeys

$$\dot{\mathbf{z}}=A_s\mathbf{z}+{\mathbf{b}}_s,A_s=\left[\begin{array}{ccc}0& I& 0\\ -M^{-1}{K}_s& -M^{-1}{C}_s& M^{-1}{\mathbf{f}}_s\\ 0& 0& 0\end{array}\right],$$

(A1)

where the trailing zero row makes the mesh phase, advancing at constant rate, a neutral (integrator) coordinate rather than an external clock. The single harmonic STE excitation ($h_1=5\mu$m, the mesh-fundamental amplitude of Section 4.4) enters through ${\mathbf{f}}_s$; higher harmonics are appended as further neutral oscillator pairs when required. Absorbing the drift ${\mathbf{b}}_s$ by one additional augmented coordinate yields a homogeneous $\dot{\mathbf{z}}={\tilde{A}}_s\mathbf{z}$, so a single segment step is one action of the propagator ${\Phi }_s\left(\tau \right)=exp\left({\tilde{A}}_s\tau \right)$.

Appendix A.2. Schur Ordering and Sylvester Decoupling

Forming $exp\left({\tilde{A}}_s\tau \right)$ directly is poorly conditioned: ${\tilde{A}}_s$ carries a defective eigenvalue at the origin (the rigid-rotation Jordan block coupled to the phase and carry coordinates) alongside lightly damped oscillatory modes whose eigenvectors are nearly parallel near resonance. We therefore never diagonalize the full operator. A real Schur factorization ${\tilde{A}}_s=QT{Q}^{\top }$ is reordered so that the neutral subspace — the rigid-rotation Jordan block together with the constant-drift and mesh-phase carry coordinates — is collected in the leading diagonal block $T_{11}$, while the remaining lightly damped modes appear as $2\times 2$ real blocks along $T_{22}$:

$$T=\left[\begin{array}{cc}{T}_{11}& T_{12}\\ 0& T_{22}\end{array}\right],T_{22}=blockdiag\left(R_1,\dots ,R_m\right),R_j=\left[\begin{array}{cc}-{\zeta }_j{\omega }_j& {\omega }_{d,j}\\ -{\omega }_{d,j}& -{\zeta }_j{\omega }_j\end{array}\right].$$

(A2)

The off-diagonal coupling $T_{12}$ is annihilated by solving the Sylvester equation $T_{11}X-X{T}_{22}=-T_{12}$ for the block X; the similarity $blockdiag(I,·)$ with X in the $(1,2)$ position block-diagonalizes T exactly. Because $T_{11}$ (spectrum at the origin) and $T_{22}$ (spectrum in the open left half-plane) have disjoint spectra, the Sylvester solve is well posed, and the neutral subspace is isolated from the oscillatory blocks without ever inverting the near-degenerate modal eigenvector matrix. This is what keeps the change of basis benign: the probed conditioning of the segment transform $V_1$ is $\kappa \left(V_1\right)\approx 8.1\times {10}^5$ in the free (backlash) segment — where the two flanks are uncoupled and the oscillatory modes crowd together — but drops to $\kappa \left(V_1\right)\approx {10}^3$ in the contact segments. Even the worst-conditioned free block reconstructs the propagator to a relative residual of $1.5\times {10}^{-11}$, well below every physical tolerance in the study.

Appendix A.3. Closed-Form Block Exponential

On the decoupled form the exponential is available in closed form block by block, so no scaling-and-squaring is needed. Each oscillatory $2\times 2$ block integrates to a damped rotation,

$$exp\left(R_j\tau \right)=e^{-{\zeta }_j{\omega }_j\tau }\left[\begin{array}{cc}cos{\omega }_{d,j}\tau & sin{\omega }_{d,j}\tau \\ -sin{\omega }_{d,j}\tau & cos{\omega }_{d,j}\tau \end{array}\right],$$

(A3)

while the neutral block, a nilpotent-plus-carry Jordan form, exponentiates to a low-order polynomial in $\tau$ (degree at most two for the rigid-rotation-plus-drift coupling). The segment displacement therefore has the exact separated form

$$\delta \left(\tau \right)=\underset{\mathrm{damped}\sin \mathrm{usoids}}{\underbrace{{\sum }_j{e}^{-{\zeta }_j{\omega }_j\tau }\left({\mathbf{a}}_{j}cos{\omega }_{d,j}\tau +{\mathbf{b}}_{j}sin{\omega }_{d,j}\tau \right)}}+\underset{\mathrm{low}-\mathrm{order}\mathrm{polynomial}}{\underbrace{{\mathbf{c}}_0+{\mathbf{c}}_1\tau +{\mathbf{c}}_2{\tau }^2}},$$

(A4)

with coefficient vectors fixed by the initial state of the segment. Equation (A4) is evaluated at machine cost and is the object monitored for a flank crossing; the same coefficients give $\dot{\delta }\left(\tau \right)$ analytically for the Newton step below.

Appendix A.4. Scalar-Newton Event Location

A contact-state switch occurs when a scalar gap function $g\left(\tau \right)={\mathbf{e}}^{\top }\delta \left(\tau \right)-{\delta }_{\pm }$ (relative flank displacement minus the half-backlash) crosses zero. Because (A4) is closed form, g and $g^{\prime }$ are evaluated without integration, and the crossing time ${\tau }^{★}$ is found by scalar Newton iteration seeded from the sign-change bracket, with a bisection safeguard. This is the per-segment inner loop of Algorithm 1: propagate exactly to a candidate horizon, detect the earliest sign change of g, refine ${\tau }^{★}$ by scalar Newton, apply the switch, and repeat. Each refinement is a handful of scalar evaluations of (A4) rather than matrix work, which is why the spectral variant spends only 74 matvec per cycle plus 124 scalar evaluations per cycle at the strong-rattle operating point, with the flexible configuration exhibiting $2.7$–$2.9$ switches per mesh cycle.

At the located switch the flow is continuous but its Jacobian is not; the sensitivity of the trajectory across the transition is carried by the saltation matrix, which maps the pre-event to the post-event variation and is the correct linearization for the hybrid flow [15,16]. The same saltation factors assemble the monodromy matrix used for the Floquet analysis of Section 4.3; here they simply confirm that the exact propagator composed across events preserves the neutral pair at unity expected of the rigid-rotation mode.

Algorithm 1:Spectral closed-form event location (V3) for one flow segment: locating the event is scalar work (0 matvec) and only applying the exact jump costs 2 matvec, cutting the strong-rattle event path from 371 matvec/cycle (V1 dyadic subdivision) to 74.

Input: state $x\in {\mathbb{R}}^n$, segment eigenbasis W (precomputed Schur + Sylvester decouple), block generator, contact-state s, remaining time $t_{\mathrm{rem}}$ Output: event time $t^{★}$, jumped state x, new contact-state 1: $y\leftarrow W^{-1}x$▹matvec 1: to spectral coords 1: – everything below is scalar / vector on the decoupled blocks:0 matvec– 2: $\delta \left(t\right)\leftarrow g_w·blockexp(y,t)-e\left[s\right]$▹ closed form: damped sinusoids + low-order poly 3: ${\delta }_0\leftarrow \delta \left(0\right)$;   $B\leftarrow {max}_{t\in [0,t_{\mathrm{rem}}]}\left|\dot{\delta }\left(t\right)\right|$▹ analytic derivative bound 4: if$margin\left({\delta }_0\right)>t_{\mathrm{rem}}·B$then▹ prescreen: no crossing possible 5:     $z\leftarrow blockexp(y,t_{\mathrm{rem}})$;   $x\leftarrow Wz$ ▹ matvec 2: advance full segment, skip 6:     return $(t^{★}=\infty ,x,s)$ 7: end if 8: $t^{★}\leftarrow \infty$ 9: for each monotone sub-interval $[a,b]$ of $\delta$ (from block frequencies) do 10:     if $\delta \left(a\right)\delta \left(b\right)\le 0$ then 11:         $\tau \leftarrow SafeNewton\left(\delta ,\dot{\delta },[a,b]\right)$ ▹ scalar; bracketed, quadratic 12:         $t^{★}\leftarrow min(t^{★},\tau )$ 13:     end if 14: end for 15: $t^{★}\leftarrow \Delta t⌈t^{★}/\Delta t⌉$▹ quantize to unit time grid 16: $z\leftarrow blockexp(y,t^{★})$▹ closed-form advance, 0 matvec 17: $x\leftarrow Wz$▹matvec 2: back to physical coords 18: $s^{\prime }\leftarrow classify\left(\delta \left(t^{★}\right),\dot{\delta }\left(t^{★}\right)\right)$▹ engage / release / reverse 19: return$(t^{★},x,s^{\prime })$

V1: dyadic descent plus bisection.

The first working locator treated each candidate crossing as an opaque sign change and isolated it by dyadic descent followed by bisection on the guard function. Each guard evaluation at a trial time requires re-propagating the state to that time, i.e. one matvec, so the geometric refinement of the bracket accumulates a logarithmic number of matvec per event. Aggregated over a full mesh cycle this reaches 371 matvec/cycle, yielding only a $5\times$ real-time factor. The scheme is robust—it inherits the guaranteed sign-change convergence of bisection—but blind to the structure of the flow it is bracketing.

V2: certificate plus Hermite guidance.

V2 replaces blind bracketing with a two-part strategy. A cheap monotonicity/convexity certificate on the guard over the current segment rules out spurious crossings without any propagation, and where a crossing is certified a Hermite interpolant built from the guard value and its time derivative at the segment endpoints supplies a high-order first guess, so that at most one or two corrective matvec suffice to reach tolerance. This halves the work to 193 matvec/cycle and lifts the factor to $8\times$. The residual cost is structural: each corrective step still propagates the full state vector to test the guard.

V3: closed-form spectral event location.

The final locator exploits the fact that, within a segment, the exact flow is a spectral (eigenmode) superposition. The guard is then a fixed linear combination of scalar modal exponentials, and its zero can be sought directly in the scalar spectral representation—no state propagation, no matvec—until the crossing time is known, at which point a single matvec advances the state across the boundary. Locating is thus reduced to scalar root-finding on a closed-form expression, approaching the bound of one scalar solve per event. On the strong-rattle point V3 spends 74 matvec/cycle plus 124 scalar guard evaluations (still $8\times$, now locator-bound rather than location-bound); on the sparse point ($\Omega =1400$) it needs only $29.1$ matvec and $14.6\mu \mathrm{s}$ per cycle for a $16\times$ factor, and the crossing times agree with a brute-force reference to $7\times {10}^{-15}$. The scalar evaluations are effectively free relative to matvec, so V3’s cost is set by the number of segments that must each be crossed with one matvec, not by the search within them.

Segment fusion.

With location made effectively free, the residual matvec count equals the number of segments, and adjacent segments that share the same active contact state can be fused into a single exact propagation. Algorithm 2 composes the per-segment spectral maps over a maximal run of like-state segments and applies the composite in one matvec, eliminating the boundary matvec at every fused interface. On the near-resonance flexible configuration this brings the residual count down to 52 matvec, and the fused trajectory reproduces the segment-by-segment spectral result to $4.0\times {10}^{-15}$—machine-precision agreement, confirming that fusion is an exact regrouping rather than an approximation. This is the step that carries the most complex configuration across the real-time line (Figure 3, $171\mu \mathrm{s}$, $1.53\times$).

Appendix C.1. Residual Flexibility

Let $\{{\omega }_i,{\varphi }_i\}$ be the free-interface normal modes and $w_{c,i}={\varphi }_i^{\top }{\mathbf{b}}_c$ the modal participation of the mesh-contact load vector ${\mathbf{b}}_c$. Retaining the first $N_k$ modes and treating all higher modes as massless—their inertia being negligible at the mesh excitation frequency—each dropped mode contributes a static compliance $w_{c,i}^2/{\omega }_i^2$ at the contact coordinate. The summed high-mode residual flexibility is

$$c_{\mathrm{res}}=\sum _{i>N_k}{\displaystyle \frac{w_{c,i}^2}{{\omega }_i^2}},$$

(A5)

the standard mode-acceleration / modal-truncation-augmentation correction [21]; it is the projection onto the mesh degree of freedom of the higher-order residual attachment modes used in free-interface component-mode synthesis [22,24]. The truncated model sees the contact spring k in series with this residual compliance, so the effective mesh stiffness felt by the retained subspace is

$$k_{\mathrm{eff}}={\displaystyle \frac{k}{1+c_{\mathrm{res}}k}}.$$

(A6)

Evaluating $c_{\mathrm{res}}$ for the $N_k=20$ truncation gives $c_{\mathrm{res}}k=1.9\times {10}^{-3}$, i.e. a $0.19\%$ softening of the mesh stiffness. Because the correction enters $k_{\mathrm{eff}}$ statically rather than through additional states, the reduced operator stays at $dim41$ and the per-step matvec cost is identical to the uncorrected $N=20$ model.

Appendix C.2. Fidelity

The augmentation is a scalar stiffness correction, yet it recovers essentially all of the fidelity lost to truncation. On the near-resonance transmission-error trajectory, the displacement RMS deviation against a common high-fidelity reference is $\delta -\mathrm{RMS}=0.0165$ mm for the MTA-corrected $N=20$ model, statistically indistinguishable from $0.0167$ mm for the full $N=40$ ($dim81$) model and $0.0164$ mm for the uncorrected simple $N=20$ truncation—the three agree to better than $2\%$. The $0.19\%$ stiffness shift therefore sits below the residual chaotic dispersion of the trajectories themselves and does not degrade the solution, while enabling the $dim81\to dim41$ reduction that, combined with segment fusion, carries the near-resonance case from $0.20\times$ to $1.53\times$ real time (Figure 3). The complete pre-processing is summarised in Alg. 3; its output $(k_{\mathrm{eff}},{\{{\omega }_i,{\varphi }_i\}}_{i\le N_k})$ feeds directly into the spectral event detector (Alg. 1) and the fused propagator (Alg. 2).

Algorithm 2:Segment-run fusion fast cycle: certified event-free runs skip r segments in one matvec, while event-bearing segments fall back to the spectral single-segment path, cutting mean cost from $124\to 52$ matvecs per cycle.

Precompute (once per parameter set), for each contact-state c and level $r\in \{1,2,4\}$:    ${\displaystyle {\Phi }_{c,r}^{\mathrm{fused}}\leftarrow {\displaystyle \prod _{s=0}^{r-1}}{\Phi }_{c,s}}$,   ${\Phi }_{c,s}=W\mathrm{blockexp}\left(\Delta t_s\right)W^{-1}$▹r consecutive propagators collapsed to one matvec 1: $\mathrm{pos}\leftarrow 0$;   $x\leftarrow x_0$▹ state at start of the segment grid 2: while$\mathrm{pos}<N$do 3:     $({\delta }_0,{\dot{\delta }}_0)\leftarrow (e_{\delta }^{\top }x,e_{\dot{\delta }}^{\top }x)$;   $c\leftarrow ContactState\left({\delta }_0\right)$▹ gap and gap-rate at $\mathrm{pos}$ 4:     $m_0\leftarrow Margin({\delta }_0,c)$▹ signed distance to the nearest engagement boundary 5:     $r\leftarrow LargestAlignedRun(\mathrm{pos},\{4,2,1\left\}\right)$▹ try coarsest level first, step down 6:     while $r\ge 1$ do 7:         $\Delta t_{\mathrm{run}}\leftarrow {\sum }_{s=0}^{r-1}\Delta t_{\mathrm{pos}+s}$ 8:         if $m_0>\kappa \Delta t_{\mathrm{run}}\left|{\dot{\delta }}_0\right|+e_{\mathrm{span}}$then▹prescreen: linear reach + curvature slack cannot cross the wall 9:            $y\leftarrow {\Phi }_{c,r}^{\mathrm{fused}}x$▹one matvec advances r segments 10:            $({\delta }_r,m_r)\leftarrow (e_{\delta }^{\top }y,Margin(e_{\delta }^{\top }y,c))$ 11:            if $ContactState\left({\delta }_r\right)=c\wedge min(m_0,m_r)>g_{\mathrm{guard}}$then▹ endpoint certificate: same state, both margins clear 12:                $x\leftarrow y$;   $\mathrm{pos}\leftarrow \mathrm{pos}+r$▹ certified event-free run accepted 13:                goto next 14:            end if 15:         end if 16:         $r\leftarrow r/2$▹ certificate failed: shorten the run and retry 17:     end while 18:     $x\leftarrow SpectralEvent(x,\mathrm{pos})$▹ fallback: resolve the event-bearing segment (Alg. 1) 19:     $\mathrm{pos}\leftarrow \mathrm{pos}+1$ 20:     next: 21: end while

Algorithm 3:Modal-truncation augmentation (MTA): folding the truncated high modes in quasi-statically as a residual-flexibility softening $k\to k/(1+c_{\mathrm{res}}k)$ drops DOF into cache without losing static compliance, at a measured mesh-stiffness bias of only $0.19$ ($c_{\mathrm{res}}k=1.9e-3$).

Require: full modal basis of N modes: frequencies ${\left\{{\omega }_i\right\}}_{i=1}^N$; per-segment mesh projections $\left\{w_{c,i}\left(s\right)\right\}$; linear mesh stiffness $k\left(s\right)$; retained count $N_k<N$ Ensure: reduced dynamic system of dimension $2N_k+1$ per segment, statically exact at the mesh Ensure: 1: keep the $N_k$ lowest modes as dynamic DOF ▹ these stay in the time integrator 2: for each mesh segment s do 3:     ${\displaystyle c_{\mathrm{res}}\left(s\right)\leftarrow \sum _{i>N_k}\frac{w_{c,i}{\left(s\right)}^2}{{\omega }_i^2}}$▹ quasi-static tail of truncated modes 4:     ${\displaystyle k_{\mathrm{eff}}\left(s\right)\leftarrow \frac{k\left(s\right)}{1+c_{\mathrm{res}}\left(s\right)k\left(s\right)}}$▹softening: never stiffer than k 5:     assemble reduced segment matrix ${\mathbf{A}}_s\in {\mathbb{R}}^{(2N_k+1)\times (2N_k+1)}$ using $k_{\mathrm{eff}}\left(s\right)$ 6: end for 7: return$\left\{{\mathbf{A}}_s\right\}$▹dimension collapses; high-mode static flexibility retained 7: 8: Measured: residual correction $c_{\mathrm{res}}k=1.9e-3$ ($0.19$) ▹ truncation bias below noise floor

Newton shooting on the cycle map.

A periodic orbit of period k mesh cycles is a fixed point of the cycle map $P\left(x\right)$, evaluated by the exact segment propagator of Appendix A. Newton’s update solves $(I-M)\Delta x=P\left(x\right)-x$, where the monodromy $M=\mathrm{d}P/\mathrm{d}x$ is accumulated analytically along the orbit: within each affine segment the fundamental solution is the same $exp\left(A_{\sigma }\Delta t\right)$ that propagates the state, and at each backlash crossing the chain is bridged by the saltation matrix [15,16]. For the present model the vector field is continuous across the switching surface up to the mesh-damping term, so the saltation reduces to $S=I\mp c_m(\mathrm{inv}m\odot w)w^{\mathsf{T}}$; its numerical effect on the multipliers is below ${10}^{-5}$ here, but retaining it is what makes the monodromy exact for the sequence and measurably improves Newton convergence.

Figure A1. End-to-end pipeline. Offline (grey), a single deformation-compatible meshing model is decoupled per (segment, contact-state) into an 18$\mathrm{M}$$\mathrm{B}$ spectral table (GT5B). Online (amber), each mesh cycle sweeps its 60 segments through a four-step loop — project into spectral coordinates ($W^{-1}x$, one matvec), test in closed form whether the deflection $\delta$ reaches the backlash edge $\pm b$, then either block-exponential jump the whole segment or locate the crossing time $t^{★}$ and switch contact state, and finally map back ($Wz$, one matvec) — so a raw cycle costs 120 matvec, cut to 52 by segment fusion. The resulting super-real-time state feeds the 56-core branch (teal) that extracts bifurcation, chaos, Floquet, and NVH structure.

Figure A1. End-to-end pipeline. Offline (grey), a single deformation-compatible meshing model is decoupled per (segment, contact-state) into an 18$\mathrm{M}$$\mathrm{B}$ spectral table (GT5B). Online (amber), each mesh cycle sweeps its 60 segments through a four-step loop — project into spectral coordinates ($W^{-1}x$, one matvec), test in closed form whether the deflection $\delta$ reaches the backlash edge $\pm b$, then either block-exponential jump the whole segment or locate the crossing time $t^{★}$ and switch contact state, and finally map back ($Wz$, one matvec) — so a raw cycle costs 120 matvec, cut to 52 by segment fusion. The resulting super-real-time state feeds the 56-core branch (teal) that extracts bifurcation, chaos, Floquet, and NVH structure.

The neutral direction.

The gear pair is rotationally symmetric: the direction $n_r$ that advances both wheels by a common rigid rotation leaves the mesh deflection unchanged ($w·n_r=0$, hence $M{n}_r=n_r$), so $(I-M)$ is structurally rank-deficient — and, physically, $P\left(x\right)-x$ has a genuine nonzero component along $n_r$, the nominal per-cycle advance. Periodicity must therefore be imposed only on the complement: we regularize with $J_{\mathrm{reg}}=(I-M)+\sigma \widehat{n}{\widehat{n}}^{\mathsf{T}}$ and project the $\widehat{n}$-component out of the update. With this, Newton contracts cleanly even onto unstable orbits that no forward simulation could ever exhibit: the strong-rattle period-1 orbit converges from residual $2.1\times {10}^2$ to $1.9\times {10}^{-9}$ within ten Newton iterations, a four-event rattling orbit with $\delta$ spanning both flanks. In the sparse-rattle regime the period-1 solve does not converge but the period-2 solve does — direct evidence that a flip (period-doubling) has already occurred there, consistent with its multiplier ${\left|\lambda \right|}_{max}=1.92$.

Ground truth versus Floquet: the grazing caveat.

Because shooting delivers the orbit and its monodromy independently, the linearized prediction can be tested against the true map. Kicking each converged orbit by ${10}^{-5}$ (in reduced norm) and iterating the full nonlinear map gives per-return growth factors of $2.61$ for the sparse period-2 orbit — in reasonable accord with its unstable multiplier — but $1.51$ for the strong-rattle period-1 orbit whose dominant multiplier is $\left|\lambda \right|=0.972$, i.e. nominally stable. The first return already amplifies the kick from ${10}^{-5}$ to $6.3\times {10}^{-2}$: the perturbation does not evolve along the linearized flow at all, it changes which backlash crossings occur, and the orbit jumps to a different event sequence. This grazing-type sensitivity [40,44,45,46] is invisible to any monodromy computed along the unperturbed sequence, however carefully the saltation is handled. The practical rule we draw: in a nonsmooth system, Floquet multipliers certify stability only within the margin by which the event sequence itself is robust; near grazing, the Lyapunov exponent measured on the true map is the only honest sensitivity.

Design sensitivity: linearize per segment, not per cycle.

The same lesson governs parameter derivatives. A Fréchet derivative of the cycle map with respect to a design parameter (a stiffness scale, a misalignment) is unusable in the chaotic regime: the phase drift it must represent grows to $O\left(1\right)$ within a cycle and the linearization leaves its validity ball. Differentiating at the segment level and re-fusing the perturbed propagators, by contrast, is accurate to $O\left(\Delta {\epsilon }^2\right)$ — the segment is short enough for the linearization to hold, and the fusion recomposes exactly. Design sensitivity in this framework therefore inherits the same architecture as the solver itself: everything exact happens per segment; cycles are only ever products.

Hardware.

All timings were measured on a dual-socket workstation carrying two Intel Xeon E5-2680 v4 processors (Broadwell-EP, 14 cores each) for 28 physical cores and 56 hardware threads in total, at $2.4$ GHz base ($3.3$ GHz turbo), organised as two NUMA nodes. The per-core cache hierarchy is 32 KB L1 data, 32 KB L1 instruction, and 256 KB L2, with a 35 MB L3 shared per socket; main memory is 32 GB of DDR4 (WSL2 exposes $\approx 26$ GB to the runtime). The cores implement AVX2 and FMA but not AVX-512. Two features of this machine are load-bearing for the argument. First, the 32 KB L1 data cache is exactly the threshold identified in Section 3.3: the dim-81 system operator ($81\times 81\times 8=52$ KB) overflows it, which is why the matvec cost jumps at that dimension and why the fidelity knee sits near DOF $\approx 22$. Second, the dual-socket NUMA topology sharpens the depth-versus-breadth wall of Figure 4: the 56-thread ensemble does not see 56 independent memory paths but two socket-level memory controllers, so independent trajectories contend for a shared and partly cross-socket bandwidth — the mechanism behind the measured $16.6\times$ (rather than $56\times$) ensemble speedup and its $\sim 30\%$ parallel efficiency. Figure A2 draws the hierarchy explicitly and places every working set of the solver on it: the reduced dim-41 operator (13 KB) inside each core’s L1, the raw dim-81 operator (52 KB) spilling into L2, the N40 spectral table (18 MB) resident in one socket’s L3, and the 8-table ensemble (144 MB) in DRAM behind the measured ≈32 GB/s aggregate-bandwidth ceiling of Figure 5b.

Figure A2. Where the data lives and how it moves: the dual-socket memory hierarchy with the solver’s working sets mapped onto it. The reduced dim-41 operator (13 KB) is L1-resident per core; the raw dim-81 operator (52 KB) overflows L1 into L2; the full N40 spectral table (18 MB) is L3-resident on one socket; and the 8-table ensemble (144 MB) spills to DRAM, whose measured bandwidth (Fig. Figure 5b) saturates at ≈32 GB/s for the whole machine. Each matvec streams one operator column-wise through the FMA units, so a raw-N40 mesh cycle moves ≈6.3 MB of table — the kernel is a bandwidth workload, and every optimization in this paper (MTA reduction, segment fusion, spectral tables) is at bottom a scheme to move this traffic up the hierarchy.

Figure A2. Where the data lives and how it moves: the dual-socket memory hierarchy with the solver’s working sets mapped onto it. The reduced dim-41 operator (13 KB) is L1-resident per core; the raw dim-81 operator (52 KB) overflows L1 into L2; the full N40 spectral table (18 MB) is L3-resident on one socket; and the 8-table ensemble (144 MB) spills to DRAM, whose measured bandwidth (Fig. Figure 5b) saturates at ≈32 GB/s for the whole machine. Each matvec streams one operator column-wise through the FMA units, so a raw-N40 mesh cycle moves ≈6.3 MB of table — the kernel is a bandwidth workload, and every optimization in this paper (MTA reduction, segment fusion, spectral tables) is at bottom a scheme to move this traffic up the hierarchy.

Software.

The runtime kernel is written in Rust (rustc/cargo $1.96.1$) and compiled with the release profile lto = true, codegen-units = 1, panic = "abort", and RUSTFLAGS = "-C target-cpu=native" so the segment matvec is vectorised to the machine’s AVX2/FMA width. Two low-level details matter for the reported numbers. Denormal floating-point operands — which arise in the near-zero modal coordinates produced by the inverse spectral basis $W^{-1}$ and which otherwise trigger a $\sim 100\times$ microcode slowdown — are flushed to zero by setting the flush-to-zero and denormals-are-zero bits of the MXCSR control register (ldmxcsr | 0x8040) once at startup; and the per-segment propagator is written as a column-major saxpy over non-aliasing slices so that LLVM emits a clean vectorised inner loop. Thread-level parallelism for the ensemble sweeps (Figure 9, Figure 10) uses std::thread::scope with one worker per operating point over a shared, read-only spectral-table set. The offline table generation — the Schur ordering and Sylvester decoupling of Appendix A, the residual-flexibility MTA correction of Appendix C, and the segment-fusion products — is performed once in Julia ($1.12.6$) and serialised to a binary table that the Rust kernel memory-maps at load; no factorisation happens on the timed path. All reported wall-clock figures are single-run measurements on an otherwise idle machine with the frequency governor left at its default; because the kernel is memory-bandwidth-bound rather than clock-bound (its Roofline position is set by arithmetic intensity, not turbo residency), run-to-run variation is dominated by cache and NUMA placement rather than by core frequency.