SPDX-License-Identifier: AGPL-3.0-or-later

Commercial license available

© Concepts 1996–2026 Miroslav Šotek. All rights reserved.

© Code 2020–2026 Miroslav Šotek. All rights reserved.

ORCID: 0009-0009-3560-0851

Contact: www.anulum.li | protoscience@anulum.li

scpn-quantum-control — Pipeline Performance Benchmarks

Pipeline Performance Benchmarks

Every module in scpn-quantum-control is verified as wired into the pipeline (not decorative) by tests/test_pipeline_wiring_performance.py (155 tests) and per-module pipeline tests embedded in each test file. This page documents the measured wall-time performance for every subsystem.

Test infrastructure

pytest tests/test_pipeline_wiring_performance.py -v -s   # 155 tests, prints benchmarks
pytest tests/test_rust_path_benchmarks.py -v -s           # 68 tests, Rust parity + timing

S2 Scaling Protocol Boundary

The S2 quantum-advantage scaling track starts from an explicit no-claim protocol manifest:

PYTHONDONTWRITEBYTECODE=1 python scripts/export_s2_scaling_protocol.py

The manifest requires classical ODE, dense exact diagonalisation, sparse eigensolver, MPS/tensor-network, and Aer/statevector baselines before any advantage-language figure can be promoted. Hardware rows are optional until credits and a preregistered campaign exist; missing hardware must degrade to a classical scaling/spoofability study rather than a broad quantum-advantage claim.

Rows are validated before promotion with:

PYTHONDONTWRITEBYTECODE=1 python scripts/validate_s2_scaling_rows.py \
  data/s2_advantage_scaling/s2_scaling_rows_schema_fixture_2026-05-06.json

The validator fails if required baselines are absent for any measured size, row keys are missing, baseline labels are unknown, statuses are invalid, successful rows omit non-negative wall-time and memory measurements, or skipped/failed rows omit explanatory notes. It also rejects duplicate (n_qubits, baseline) rows and malformed provenance payloads (metric_payload, command, machine, dependencies, git_commit, and notes). This prevents a partial or ambiguous matrix from being promoted as a complete preregistered comparison.

A lightweight protocol-compliant rehearsal harness is available before the full S2 campaign:

scpn-bench s2-scaling-lite
PYTHONDONTWRITEBYTECODE=1 python scripts/bench_s2_scaling_lite.py --sizes 4,6

It emits required baseline rows for the selected small sizes, measures classical ODE, gated dense exact diagonalisation, gated sparse eigensolver rows, and small-size MPS/TN spoofability diagnostics, and gated Aer/statevector rows where cheap, marks remaining heavier baselines as explicit skipped rows, validates the result against the S2 protocol, and records hardware_submission=false plus advantage_claim=false.

The lite harness records tracemalloc peak bytes for measured rows and exposes explicit size gates:

PYTHONDONTWRITEBYTECODE=1 python scripts/bench_s2_scaling_lite.py \
  --sizes 4,6,8 \
  --max-dense-qubits 6 \
  --max-sparse-qubits 8 \
  --max-tn-qubits 6 \
  --max-statevector-qubits 8

The claim-boundary report is regenerated by the same CLI command and can also be run directly:

PYTHONDONTWRITEBYTECODE=1 python scripts/report_s2_scaling_claim_boundary.py

It records allowed claims, forbidden claims, remaining blockers, validation state, and the protocol claim boundary for the current rows.

docs/campaigns/s2_scaling_readiness_index_2026-05-06.md summarises the current lite baseline rows, allowed claims, forbidden claims, full-campaign blockers, and hardware boundary.

S3 Pulse / Ansatz Design Readiness

The S3 ML-augmented pulse / ansatz track starts from a deterministic no-QPU candidate-ranking gate:

scpn-bench s3-design-ready

The gate ranks structured Kuramoto-XY ansatz candidates by depth, size, and two-qubit gate proxies, and hypergeometric pulse-schedule candidates by analytic infidelity and pulse-count proxies. It records hardware_submission=false and ml_training_performed=false; the output is a schema and claim-boundary layer for later surrogate training, not evidence of a learned optimiser.

Generated artefacts:

• data/s3_pulse_ansatz_design/s3_design_readiness_2026-05-06.json
• data/s3_pulse_ansatz_design/s3_design_readiness_2026-05-06.md

The first surrogate rehearsal is also no-QPU:

scpn-bench s3-design-surrogate

It expands the candidate grid across small system sizes, trains a closed-form ridge linear surrogate on proxy scores, and reports held-out plus per-family metrics. It remains a rehearsal over deterministic proxy labels.

Promoted ansatz candidates can then be checked against exact no-QPU observables:

scpn-bench s3-ansatz-observables

The output records exact statevector energy expectation, dense exact ground energy, energy error, and a synchronisation proxy for the promoted ansatz candidates. This remains an observable validation, not VQE optimisation.

Provider pulse feasibility is checked by metadata-only probes:

scpn-bench s3-pulse-feasibility

The probe compares schedule duration, sample spacing, pulse count, qubit count, and declared pulse or native-XY support. It does not open a provider session or submit hardware jobs.

Human-reviewable hardware-job dossiers are exported with:

scpn-bench s3-hardware-dossiers

These dossiers document the promoted ansatz and pulse follow-up routes, QPU budget state, prerequisites, falsification boundaries, and reproducibility artefacts. They do not authorise execution.

S6 Quantum-Kuramoto Package Boundary

The S6 decoupling track is deliberately gated before any package skeleton is created:

scpn-bench s6-split-audit
scpn-bench s6-boundary-review
scpn-bench s6-api-contract

The split audit inventories reusable, needs-review, and SCPN-specific modules. The boundary review selects a conservative future public surface and keeps package creation blocked until config, provenance, and analysis-dependent rows are separated. The API contract then checks that each proposed quantum_kuramoto.* export has a valid package-local name, imports from the current package, avoids deferred runner/provenance surfaces, and records SCPN-specific rows as warnings rather than silently promoting them.

S7 Logical DLA Parity Roadmap

The S7 fault-tolerant track is gated as a roadmap and resource estimate before any logical-DLA simulation or hardware claim is allowed:

scpn-bench s7-logical-dla-roadmap

The gate regenerates data/s7_logical_dla_parity/logical_dla_parity_roadmap_2026-05-20.json and docs/logical_dla_parity.md. The flat surface-code rows cover 16 oscillators at distances 3, 5, and 7; the MS-QEC cross-check records that the default hierarchical bundle has lower qubit overhead but a non-viable logical rate, so theory review remains mandatory.

S8 Adaptive Branching Readiness

The S8 dynamic-circuit track is gated as a no-submit readiness artefact:

scpn-bench s8-adaptive-branching-readiness

The gate regenerates data/s8_adaptive_branching/adaptive_branching_readiness_2026-05-20.json and docs/adaptive_branching.md. It records the local-order threshold, DLA-parity leakage, and chimera-cluster branch policies plus the backend features needed before live execution. Missing dynamic-circuit support keeps the current readiness state blocked.

S9 Quantum Thermodynamics Readiness

The S9 thermodynamics track is gated as a no-submit calibrated-protocol readiness artefact:

scpn-bench s9-quantum-thermo-readiness

The gate regenerates data/s9_quantum_thermo/quantum_thermo_readiness_2026-05-20.json and docs/quantum_thermo.md. It records a five-point K-sweep for entropy production, heat current, and irreversibility diagnostics while keeping thermodynamic peak claims blocked until theory review, classical reference, and raw-count execution are complete.

S10 Analog-Native Kuramoto Readiness

The S10 analog-native track is gated as a no-submit primitive-accounting and provider-readiness artefact:

scpn-bench s10-analog-native-readiness

The gate regenerates data/s10_analog_native/analog_native_readiness_2026-05-20.json and docs/analog_native_readiness.md. It records native-coupler accounting against a matched declared-tolerance digital Trotter baseline and keeps provider execution plus analog-advantage claims blocked until SDK construction, calibrated unit constraints, and raw execution records are complete.

S11 Quantum Sensing Readiness

The S11 sync-order sensing track is gated as a no-submit QFI/classical-Fisher readiness artefact:

scpn-bench s11-quantum-sensing-readiness

The gate regenerates data/s11_quantum_sensing/quantum_sensing_readiness_2026-05-20.json and docs/quantum_sensing.md. It records a finite K-grid QFI scan, sync-order classical Fisher proxy, and optimal readiness row while keeping hardware submission and sensing-advantage claims blocked until preregistration, shot budget, and raw-count uncertainty evidence are complete.

Phase 3 State/Layout DLA Randomisation

The state/layout mechanism-separation control has a dedicated fail-closed IBM submitter:

python scripts/phase3_state_layout_dla_ibm.py --backend ibm_marrakesh
python scripts/phase3_state_layout_dla_ibm.py --backend ibm_marrakesh --submit --confirm-budget

The first command performs live readiness only. It selects three connected four-qubit windows before outcome data exists, builds the preregistered 495-circuit matrix, live-transpiles every circuit with fixed initial layouts, records depth/gate summaries, and writes a readiness artefact under data/phase3_state_layout_dla/. The second command is the only submitting path and remains blocked unless the depth/gate guards and the 20-minute QPU ceiling pass.

docs/campaigns/s3_design_readiness_index_2026-05-06.md records allowed claims, forbidden claims, and the follow-up path for S3.

Hardware: ML350 Gen8, 2× Xeon E5-2650v2, 128 GB RAM, Ubuntu 24.04. Python: 3.12.3 with Qiskit 1.4.5, Aer 0.17.2. Rust: scpn-quantum-engine 0.2.0 (PyO3 + rayon).

---

1. Bridge Layer

The bridge compiles SCPN coupling matrices (K_nm) into quantum objects.

Knm to Hamiltonian

System size	Compilation time	Output
L=2 (4×4 Hilbert)	<0.1 ms	SparsePauliOp, 2 qubits
L=4 (16×16 Hilbert)	<0.1 ms	SparsePauliOp, 4 qubits
L=8 (256×256 Hilbert)	<0.1 ms	SparsePauliOp, 8 qubits
L=16 (65536×65536 Hilbert)	~6.7 ms	SparsePauliOp, 16 qubits, 256 Pauli terms

Rust path: build_xy_hamiltonian_dense matches Qiskit SparsePauliOp.to_matrix() to machine precision (atol=1e-10). Dense matrix construction in Rust takes 0.02 ms for 3-qubit systems.

Knm to Ansatz

System size	Reps	Parameters	Time
L=2	2	8	<0.1 ms
L=4	2	16	<0.1 ms
L=8	2	32	<0.1 ms
L=16	1	32	<0.1 ms

The Knm-informed ansatz uses coupling topology to determine entanglement gates, producing fewer parameters than generic two_local (12 vs 18 for 3 qubits). Benchmark: knm_informed E=-3.19 beats two_local E=-2.68 at equal iterations.

Knm Construction (Rust)

Function	Input	Time	Speedup
build_knm(16, 0.45, 0.3)	16×16 matrix	0.02 ms	4.7× vs Python

Rust build_knm includes paper27 overrides (L1-L2, L3-L5, L1-L16 boosted couplings) — exact parity with Python build_knm_paper27.

---

2. Phase Solvers

QuantumKuramotoSolver

The core solver maps Kuramoto dynamics to Trotterised XY Hamiltonian evolution. Trajectory sampling uses explicit time boundaries rather than rescaled labels: for non-divisible horizons, the last evolution interval is shortened so the state reaches the reported final time exactly.

Operation	System	Time	Output
run(t_max=0.3, dt=0.1)	4 qubits	16.5 ms	R trajectory: 0.806 → 0.796 → 0.766
evolve(t=0.5, trotter_steps=3)	4 qubits	~5 ms	QuantumCircuit, depth ~45
energy_expectation(sv)	4 qubits	<1 ms	float

Quantum-classical agreement: R(quantum)=0.702 vs R(classical)=0.700 at t=0.2, dt=0.1, trotter_per_step=5 — 0.3% deviation.

Trotter convergence: Error decreases as O(t²/n) for first-order, O(t³/n²) for second-order Suzuki-Trotter. At 4 qubits, second-order produces strictly lower Frobenius error than first-order at equal step count.

PhaseVQE

Operation	System	Time	Output
solve(maxiter=20, seed=42)	2 qubits	~150 ms	E=-3.94, exact=-3.94
solve(maxiter=30, seed=0)	3 qubits	~200 ms	E, exact_energy, gap, params

Variational principle verified: VQE energy >= exact ground energy (within 0.5 tolerance for short optimisation).

VarQITE (Imaginary Time Evolution)

Operation	System	Time	Output
varqite_ground_state(tau=0.5, n_steps=5)	3 qubits	196.8 ms	E=-4.783 vs exact=-4.783

0.0% error — VarQITE achieves exact ITE convergence on 3-qubit system. Energy trajectory: -4.753 → -4.783 (monotonic decrease).

QuantumUPDESolver (Trotter UPDE)

Operation	System	Time	Output
run(n_steps=5, dt=0.05)	4 qubits	20.5 ms	R: 0.806 → 0.804 → 0.796 → 0.783 → 0.765 → 0.743
step(dt=0.1)	3 qubits	~4 ms	R_global, theta

Second-order Trotter (trotter_order=2) passes through correctly to underlying solver. Reset reinitialises state exactly (first step after reset matches first step from fresh solver).

Real-Time Feedback Controller

Operation	System	Backend	Machine	Time	Output
RealtimeSyncFeedbackController.run(5, seed=20260429)	3 qubits, 128 shots/update	Qiskit statevector + NumPy/PyO3 policy fallback	ASRock H510 Pro BTC+, i5-11600K, Ubuntu 24.04.4	68.147 ms	actions=synchronise×5, R_live=0.446→0.372
build_monitored_feedback_circuit(n_rounds=2)	3 system qubits + 1 monitor	Qiskit dynamic circuit	same as above	tested in tests/test_realtime_feedback.py	mid-circuit monitor measurement, conditional reset, conditional correction
scpn-bench s1-feedback	2 qubits, 4 bounded cross-shot scheduler steps, 32 shots/update	FeedbackRunner + RealtimeControllerScheduler + Qiskit statevector controller	current local machine when regenerated	median 21.892 ms, p95 25.570 ms	no-QPU latency, mean final R_live=0.562270, mean final scale=1.094638
scpn-bench s1-feedback	3 qubits, 4 bounded cross-shot scheduler steps, 64 shots/update	same	same	median 18.717 ms, p95 23.680 ms	no-QPU latency, mean final R_live=0.561569, mean final scale=1.095058
scpn-bench s1-feedback	4 qubits, 4 bounded cross-shot scheduler steps, 64 shots/update	same	same	median 25.934 ms, p95 30.455 ms	no-QPU latency, mean final R_live=0.553217, mean final scale=1.100070

Command provenance:

python -c "import time, numpy as np; from scpn_quantum_control.control import RealtimeFeedbackConfig, RealtimeSyncFeedbackController; K=np.array([[0.0,0.35,0.2],[0.35,0.0,0.25],[0.2,0.25,0.0]], dtype=np.float64); omega=np.array([0.1,0.4,0.7], dtype=np.float64); cfg=RealtimeFeedbackConfig(measurement_shots=128, target_r=0.7); start=time.perf_counter(); controller=RealtimeSyncFeedbackController(K, omega, config=cfg); steps=controller.run(5, seed=20260429); elapsed=(time.perf_counter()-start)*1000; print(f'elapsed_ms={elapsed:.3f}'); print('actions=' + ','.join(step.action for step in steps)); print('r_live=' + ','.join(f'{step.r_live:.3f}' for step in steps)); print('next_scale=' + ','.join(f'{step.next_coupling_scale:.3f}' for step in steps))"
scpn-bench s1-feedback
scpn-bench s1-feedback-ready

The s1-feedback benchmark is a control-plane benchmark only. It deliberately excludes provider round-trip latency, IBM queue time, and QPU execution time; an IBM follow-up requires a separate preregistered budget and explicit approval.

Kuramoto Variant Trajectories

Operation	System	Backend	Machine	Time	Output
Run higher-order, monitored, and PT-symmetric variants	4 oscillators, 64 Euler steps per variant, 4 anchored triadic terms	Rust PyO3 higher_order_kuramoto_trajectory, monitored_kuramoto_trajectory, pt_symmetric_kuramoto_trajectory	ASRock H510 Pro BTC+, i5-11600K, Ubuntu 24.04.4	2.205 ms	final R=0.514092,0.775172,0.606729; PT norm final=1.000000

Command provenance:

.venv-linux/bin/python -c "import time, numpy as np; from scpn_quantum_control.phase import HigherOrderKuramotoSpec, MonitoredKuramotoSpec, PTSymmetricKuramotoSpec, build_triadic_ring_terms, simulate_higher_order_kuramoto, simulate_monitored_kuramoto, simulate_pt_symmetric_kuramoto; K=np.array([[0.0,0.45,0.0,0.45],[0.45,0.0,0.45,0.0],[0.0,0.45,0.0,0.45],[0.45,0.0,0.45,0.0]], dtype=np.float64); omega=np.array([0.0,0.6,1.2,2.4], dtype=np.float64); theta0=np.array([0.0,0.8,2.0,4.2], dtype=np.float64); edges, weights=build_triadic_ring_terms(4, 0.25); start=time.perf_counter(); higher=simulate_higher_order_kuramoto(HigherOrderKuramotoSpec(K, omega, edges, weights, theta0=theta0), dt=0.02, n_steps=64); monitored=simulate_monitored_kuramoto(MonitoredKuramotoSpec(K, omega, target_r=0.85, monitor_gain=1.1, measurement_strength=0.25, theta0=theta0), dt=0.02, n_steps=64); pt=simulate_pt_symmetric_kuramoto(PTSymmetricKuramotoSpec(K, omega, np.array([0.08,-0.08,0.04,-0.04], dtype=np.float64), theta0=theta0), dt=0.02, n_steps=64); elapsed=(time.perf_counter()-start)*1000; print(f'elapsed_ms={elapsed:.3f}'); print('backends=' + ','.join([higher.backend, monitored.backend, pt.backend])); print('final_r=' + ','.join(f'{x.final_r:.6f}' for x in [higher, monitored, pt])); print('peak_r=' + ','.join(f'{x.peak_r:.6f}' for x in [higher, monitored, pt])); print('hyperedges=' + str(higher.diagnostics['n_hyperedges'])); print('pt_norm_final=' + f\"{pt.diagnostics['pt_norm'][-1]:.6f}\")"

Automated Witness Discovery

Operation	System	Backend	Machine	Time	Output
Run Bayesian/bandit Kuramoto witness discovery	4 oscillators, 20 evaluated candidates, 48 Euler steps/candidate	Rust PyO3 kuramoto_witness_candidate_features + NumPy RBF UCB/scoring	ASRock H510 Pro BTC+, i5-11600K, Ubuntu 24.04.4	10.145 ms	best score=1.399925; R=0.816060; corr=0.554605; Fiedler=0.909005; source=bayesian_ucb

Command provenance:

.venv-linux/bin/python -c "import time, numpy as np; from scpn_quantum_control.analysis import WitnessDiscoverySpec, discover_kuramoto_witnesses; K=np.array([[0.0,0.5,0.2,0.0],[0.5,0.0,0.4,0.1],[0.2,0.4,0.0,0.3],[0.0,0.1,0.3,0.0]], dtype=np.float64); omega=np.array([0.0,0.35,0.7,1.05], dtype=np.float64); theta0=np.array([0.0,0.7,1.4,2.8], dtype=np.float64); spec=WitnessDiscoverySpec(dt=0.025,n_steps=48,n_initial=8,n_iterations=4,batch_size=3,pool_size=32,seed=20260429,correlation_threshold=0.25,fiedler_threshold=0.2); start=time.perf_counter(); result=discover_kuramoto_witnesses(K, omega, theta0=theta0, spec=spec); elapsed=(time.perf_counter()-start)*1000; print(f'elapsed_ms={elapsed:.3f}'); print('backend=' + result.backend); print('evaluations=' + str(len(result.evaluations))); print('best_score=' + f'{result.best.score:.6f}'); print('best_final_r=' + f'{result.best.final_r:.6f}'); print('best_corr=' + f'{result.best.mean_correlation:.6f}'); print('best_fiedler=' + f'{result.best.fiedler_value:.6f}'); print('best_source=' + result.best.source.value); print('candidate=' + ','.join(f'{v:.6f}' for v in result.best.candidate.as_array()))"

Application Benchmark Plugins

Operation	System	Backend	Machine	Time	Output
Run application plugin benchmark suite	EEG alpha PLV 8-channel, FEP 6-node workflow, ITER MHD 8-mode, IEEE 5-bus grid	scipy.stats.spearmanr + NumPy domain scorers; FEP path through scpn_quantum_control.fep	ASRock H510 Pro BTC+, i5-11600K, Ubuntu 24.04.4	117.635 ms	plugins=eeg_alpha,friston_fep,plasma_iter_mhd,power_grid_ieee5; datasets=eeg_alpha_plv_8ch,friston_fep_6node,iter_mhd_8mode,ieee5bus_power_grid; n=8,6,8,5

The plasma plugin routes the packaged iter_mhd_8mode artifact as a curated reference. The power-grid plugin routes the packaged ieee5bus_power_grid artifact as a curated reference. Direct iter_benchmark() calls require allow_synthetic_reference=True for the built-in synthetic reference, and those results are marked publication_safe=False; direct power_grid_benchmark() calls require allow_builtin_reference=True for the built-in IEEE 5-bus reference and return provenance metadata.

Command provenance:

.venv-linux/bin/python -c "import time; from scpn_quantum_control.applications import run_application_benchmark_suite; start=time.perf_counter(); results=run_application_benchmark_suite(); elapsed=(time.perf_counter()-start)*1000; print(f'elapsed_ms={elapsed:.3f}'); print('plugins=' + ','.join(r.plugin_name for r in results)); print('datasets=' + ','.join(r.dataset_id for r in results)); print('domains=' + ','.join(r.domain for r in results)); print('backends=' + '|'.join(r.backend for r in results)); print('n_oscillators=' + ','.join(str(r.n_oscillators) for r in results)); print('key_metrics=' + '|'.join(','.join(f'{k}={v:.6f}' for k, v in sorted(r.metrics.items())[:3]) for r in results))"

Analog Kuramoto Backend Compiler

Operation	System	Backend	Machine	Time	Output
Compile one problem for neutral_atoms, circuit_qed, and continuous_variable	4 oscillators, 5 non-zero couplers	NumPy/PyO3 analog term kernel + serialisable native schemas	ASRock H510 Pro BTC+, i5-11600K, Ubuntu 24.04.4	0.610 ms	schemas=native_ahs_v1, exchange_resonator_v1, cv_gaussian_schedule_v1
Compile one hybrid digital-analog execution plan	4 oscillators, 5 non-zero couplers split into 2 analog + 3 digital residual couplers	Rust PyO3 partition kernel + Qiskit Trotter residual compiler	ASRock H510 Pro BTC+, i5-11600K, Ubuntu 24.04.4	6.275 ms	schema=hybrid_digital_analog_v1, Rust partition=True, digital depth=1
Certify and witness one DLA-protected logical synchronisation memory	4 logical oscillators, distance-3 repetition blocks, 12 physical qubits	Rust PyO3 protected-memory mask + probability metrics	ASRock H510 Pro BTC+, i5-11600K, Ubuntu 24.04.4	0.390 ms	DLA dimension=8,388,606; protected=0.950; sync=0.950
Simulate one DLA-protected scar memory revival	4 logical oscillators, distance-3 repetition blocks, 12 physical qubits	NumPy diagonal phase evolution + Rust PyO3 trajectory protected-sector metrics	ASRock H510 Pro BTC+, i5-11600K, Ubuntu 24.04.4	9.478 ms	scar states=2; revival=1.000000; midcycle=0.000000; protected=1.000000; parity leakage=0.000000
Run real-data synchronisation forecasting suite	IBM Heron r2 4-oscillator hardware trace + IEEE 5-bus topology replay	Stored hardware exact_R baseline + Rust PyO3 kuramoto_trajectory for topology replay + NumPy train-window calibration	ASRock H510 Pro BTC+, i5-11600K, Ubuntu 24.04.4	17.142 ms	hardware held-out MSE 0.104217 -> 0.000104; IEEE replay backend=rust:kuramoto_trajectory; suite passes

Command provenance:

python -c "import time, numpy as np; from scpn_quantum_control.hardware.analog_kuramoto import compile_analog_kuramoto; K=np.array([[0.0,0.5,-0.25,0.125],[0.5,0.0,0.2,0.0],[-0.25,0.2,0.0,0.3],[0.125,0.0,0.3,0.0]], dtype=np.float64); omega=np.array([0.1,-0.2,0.3,0.0], dtype=np.float64); start=time.perf_counter(); programs=[compile_analog_kuramoto(K, omega, platform=p, duration=1.25) for p in ('neutral_atoms','circuit_qed','continuous_variable')]; elapsed=(time.perf_counter()-start)*1000; print(f'elapsed_ms={elapsed:.3f}'); print('platforms=' + ','.join(p.platform.value for p in programs)); print('couplers=' + ','.join(str(p.n_couplers) for p in programs)); print('schemas=' + ','.join(p.payload['schema'] for p in programs))"
.venv-linux/bin/python -c "import time, numpy as np; import scpn_quantum_engine as e; from scpn_quantum_control.hardware.hybrid_digital_analog import compile_hybrid_digital_analog; K=np.array([[0.0,0.8,-0.4,0.1],[0.8,0.0,0.3,0.0],[-0.4,0.3,0.0,-0.2],[0.1,0.0,-0.2,0.0]], dtype=np.float64); omega=np.array([0.1,-0.2,0.05,0.3], dtype=np.float64); start=time.perf_counter(); program=compile_hybrid_digital_analog(K, omega, platform='circuit_qed', duration=1.25, max_analog_couplers=2, trotter_steps=3); elapsed=(time.perf_counter()-start)*1000; print(f'elapsed_ms={elapsed:.3f}'); print('rust_partition=' + str(hasattr(e, 'hybrid_coupling_partition'))); print('analog=' + str(program.n_analog_couplers)); print('digital=' + str(program.n_digital_couplers)); print('schema=' + program.payload['schema']); print('digital_depth=' + str(program.digital_circuit.depth()))"
.venv-linux/bin/python -c "import time, numpy as np, scpn_quantum_engine as e; from scpn_quantum_control.qec import DLAProtectedSubspaceSpec, certify_dla_protected_subspace, evaluate_dla_protected_memory; spec=DLAProtectedSubspaceSpec(n_logical=4, code_distance=3, target_parity=0); probs=np.zeros(spec.hilbert_dim, dtype=np.float64); probs[0]=0.55; probs[(1 << spec.n_physical)-1]=0.40; probs[0b000000000111]=0.03; probs[0b000000010111]=0.02; start=time.perf_counter(); cert=certify_dla_protected_subspace(spec); result=evaluate_dla_protected_memory(probs, spec=spec); elapsed=(time.perf_counter()-start)*1000; print(f'elapsed_ms={elapsed:.3f}'); print('rust_mask=' + str(hasattr(e, 'dla_protected_memory_mask'))); print('rust_metrics=' + str(hasattr(e, 'dla_protected_memory_metrics'))); print('dla_dim=' + str(cert.physical_dla_dimension)); print('protected=' + f'{result.protected_weight:.3f}'); print('sync=' + f'{result.sync_weight:.3f}'); print('passes=' + str(result.passes))"
.venv-linux/bin/python -c "import time; from scpn_quantum_control.qec import DLAProtectedScarSpec, build_dla_protected_scar_prototype, simulate_dla_protected_scar_memory; start=time.perf_counter(); spec=DLAProtectedScarSpec(); prototype=build_dla_protected_scar_prototype(spec); result=simulate_dla_protected_scar_memory(prototype); elapsed=(time.perf_counter()-start)*1000; print(f'elapsed_ms={elapsed:.3f}'); print('backend=' + result.backend); print('n_physical=' + str(spec.memory_spec.n_physical)); print('scar_states=' + str(prototype.n_scar_states)); print('revival=' + f'{result.final_revival_fidelity:.6f}'); print('midcycle=' + f'{result.midcycle_survival:.6f}'); print('protected=' + f'{result.min_protected_weight:.6f}'); print('parity_leakage=' + f'{result.max_parity_leakage:.6f}'); print('passes=' + str(result.passes))"
.venv-linux/bin/python -c "import time; from scpn_quantum_control.forecasting import run_real_data_sync_forecast_suite; start=time.perf_counter(); results=run_real_data_sync_forecast_suite(include_topology_replay=True); elapsed=(time.perf_counter()-start)*1000; print(f'elapsed_ms={elapsed:.3f}'); print('passes=' + str(all(r.passes for r in results))); [print(f'{r.dataset.name}: baseline_mse={r.baseline.holdout_mse:.9f}; calibrated_mse={r.calibrated.holdout_mse:.9f}; improvement_pct={100*r.improvement_fraction:.2f}; backend={r.baseline.backend}; source={r.dataset.source_kind}') for r in results]"

Adiabatic State Preparation

Operation	System	Time	Output
adiabatic_ramp(K_target=3.0, T=5.0, n_steps=15)	3 qubits	54.4 ms	min_gap=0.0012 at K=2.80

This dense finite-size scan records the observed gap minimum and fidelity loss for the selected schedule. The reported K value is a finite-size gap minimum for this scan, not a standalone transition-location proof.

Floquet-Kuramoto (Time Crystal)

Operation	System	Time	Output
floquet_evolve(K=1.0, amp=0.5, freq=2.0)	2 qubits	0.6 ms	R(t), subharmonic ratio
scan_drive_amplitude(5 amplitudes)	2 qubits	~3 ms	subharmonic ratio per amplitude

DTC candidate detection via subharmonic_ratio > threshold. Drive signal oscillates between K_base(1-amp) and K_base(1+amp) as expected.

---

3. Hardware Layer

HardwareRunner (Simulator)

Operation	System	Time	Output
connect()	AerSimulator	~50 ms	Backend ready
transpile(GHZ circuit)	4 qubits	~10 ms	ISA circuit, depth ~10
run_sampler(shots=1000)	4 qubits	~100 ms	counts dict
circuit_stats()	—	<1 ms	depth, n_qubits, ECR count

Fractional gates: With use_fractional_gates=True, Kuramoto circuit depth reduces from ~80 to ~60 (25% reduction) for 4 qubits, 2 Trotter steps. RZZ gates remain native instead of decomposing to ECR+RZ.

Noise Model (Heron r2)

Operation	System	Time	Output
heron_r2_noise_model()	—	<1 ms	NoiseModel
Noisy Bell pair (10k shots)	2 qubits	~150 ms	non-ideal counts
Noisy Kuramoto R comparison	3 qubits	1349 ms	R_clean=0.734, R_noisy=0.734

At default Heron r2 parameters (CZ error=0.005), noise degradation is minimal (R_clean ≈ R_noisy). Higher CZ error (0.1) produces measurable non-ideal counts.

Trapped-Ion Backend

Operation	System	Time	Output
transpile_for_trapped_ion(allow_proxy_basis=True)	4 qubits	~5 ms	All-to-all connectivity, no SWAPs

Kuramoto circuits transpile without SWAP gates under the all-to-all representative model. The output circuit is metadata-labelled as a CX-basis proxy for native MS/RXX-style trapped-ion operations; unitarity preservation is verified by Operator equivalence.

Circuit Depth Regression

System	Trotter reps	Transpiled depth	Gate count
2q, 1 rep	1	<50	<100
4q, 1 rep	1	<100	<300
4q, 3 reps	3	134	~200
8q, 1 rep	1	<300	<600
16q, 1 rep	1	<1000	~1500

Depth scales sub-linearly with reps (3 reps < 4× depth of 1 rep due to gate cancellation in transpilation).

QASM Export

Operation	System	Time	Output
export_trotter_qasm(K, omega, t=0.5, reps=3)	4 qubits	3.4 ms	1903 chars, 48 gates

Exports OpenQASM 3.0 with qubit declarations and gate definitions.

---

4. Error Mitigation

Zero-Noise Extrapolation (ZNE)

Operation	System	Time	Output
Fold at scales [1,3,5] + extrapolate	3 qubits	34.4 ms	R_ZNE estimate

Folded circuits preserve unitarity (norm=1.0 at all odd scales). Fit residual >= 0. On noiseless simulator, all scale values are identical (folding is identity).

Probabilistic Error Cancellation (PEC)

Operation	System	Time	Output
pauli_twirl_decompose(0.05)	1 qubit	<0.01 ms	4 coefficients
pec_sample(circuit, p=0.05, n=200)	1 qubit	160.9 ms	mitigated =-1.07 (ideal -1.0)

Rust path: pec_coefficients(p) matches Python pauli_twirl_decompose(p) to machine precision (atol=1e-10). Rust pec_sample_parallel(100k samples) takes 49-91 ms using rayon parallelism.

Quasi-probability invariant: identity coefficient > 1, error coefficients < 0, sum = 1.0 (trace preservation).

---

5. Quantum Error Correction

ControlQEC (Surface Code)

Operation	System	Time	Output
ControlQEC(distance=3)	18 data qubits	<0.1 ms	Decoder ready
get_syndrome() + decode()	d=3	0.6 ms	correction vector

Below-threshold correction: >80% success at p=0.01. Above-threshold: significant failure at p=0.3. Zero-error syndrome is all-zero (verified).

FaultTolerantUPDE (Repetition Code)

Operation	System	Time	Output
build_step_circuit(dt=0.1)	2 osc, d=3	<0.1 ms	10-qubit circuit
step_with_qec(dt=0.1)	3 osc, d=3	0.3 ms	syndromes, errors_detected

Qubit layout: n_osc × (2d-1) physical qubits. Contains RZZ (transversal coupling), CX (encoding + syndrome), RZ (field terms).

SurfaceCodeUPDE

Operation	System	Time	Output
SurfaceCodeUPDE(n_osc=4, code_distance=3)	4 oscillators	<1 ms	Resource model

Total physical qubits = n_osc × (2d²-1). For d=3: 4 × 17 = 68 physical qubits.

---

6. QSNN (Quantum Spiking Neural Network)

QuantumSynapse

Operation	Time	Output
apply(circuit, pre, post)	<0.01 ms	CRy gate appended

theta = pi × (w - w_min) / (w_max - w_min). Effective weight = sin²(theta/2). Pre=|1> → post rotates; pre=|0> → post stays |0> (controlled rotation).

QuantumLIFNeuron

Operation	Time	Output
step(input_current=1.5)	~1 ms	spike ∈

Membrane equation: v(t+1) = v(t) - (dt/tau)(v(t) - v_rest) + RIdt. Quantum mapping: P(spike) = sin²(theta/2) where theta encodes membrane potential.

QuantumSTDP

Operation	Time	Output
update(syn, pre=1, post=1)	<0.01 ms	weight updated

Hebbian LTP: pre+post fire → weight increases. LTD: pre fires, post doesn't → weight decreases. No pre spike → no change (verified).

QSNNTrainer

Operation	System	Time	Output
train(X, y, epochs=3)	2×2 layer	47.6 ms	loss history

Parameter-shift gradient: g = (L(+pi/2) - L(-pi/2)) / 2. Gradient sign flips for opposite targets (antisymmetry). Zero learning rate → zero weight change (verified to 1e-14). Forward probabilities bounded [0,1].

SNNQuantumBridge

Operation	System	Time	Output
forward(spike_history)	4→3 neurons	2.2 ms	output currents

Spike-to-rotation: firing_rate × pi ∈ [0, pi]. Higher rate → larger angle (monotonic). Measurement-to-current: P(|1>) × scale.

---

7. Identity Layer (Arcane Sapience)

IdentityAttractor

Operation	System	Time	Output
solve(maxiter=30, seed=42)	3 qubits	108.4 ms	E_0=-4.749, gap=1.383

Robustness gap = E_1 - E_0. Gap=1.383 provides strong identity protection. Eigenvalues sorted ascending. Variational bound: E_vqe >= E_exact. Stronger coupling → larger gap (verified).

Identity Fingerprint

Operation	System	Time	Output
identity_fingerprint(K, omega)	4 qubits	~150 ms	commitment (SHA-256 hex)

Returns dict with commitment, spectral data (fiedler, eigenvalues), ground_energy, n_parameters. Different K → different commitment. Spectral data deterministic.

Challenge-Response Protocol

Operation	System	Time	Output
prove_identity(K, challenge)	3 qubits	<1 ms	response bytes
verify_identity(K, challenge, response)	3 qubits	<1 ms	True/False

Wrong K produces wrong response → verification fails. Different challenges → different responses (no replay).

Robustness Certificate

Operation	System	Time	Output
compute_robustness_certificate(K, omega)	3 qubits	0.9 ms	gap=1.383, P_transition=5.2e-5

P_transition = 5.2×10⁻⁵ — probability of identity confusion under noise. Fidelity at depth: deeper circuits → lower fidelity (decoherence monotonicity).

---

8. Cryptographic Layer

Key Hierarchy

Operation	System	Time	Output
key_hierarchy(K, phases, R, nonce)	4 layers	0.11 ms	master (32 bytes) + 4 layer keys

All layer keys unique. Master key differs from all layer keys. Same inputs → same keys (deterministic). Different R or nonce → different keys. verify_key_chain() detects tampered master, tampered layer keys, wrong nonce.

Topology Commitment

Operation	System	Time	Output
topology_commitment(K)	4×4 matrix	<0.1 ms	32-byte SHA-256

Deterministic hash of coupling topology. Combined pipeline (hierarchy + fingerprint + commitment): 0.46 ms.

SCPN-QKD Protocol

Operation	System	Time	Output
scpn_qkd_protocol(K, omega, alice, bob)	4 qubits	692 ms	QBER, raw keys, Bell

QBER ∈ [0, 1]. Ground energy < 0. Raw key shapes match qubit allocation. Secure key length >= 0.

Evolving Key Phases

Operation	System	Time	Output
evolve_key_phases(K, omega, theta_0, t=0.5)	4 layers	~1 ms	(n_layers, n_samples) trajectory

Kuramoto ODE integration via solve_ivp(RK45). Initial condition preserved at t=0. All values finite. ODE failure → RuntimeError with message.

---

9. Analysis Layer

Finite-Size Scaling

Operation	System	Time	Output
finite_size_scaling(sizes=[2,3,4])	3 sizes	0.8 ms	K_c per size + extrapolation

K_c values finite. gap_min > 0. Extrapolation via BKT or power-law fit.

H1 Persistence

Operation	System	Time	Output
scan_h1_persistence(omega, n_points=10)	4 osc	14.9 ms	K_critical, p_h1

K_critical > 0. p_h1 ∈ [0, 1]. Vortex densities bounded. K values sorted.

OTOC Synchronisation Probe

Operation	System	Time	Output
otoc_sync_scan(K, omega, n_K=6, n_t=8)	3 qubits	7.6 ms	Lyapunov, R_classical

R_classical bounded [0, 1]. Lyapunov values finite. OTOC detects transition: True.

Berry Phase

Operation	System	Time	Output
berry_phase_scan(omega, T, k_range)	3 qubits	6.6 ms	curvature peak at K=0.75

Fidelity ∈ [0, 1]. Spectral gap > 0. Curvature finite. Fidelity susceptibility chi_F peaks near BKT transition (max chi_F = 0.005).

Loschmidt Echo / DQPT

Operation	System	Time	Output
loschmidt_quench(K_i=0.5, K_f=3.0)	3 qubits	0.8 ms	3 cusps detected

|G(0)| = 1 exactly. Rate function r(0) = 0. Times monotonic. No-quench: |G(t)| = 1 for all t. Large quench: amplitude oscillations.

Krylov Complexity

Operation	System	Time	Output
krylov_complexity(H, Z0, t_max=5.0)	3 qubits	155 ms	peak K(t) = 3.031

K(0) = 0 (operator starts in first basis element). K(t) >= 0. K(t) <= d² (bounded by Hilbert space dimension). Lanczos b_n decay for finite dimension.

Rust path: lanczos_b_coefficients produces same coefficients as Python (verified to atol=1e-6 on first few b_n).

Entanglement Entropy

Operation	System	Time	Output
entanglement_at_coupling(omega, T, K=2.0)	4 qubits	0.3 ms	S=0.928, gap=0.224

S ∈ [0, log₂(d)] where d = 2^(n/2). Schmidt gap ∈ [0, 1]. Weak coupling → S ≈ 0 (product state). Strong coupling → S > 0. Schmidt gap closes near BKT.

QFI Criticality

Operation	System	Time	Output
qfi_vs_coupling(omega, T, k_range)	3 qubits	8.5 ms	peak QFI=0.225 at K=3.07

QFI >= 0. Total QFI >= max single-generator QFI. Peak at K=3.07 confirms criticality-enhanced quantum correlations.

Quantum Speed Limit

Operation	System	Time	Output
compute_qsl(K, omega, t=1.0)	3 qubits	10.4 ms	tau_MT, tau_ML bounds

Mandelstam-Tamm bound tau_MT >= 0. Margolus-Levitin bound tau_ML >= 0. Actual time tau_actual >= both bounds (QSL is a lower bound).

Spectral Form Factor

Operation	System	Time	Output
compute_sff(K, omega, n_times=20)	4 qubits	1.2 ms	r_bar=0.488, gap=1.132

K(t=0) = 1 exactly (trace identity). SFF ∈ [0, 1]. Times monotonic. Level spacing ratio r_bar = 0.488 (near GOE Wigner-Dyson 0.536 — quantum chaotic).

Magic (Non-stabilizerness)

Operation	System	Time	Output
magic_vs_coupling(omega, T, k_range)	3 qubits	~5 ms	SRE peak

SRE (stabiliser Renyi entropy) M₂ >= 0. Weak coupling → M₂ ≈ 0 (stabiliser ground state). Strong coupling → M₂ > 0 (magic resource). Berry curvature F_μν is antisymmetric (traceless).

Lindblad NESS

Operation	System	Time	Output
compute_ness(omega, T, K=2.0, gamma=0.1)	2 qubits	~1 ms	R_ness, purity

Purity ∈ [1/d, 1]. R_ness ∈ [0, 1]. gamma=0 → NESS = ground state (R_ness ≈ R_ideal). Purity decreases with noise (generally).

Hamiltonian Learning

Operation	System	Time	Output
measure_correlators + learn_hamiltonian	3 qubits	34.6 ms	loss=0, corr_error=0

Correlator matrix symmetric, zero diagonal, bounded [-2, 2]. Learned K symmetric, non-negative. Perfect recovery for 3-qubit system (loss=0). Self-consistent: true K as init → near-zero error.

Hamiltonian Self-Consistency

Operation	System	Time	Output
self_consistency_from_exact(K, omega)	2 qubits	10.9 ms	Frobenius=1.81, loss=0

2-qubit inverse problem is degenerate: loss=0 but Frobenius error=1.81 because different K values produce identical correlators.

XXZ Phase Diagram

Operation	System	Time	Output
anisotropy_phase_diagram(3δ × 6K)	3 qubits	36.1 ms	K_c(Δ=0)=0.5, K_c(Δ=0.5)=1.2

XY (Δ=0) and Heisenberg (Δ=1) produce different gap structure. All gaps > 0.

QRC Phase Detector

Operation	System	Time	Output
qrc_phase_detection(8 train, 2 test)	3 qubits	39.3 ms	accuracy=100%, 36 features

Self-probing: reservoir features from ground state observables. Linear readout achieves perfect phase classification on well-separated data.

---

10. Application Layer

Quantum Reservoir Computing

Operation	System	Time	Output
reservoir_ridge_regression(12 samples)	3 qubits	33.9 ms	MSE=0.022

Feature matrix has non-trivial rank (expressive reservoir). Higher weight → more features. Ridge regression produces actionable predictions.

Quantum Kernel

Operation	System	Time	Output
compute_kernel_matrix(5 samples)	3 qubits	16.1 ms	PSD Gram matrix

Mercer conditions verified: symmetric, PSD (min eigenvalue=0.028 > 0), diagonal=1. K(x,x) = 1. Close inputs → high overlap (>0.95). Different Knm → different kernel.

Disruption Classifier

Operation	System	Time	Output
run_disruption_benchmark(10+5, allow_synthetic=True)	3 qubits	297 ms	accuracy=80%, source_mode=synthetic

Kernel Gram matrix symmetric + PSD. Binary predictions. Accuracy bounded [0, 1]. This benchmark is generated-data only and is not publication-safe without measured plasma diagnostics.

Quantum Disruption (ITER)

Operation	System	Time	Output
predict(features)	5 qubits	4.6 ms	risk=0.495
DisruptionBenchmark(20+10, allow_synthetic=True, 2 epochs)	5 qubits	11.9 s	accuracy=70%, source_mode=synthetic

Feature normalisation clamps to [0, 1]. Prediction deterministic for same params. Circuit depth > 0. Training updates parameters. The benchmark dataset is generated and is not publication-safe without measured plasma diagnostics or from_fusion_core_shot() inputs.

FMO Photosynthetic Benchmark

Operation	System	Time	Output
fmo_benchmark(K, omega, allow_builtin_reference=True)	7 sites	1.4 ms	topology ρ=0.304, source_mode=builtin_literature_reference

SCPN vs FMO topology correlation ρ=0.304 (weak positive). FMO self-comparison: ρ=1.0. FMO coupling: symmetric, non-negative, zero diagonal, 7×7. The packaged Adolphs-Renger reference requires explicit opt-in and is marked publication_safe=False; publication-safe FMO claims must pass measured fmo_coupling and fmo_frequencies arrays.

Quantum Advantage Scaling — Exact-Simulation Crossover Scope

These rows describe only the exact Hilbert-space crossover benchmark: classical exact diagonalisation and Python/SciPy ODE vs statevector Qiskit wall time for the same order-parameter trajectory observable. They are not evidence for broad observable-level quantum advantage.

Operation	System	Time	Output
run_scaling_benchmark(sizes=[3,4])	3-4 qubits	101 ms	timing comparison

n=3: classical=23 ms, quantum=11 ms (quantum wins). n=4: classical=26 ms, quantum=34 ms (classical wins). Crossover near n=4.

Quantum Advantage Scaling — Provenance Run

Measured on 2026-04-29 on local host aaarthuus (Intel i5-11600K, 6 cores / 12 threads, Linux 6.17.0-22-generic). Backend versions: Python 3.12.3, NumPy 2.4.4, SciPy 1.17.1, Qiskit 2.4.0.

Command:

python -c "from scpn_quantum_control.benchmarks.quantum_advantage import run_scaling_benchmark; results = run_scaling_benchmark(sizes=[3,4,5], t_max=0.2, dt=0.1); print('n,t_classical_ms,t_quantum_ms,crossover'); [print(f'{r.n_qubits},{r.t_classical_ms:.3f},{r.t_quantum_ms:.3f},{r.crossover_predicted}') for r in results]"

Backend	Machine	Command	Dependency	Commit	n	Classical exact evolution + diagonalisation	Qiskit statevector Trotter	Predicted crossover
python-scipy-qiskit-local	Intel i5-11600K / Linux 6.17.0-22-generic	python -c "from scpn_quantum_control.benchmarks.quantum_advantage import run_scaling_benchmark; results = run_scaling_benchmark(sizes=[3,4,5], t_max=0.2, dt=0.1); print('n,t_classical_ms,t_quantum_ms,crossover'); [print(f'{r.n_qubits},{r.t_classical_ms:.3f},{r.t_quantum_ms:.3f},{r.crossover_predicted}') for r in results]"	Python 3.12.3; NumPy 2.4.4; SciPy 1.17.1; Qiskit 2.4.0	638bba955e84044f887b537bc0dce46234cb80fc	3	53.197 ms	68.128 ms	6
python-scipy-qiskit-local	Intel i5-11600K / Linux 6.17.0-22-generic	python -c "from scpn_quantum_control.benchmarks.quantum_advantage import run_scaling_benchmark; results = run_scaling_benchmark(sizes=[3,4,5], t_max=0.2, dt=0.1); print('n,t_classical_ms,t_quantum_ms,crossover'); [print(f'{r.n_qubits},{r.t_classical_ms:.3f},{r.t_quantum_ms:.3f},{r.crossover_predicted}') for r in results]"	Python 3.12.3; NumPy 2.4.4; SciPy 1.17.1; Qiskit 2.4.0	638bba955e84044f887b537bc0dce46234cb80fc	4	23.316 ms	53.197 ms	6
python-scipy-qiskit-local	Intel i5-11600K / Linux 6.17.0-22-generic	python -c "from scpn_quantum_control.benchmarks.quantum_advantage import run_scaling_benchmark; results = run_scaling_benchmark(sizes=[3,4,5], t_max=0.2, dt=0.1); print('n,t_classical_ms,t_quantum_ms,crossover'); [print(f'{r.n_qubits},{r.t_classical_ms:.3f},{r.t_quantum_ms:.3f},{r.crossover_predicted}') for r in results]"	Python 3.12.3; NumPy 2.4.4; SciPy 1.17.1; Qiskit 2.4.0	638bba955e84044f887b537bc0dce46234cb80fc	5	65.367 ms	66.430 ms	6

These rows are a smoke-scale provenance table, not the publication scaling claim. Sparse SciPy emitted conversion warnings during the run, which is useful context for interpreting the small-n noise.

Quantum Advantage — Classical / Rust / GPU Matrix (Gap B item 3)

Gap B item 3 now tracks the pre-QPU guardrail matrix in results/classical_rust_gpu_matrix_2026-05-03.json. It records, for matching problem sizes:

• classical order-parameter ODE (SciPy path),
• classical exact diagonalisation (CPU dense and sparse),
• exact diagonalisation on GPU when available,
• Rust ODE kernel timing.

Run:

python scripts/bench_quantum_advantage_classical_matrix.py \
  --sizes 4,6,8,10 \
  --max-dense-dim 64 \
  --max-sparse-dim 64 \
  --max-sparse-iter 20 \
  --output results/classical_rust_gpu_matrix_2026-05-03.json

Matrix rows are used as the classical-rank evidence gate before expanding frontier QPU campaigns.

Broad Quantum Advantage (Observable-Level) — Not yet claimed

No broad observable-level advantage claim is documented yet. The committed evidence must distinguish:

• exact-simulation crossover (this section),
• observable-level hardware-native/algorithm-native comparisons,
• matched-accuracy and matched-budget classical alternatives (including sparse and tensor-network paths), and
• gating path that justifies any future broad-advantage wording in captions, README content, or preprint claims.

---

11. Bridge Adapters

SSGF Adapter

Operation	System	Time	Output
W→H→encode→decode	4 oscillators	1.5 ms	R_global=0.767
SSGFQuantumLoop.quantum_step	4 oscillators	~9 ms	theta updated, R returned

Encoding: 2 gates per oscillator (Ry + Rz). Normalisation preserved. Uniform phases → R ≈ 1. Opposite phases → R ≈ 0.

SSGF Spectral Bridge

Operation	System	Time	Output
spectral_bridge_analysis(K, omega)	4 oscillators	0.2 ms	fiedler=0.872, QPE=7 bits

Fiedler > 0 for connected graph. Eigenvalues non-negative (Laplacian PSD). Disconnected graph → fiedler=0. QPE bits estimate for spectral resolution.

SSGF W Adapter

Operation	System	Time	Output
adapt_w_from_quantum(K, theta, lr=0.1)	4 oscillators	4.9 ms	max_update=0.027

W_updated symmetric, non-negative, zero diagonal. Correlators symmetric. lr=0 → no change. W changes with non-zero lr.

Orchestrator Adapter

Operation	System	Time	Output
from_orchestrator_state → to_scpn_control_telemetry	3 layers	0.07 ms	regime, R, stability

Handles both dataclass and dict payloads. Legacy field names (locks, cross_alignment, stability, regime) resolved automatically.

Orchestrator Feedback

Operation	System	Time	Output
compute_orchestrator_feedback(K, omega)	4 qubits	~0.5 ms	action, confidence, R_global

Actions: advance, hold, rollback. Confidence ∈ [0, 1]. R_global ∈ [0, 1]. Custom thresholds supported.

---

12. PGBO (Parameter-space Geometry Bridge)

Operation	System	Time	Output
compute_pgbo_tensor(K, omega)	4 qubits	6.7 ms	metric (6×6), curvature (6×6)

Quantum Fisher metric: symmetric, PSD (det >= 0). Berry curvature: antisymmetric (traceless). Parameter count: C(n,2) upper-triangle couplings.

---

13. TCBO Observer

Operation	System	Time	Output
compute_tcbo_observables(K, omega)	4 qubits	4.6 ms	p_h1, TEE, string_order

p_h1 ∈ [0, 1]. TEE finite. |string_order| <= 1. beta_0 + beta_1 ≈ 1 (connected components + loops = 1). Different coupling → different observables.

---

14. Trotter Error Analysis

Commutator Bounds

Operation	System	Time	Output
commutator_norm_bound + optimal_dt	4 qubits	<0.1 ms	gamma=5.344, dt*=0.004, n_steps=268

Equal frequencies → gamma=0 (no Trotter error). Heterogeneous frequencies → larger gamma. Second-order bounds use exact nested-commutator spectral norms for small systems and a Pauli coefficient-norm upper bound for larger systems. Optimal dt respects the epsilon target.

Trotter Error Sweep

Operation	System	Time	Output
trotter_error_sweep(3t × 3reps)	3 qubits	~3.9 ms	2D error map

Functional, non-isolated timing (median of 5, host load ≈ 5.3). The two-group product formula is measured with direct dense matrix exponentials, replacing the earlier circuit-construction path.

Error at t=0: < 1e-10. Error decreases with reps. Error increases with time. First-order error shows quadratic time scaling; the symmetric second order (order=2) shows cubic time scaling and lower error at fixed (t, reps).

---

15. Experiment Registry

Operation	Time	Output
List all experiments	0.18 ms	20 registered experiments

Every experiment has: runner as first param, docstring > 10 chars, lowercase underscore name, and no hidden experiment route. At least half accept shots parameter.

---

16. Cutting Runner (Large-Scale)

Operation	System	Time	Output
run_cutting_simulation(n=16, max=8, allow_partition_energy_estimate=True)	16 oscillators	39.3 ms	2 partitions, R=1.0, energy scope partition_local_sum
run_cutting_simulation(n=24, max=8, allow_partition_energy_estimate=True)	24 oscillators	~53 ms	3 partitions, energy scope partition_local_sum
run_cutting_simulation(n=32, max=8, allow_partition_energy_estimate=True)	32 oscillators	~60 ms	4 partitions, energy scope partition_local_sum

Partitions: ceil(n/max_partition_size). R per partition bounded [0, 1]. Combined R bounded [0, 1]. For multi-partition runs, energy is a labelled partition-local diagnostic and the omitted cross-partition coupling L1 norm is reported in the result; full-system energy is not claimed.

---

17. GUESS Symmetry-Decay ZNE (April 2026)

mitigation/symmetry_decay.py (Rust path: fit_symmetry_decay, guess_extrapolate_batch).

Operation	Input	Time	Notes
learn_symmetry_decay	5 noise scales (Rust)	< 1 µs	least-squares α fit
learn_symmetry_decay	5 noise scales (Python fallback)	< 1 µs	numpy polyfit
guess_extrapolate	single observable	< 0.1 µs	analytic correction
guess_extrapolate_batch	1,000 observables (Rust)	< 50 µs	rayon-parallel
guess_extrapolate_batch	100,000 observables (Rust)	< 2 ms	scales linearly

Pipeline test: tests/test_pipeline_wiring_performance.py::TestGUESSPipeline.

18. DynQ Topology-Agnostic Qubit Mapper (April 2026)

hardware/qubit_mapper.py (Rust path: score_regions_batch).

Operation	Input	Time
build_calibration_graph	156-qubit heavy-hex	< 5 ms
detect_execution_regions	156 qubits	< 50 ms
dynq_initial_layout (full pipeline)	156 qubits → 5-qubit layout	< 100 ms
score_regions_batch (Rust)	100 regions × 50 qubits	< 5 ms

Pipeline test: tests/test_pipeline_wiring_performance.py::TestDynQPipeline.

19. Pulse Shaping — ICI + (α,β)-Hypergeometric (April 2026)

phase/pulse_shaping.py (Rust paths: hypergeometric_envelope_batch, ici_mixing_angle_batch, ici_three_level_evolution_batch).

Operation	Input	Python	Rust	Speedup
hypergeometric_envelope	200 points	0.4 ms	0.1 ms	4×
hypergeometric_envelope	10,000 points	114.5 ms	2.6 ms	44×
ici_mixing_angle	1,000 points	0.05 ms	0.04 ms	parity-checked
ici_three_level_evolution	2,000 points	68.30 ms	0.04 ms	1,665×
build_trotter_pulse_schedule	4 qubits, K all-to-all	1.2 ms	n/a	full schedule

Verified parity (Rust vs Python): max abs diff \(4.97 \times 10^{-14}\) for \(n_\text{points} = 500\) on ici_three_level_evolution.

Pipeline test: tests/test_pipeline_wiring_performance.py::TestPulseShapingPipeline.

20. Phase 1 IBM Hardware Campaign (April 2026)

Real-world QPU runs on ibm_kingston (Heron r2, 156 q):

Sub-phase	Circuits	Wall (queue + exec)	QPU rate
Pipe cleaner	2	~0.1 s	1.0 s/circuit
Phase 1 (A/B/C)	42	44.1 s	0.65 s/circuit
Phase 1.5 (D/E)	72	56.7 s	0.55 s/circuit
Phase 2 exhaust (F/G/H/I)	138	97.5 s	0.55 s/circuit
Phase 2.5 final burn (J)	90	65.1 s	0.55 s/circuit
Total	344	~264 s	~0.55 s/circuit

Reproducible from data/phase1_dla_parity/*.json via scripts/analyse_phase1_dla_parity.py (no QPU needed for the analysis).

---

21. Measured Rust speedups vs. Python baseline

Re-run from tests/test_rust_path_benchmarks.py on 2026-04-17 (ML350 Gen8, scpn-quantum-engine 0.2.0, PyO3 0.29 + rayon 1.10). These are the only cross-language acceleration numbers we publish; they are measured, not estimated.

The Koopman rows were measured on 2026-05-12 with CPython 3.12.3 on x86_64 from a release wheel built by python -m maturin build --release --manifest-path scpn_quantum_engine/Cargo.toml --features extension-module --interpreter python. Each row reports the median of seven Python and Rust calls after a forced require_rust=True parity check against the NumPy generator.

Function	Python	Rust	Speedup
build_knm (16×16)	0.1 ms	0.01 ms	18.4×
koopman_generator (n=4, dim=16)	0.028014 ms	0.012752 ms	2.197×
koopman_generator (n=8, dim=64)	0.172246 ms	0.020980 ms	8.210×
koopman_generator (n=16, dim=256)	1.478324 ms	0.083449 ms	17.715×
kuramoto_euler (8 osc, 1 000 steps)	2.3 ms	0.25 ms	9.3×
correlation_matrix_xy (n=3)	0.7 ms	0.04 ms	19.5×
lindblad_jump_ops_coo (n=3)	0.0 ms	0.0 ms	9.2×
lindblad_anti_hermitian_diag (n=3)	0.0 ms	0.0 ms	4.7×

Standalone Rust paths (no Python parity comparison; absolute wall time):

Function	Configuration	Wall time
kuramoto_trajectory	16 osc, 10 000 steps	11.79 ms
expectation_pauli_fast	10 qubits, single-Z	0.02 ms
pec_sample_parallel	100 000 samples, 5 gates	14.04 ms
mc_xy_simulate	8 osc, 5k therm + 2k meas	1.39 ms
brute_mpc	dim=3, horizon=4	0.27 ms
parity_filter_mask	10 000 × 20-bit	0.14 ms

Cross-language outlook

These measured numbers are the verification baseline against which any future Mojo / Julia / Lean 4 investment is judged. Julia is now wired for order_parameter (see the next section); Mojo and Lean 4 are not yet benchmarked in this codebase.

Decision criteria for adopting a new acceleration backend:

1. Identify a specific compute-hot Python module that does not already have a Rust path. (Currently every module flagged in the internal Rust audit already has one — see docs/rust_engine.md.)
2. Build a minimal port and re-run the relevant section of tests/test_rust_path_benchmarks.py shape.
3. Publish the measured number in this table. Vague "X–Y faster" ranges are explicitly out of scope.

Python Orchestration Hot Paths

The table below is from code inspection plus the measured benchmark surfaces in this document. It identifies work that should move out of Python loops only after a parity test and benchmark row are added.

Current Python surface	Evidence	Target backend	Acceptance gate
analysis/koopman.py::build_koopman_generator	Python builds the dense n²×n² Koopman generator and now has a native scpn_quantum_engine.koopman_generator route.	Rust kernel returning the dense generator, with Python reconstructing canonical labels and retaining the NumPy fallback.	Parity on n=4/8/16 and a new Rust-vs-Python timing row.
benchmarks/quantum_advantage.py::quantum_benchmark	The benchmark constructs Qiskit circuits once, then repeatedly calls Statevector.evolve in a Python loop.	Vectorised simulator batch or Rust sparse statevector stepper once Qiskit semantics are matched.	Same final state/order-parameter tolerance for fixed K, omega, t_max, dt, and a provenance table row.
phase/mps_evolution.py::{dmrg_ground_state,tebd_evolution}	Python drives per-sweep and per-step tensor operations around quimb.	Prefer quimb-native vectorised calls first; only move glue to Rust if profiling shows Python overhead after tensor contraction cost is removed.	Quimb parity test with long-range-coupling caveat preserved and a TEBD timing table.
hardware/async_runner.py::_submit_blocking	The module is an I/O bridge and is explicitly exempt from the Rust-path rule; its hot path is repeated transpilation/submission orchestration, not numeric compute.	Keep in Python; batch or cache transpilation inputs before considering native code.	Hardware mock tests plus real-job provenance showing lower submission overhead without changing counts.

---

Multi-language accel chain

The dispatcher in src/scpn_quantum_control/accel/dispatcher.py forwards each compute function through an ordered chain of acceleration tiers (Rust → Julia → Python). The chain order must match measured wall-time on the runner class used in production. This section is the authoritative place where those measurements are recorded. The raw JSON is at docs/benchmarks/order_parameter_tiers.json and the reproducer is scripts/bench_order_parameter_tiers.py.

Consolidated tier benchmark (CI-canonical)

The per-primitive sections below are the original 2026-04-17 workstation snapshot (one reproducer script per primitive). The whole Kuramoto compute surface — every registered dispatcher, all three tiers — is now measured by a single consolidated runner that conforms to the unified benchmark standard (agentic-shared/BENCHMARK_STANDARD.md): warm-up plus repeats, the full P50 / P95 / P99 / mean / min / max set with throughput, one row per (operation, backend), explicit unavailable rows with their reason, and a full provenance block.

• Runner: scripts/bench_kuramoto_tiers.py (measurement core in src/scpn_quantum_control/accel/tier_benchmark.py).
• Canonical numbers: the fixed CI runner, which builds the Rust engine from source so every Rust tier is present. The workflow is .github/workflows/benchmark-kuramoto-tiers.yml (workflow_dispatch + weekly); the artefact is evidence-only because the hosted-runner CPU differs from the declared workstation.
• Side-by-side comparison: docs/tier_benchmarks.md (generated by tools/render_tier_benchmarks.py) places the CI and workstation columns next to each other for all primitives.
• Reproducibility manifest: docs/benchmarks/tiers/kuramoto_tiers_manifest.<env>.json records the host, toolchain versions (Python, NumPy, Rust engine, juliacall, rustc), commit, seed, run parameters, per-tier availability, and the fastest-backend / parity index for each primitive.

A workstation run against a stale installed engine wheel honestly marks the Rust tiers it lacks as unavailable; the CI column is the complete picture.

order_parameter(theta)

Measured 2026-04-17 on the local Linux runner (Intel i5-11600K, Python 3.12, juliacall 0.9.31, scpn-quantum-engine local build). Inner loop: 50 calls per sample × 5 samples; reported value is the per-call median.

N	Rust	Julia	Python
4	1.13 µs	11.19 µs	6.22 µs
16	0.90 µs	13.93 µs	5.92 µs
64	1.32 µs	13.32 µs	7.21 µs
256	2.97 µs	16.83 µs	12.82 µs
1024	13.12 µs	21.62 µs	26.60 µs
4096	38.10 µs	58.29 µs	123.89 µs
16384	256.50 µs	275.80 µs	465.78 µs

Read-offs:

• Rust wins at every measured N — the dispatcher places it first, which matches measurement.
• Julia is slower than Python for N ≤ 256 because the juliacall FFI crossing + Julia unboxing overhead cost more than the SIMD win on tiny arrays. The chain still places Julia before Python, which is the correct behaviour when Rust is unavailable and the caller expects some acceleration; callers with N < 256 who cannot use the Rust tier should pick the Python floor explicitly by importing _python_order_parameter directly.
• Julia beats Python from N ≥ 1024 onwards (1.23× at N = 1024, 2.13× at N = 4096, 1.69× at N = 16 384) — here the SIMD accumulator amortises the FFI cost.
• Rust's margin over Julia narrows as N grows (12.3× at N = 16, 1.08× at N = 16 384) — at very large N the bottleneck is memory bandwidth and both tiers converge to the same throughput.

Re-running

python scripts/bench_order_parameter_tiers.py

Default arguments: 7 outer repeats × 100 inner reps, sizes 4/16/64/256/1024/4096/16384. Output lands at docs/benchmarks/order_parameter_tiers.json. If the measured ordering changes (e.g. Julia overtakes Rust on a future runner with a smaller FFI overhead), update the chain comment at the top of _ORDER_PARAMETER_CHAIN in src/scpn_quantum_control/accel/order_parameter_observables.py to match.

order_parameter_gradient(theta)

Gradient of the order parameter, ∂r/∂θ_j = (1/N) sin(ψ − θ_j) with mean phase ψ. Measured 2026-06-22 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; reported value is the per-call median. Raw JSON: docs/benchmarks/order_parameter_gradient_tiers.json.

N	Rust	Julia	Python
4	1.18 µs	14.86 µs	18.27 µs
16	1.18 µs	13.76 µs	16.42 µs
64	2.12 µs	19.16 µs	22.87 µs
256	6.68 µs	37.40 µs	20.39 µs
1024	22.94 µs	41.92 µs	44.45 µs
4096	124.24 µs	284.09 µs	227.41 µs
16384	651.37 µs	904.09 µs	902.60 µs

Read-offs:

• Rust wins at every measured N (15.5× faster than Python at N = 4, 1.39× at N = 16 384), so the dispatcher places it first — matching measurement.
• Julia and the Python floor are close for the gradient: Julia is ahead for small N (the gradient's extra per-element trigonometry amortises the FFI crossing sooner than the scalar value does), but the Python floor overtakes at N = 256 and N = 4096 where NumPy's vectorised cos/sin is competitive. The chain still places Julia before the Python floor, consistent with the value chain: it is the "some acceleration when Rust is unavailable" default, and callers who cannot use the Rust tier and run mid-sized inputs should import _python_order_parameter_gradient directly.
• The gradient touches the same per-oscillator cos/sin work as the scalar value plus one fused multiply-add per element, so its tier ordering tracks the value chain (Rust → Julia → Python).

Re-run with:

python scripts/bench_order_parameter_gradient_tiers.py

If the measured ordering changes, update the chain comment at _ORDER_PARAMETER_GRADIENT_CHAIN in src/scpn_quantum_control/accel/order_parameter_observables.py to match.

order_parameter_hessian(theta)

Hessian of the order parameter, ∂²r/∂θ_i∂θ_j = a_i a_j/(N²r) − δ_ij a_j/N with alignment a_j = cos(ψ − θ_j). Symmetric, rows sum to zero. The result is an N × N matrix, so per-call cost grows quadratically and the benchmark stops at N = 2048 (~32 MB per matrix). Measured 2026-06-22 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/order_parameter_hessian_tiers.json.

N	Rust	Julia	Python
4	1.11 µs	16.17 µs	21.79 µs
16	1.56 µs	12.33 µs	22.24 µs
64	3.60 µs	75.54 µs	43.73 µs
256	42.99 µs	378.29 µs	203.12 µs
1024	2663.10 µs	7685.32 µs	8069.54 µs
2048	14237.46 µs	53647.20 µs	33803.40 µs

Read-offs:

• Rust wins at every measured N (19.6× faster than Python at N = 4, 2.4× at N = 2048); the dispatcher places it first, matching measurement.
• The Python floor beats Julia from N = 64 upward: NumPy builds the rank-one outer(a, a) matrix with a single vectorised BLAS-style kernel, whereas the juliacall crossing plus the per-element matrix fill dominate here. The chain still lists Julia before the Python floor for consistency with the value and gradient chains (the "some acceleration when Rust is unavailable" default); callers that cannot use the Rust tier and need the Hessian at N ≥ 64 should import _python_order_parameter_hessian directly.
• Cost is quadratic in N for every tier (the N × N matrix fill dominates the shared O(N) trigonometric pre-pass), so the gap to the scalar value and the vector gradient widens with N.

Re-run with:

python scripts/bench_order_parameter_hessian_tiers.py

If the measured ordering changes, update the chain comment at _ORDER_PARAMETER_HESSIAN_CHAIN in src/scpn_quantum_control/accel/order_parameter_observables.py to match.

mean_phase(theta) and mean_phase_gradient(theta)

Circular mean phase ψ = atan2(⟨sin θ⟩, ⟨cos θ⟩) — the collective phase of the complex order parameter z = r e^{iψ} — and its gradient ∂ψ/∂θ_j = cos(ψ − θ_j)/(N r) = (C cos θ_j + S sin θ_j)/(N r²), whose components sum to one. Both are O(N). Measured 2026-06-23 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/mean_phase_tiers.json.

mean_phase:

N	Rust	Julia	Python
4	0.80 µs	7.85 µs	7.28 µs
64	1.04 µs	7.39 µs	10.82 µs
1024	15.12 µs	18.48 µs	21.16 µs
16384	245.31 µs	286.42 µs	413.35 µs

mean_phase_gradient:

N	Rust	Julia	Python
4	0.78 µs	9.80 µs	11.48 µs
64	3.23 µs	12.03 µs	13.05 µs
1024	22.01 µs	36.89 µs	41.53 µs
16384	509.11 µs	596.23 µs	734.89 µs

Read-offs:

• Rust wins at every measured N for both routes; the dispatcher places it first, matching measurement.
• Julia and the Python floor are close — the mean-phase routes do the same per oscillator trigonometric pre-pass as the order parameter, so the FFI crossing and NumPy's vectorised reductions trade places much as they do there. The chain keeps Julia before the Python floor for consistency with the order-parameter chains.

Re-run with:

python scripts/bench_mean_phase_tiers.py

If the measured ordering changes, update the chain comment at _MEAN_PHASE_CHAIN and _MEAN_PHASE_GRADIENT_CHAIN in src/scpn_quantum_control/accel/mean_phase_observables.py to match.

mean_phase_hessian(theta)

Hessian ∂²ψ/∂θ_i∂θ_j = δ_ij s_j/(N r) − (s_i c_j + c_i s_j)/(N² r²) of the mean phase, with s_k = sin(ψ − θ_k), c_k = cos(ψ − θ_k). Symmetric, rows sum to zero. The result is an N × N matrix, so per-call cost grows quadratically and the benchmark stops at N = 2048. Measured 2026-06-23 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/mean_phase_hessian_tiers.json.

N	Rust	Julia	Python
4	1.54 µs	11.84 µs	28.28 µs
64	7.14 µs	43.37 µs	59.55 µs
256	44.04 µs	248.43 µs	255.57 µs
1024	1696.71 µs	5686.97 µs	9760.51 µs
2048	11849.49 µs	45918.33 µs	57458.32 µs

Read-offs:

• Rust wins at every measured N (18.4× faster than Python at N = 4, 4.8× at N = 2048); the dispatcher places it first, matching measurement.
• The mean-phase Hessian forms two rank-one outer products, so it is heavier per element than the order-parameter Hessian's single one; Rust's lead over the interpreted tiers therefore stays wider here than for the scalar and vector routes. Cost is quadratic in N for every tier.

Re-run with:

python scripts/bench_mean_phase_hessian_tiers.py

If the measured ordering changes, update the chain comment at _MEAN_PHASE_HESSIAN_CHAIN in src/scpn_quantum_control/accel/mean_phase_observables.py to match.

daido_order_parameter(theta, m) and daido_order_parameter_gradient(theta, m)

The m-th Daido order parameter r_m = |⟨exp(i m θ)⟩| detects m-cluster synchronisation (r_1 = 0 but r_m = 1 for m evenly spaced clusters), with gradient ∂r_m/∂θ_j = (m/N) sin(ψ_m − m θ_j). Both are O(N); for m = 1 they equal order_parameter and its gradient. Measured 2026-06-23 at the fixed harmonic m = 2 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/daido_order_parameter_tiers.json.

daido_order_parameter (m = 2):

N	Rust	Julia	Python
4	0.68 µs	7.57 µs	4.97 µs
64	1.25 µs	7.20 µs	5.61 µs
1024	10.04 µs	14.40 µs	21.79 µs
4096	67.80 µs	40.63 µs	96.12 µs
16384	289.61 µs	263.23 µs	413.67 µs

daido_order_parameter_gradient (m = 2):

N	Rust	Julia	Python
4	0.96 µs	8.35 µs	10.80 µs
64	2.17 µs	10.56 µs	13.02 µs
1024	20.52 µs	32.94 µs	42.35 µs
4096	129.95 µs	168.62 µs	179.64 µs
16384	579.87 µs	674.50 µs	868.64 µs

Read-offs:

• Rust wins the gradient at every N and the value through N = 1024. For the value at N ≥ 4096 Julia's BLAS-backed reduction edges ahead (40.63 vs 67.80 µs at N = 4096; 263 vs 290 µs at N = 16 384, within ~10%) — the extra m·θ multiply makes the value loop reduction-bound at large N. The chain still places Rust first: it wins every small-to-mid workload and the whole gradient route, and the large-N value gap is small. Callers running the value at N ≥ 4096 who have Julia installed can call _julia_daido_order_parameter directly.
• The value sums one complex exponential per oscillator; the gradient adds the per-oscillator sin(ψ_m − m θ_j) term, so Rust's lead is widest on the gradient.

Re-run with:

python scripts/bench_daido_order_parameter_tiers.py

If the measured ordering changes, update the chain comment at _DAIDO_ORDER_PARAMETER_CHAIN and _DAIDO_ORDER_PARAMETER_GRADIENT_CHAIN in src/scpn_quantum_control/accel/daido_observables.py to match.

daido_mode_phase(theta, m), its gradient and its Hessian

The phase of the m-th Fourier mode, ψ_m = atan2(⟨sin mθ⟩, ⟨cos mθ⟩) (the phase partner of the Daido magnitude r_m; for m = 1 it is the Kuramoto mean phase), its gradient ∂ψ_m/∂θ_j = (m/(N r_m)) cos(ψ_m − m θ_j) (components sum to m) and its Hessian H_ij = m²[δ_ij s_j/(N r_m) − (s_i c_j + c_i s_j)/(N² r_m²)] (symmetric, zero-row-sum; for m = 1 the mean-phase Hessian). The value and gradient are O(N); the Hessian is an N × N matrix, so it stops at N = 1024. Measured 2026-06-24 at the fixed harmonic m = 2 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/daido_mode_phase_tiers.json.

daido_mode_phase (m = 2):

N	Rust	Julia	Python
4	0.72 µs	9.03 µs	9.83 µs
64	1.41 µs	8.75 µs	8.39 µs
1024	12.67 µs	20.27 µs	25.79 µs
16384	357.03 µs	385.04 µs	551.87 µs

daido_mode_phase_gradient (m = 2):

N	Rust	Julia	Python
4	1.28 µs	13.50 µs	14.28 µs
64	2.87 µs	14.16 µs	16.19 µs
1024	25.64 µs	42.14 µs	51.56 µs
16384	694.88 µs	836.39 µs	1076.42 µs

daido_mode_phase_hessian (m = 2):

N	Rust	Julia	Python
4	1.26 µs	13.62 µs	24.12 µs
64	4.76 µs	20.69 µs	41.64 µs
256	37.03 µs	238.92 µs	283.44 µs
1024	1559.57 µs	6513.39 µs	9842.16 µs

Read-offs:

• Rust is first at every N for all three routes. The value and gradient are O(N) reductions where Julia closes to within ~10% on the value at N = 16384; the Hessian is O(N²) and Rust wins it decisively (4.2× faster than Julia and 6.3× faster than NumPy at N = 1024) because its fused outer-product loop avoids the materialised s⊗c temporaries. The chain places Rust first, matching measurement.

Re-run with:

python scripts/bench_daido_mode_phase_tiers.py

If the measured ordering changes, update the chain comments at _DAIDO_MODE_PHASE_CHAIN, _DAIDO_MODE_PHASE_GRADIENT_CHAIN and _DAIDO_MODE_PHASE_HESSIAN_CHAIN in src/scpn_quantum_control/accel/daido_phase.py to match.

daido_order_parameter_hessian(theta, m)

Hessian ∂²r_m/∂θ_i∂θ_j = m² (a_i a_j/(N² r_m) − δ_ij a_j/N) of the m-th Daido order parameter, with a_k = cos(ψ_m − m θ_k). Symmetric, rows sum to zero; for m = 1 it equals order_parameter_hessian. The result is an N × N matrix, so the benchmark stops at N = 2048. Measured 2026-06-23 at the fixed harmonic m = 2 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/daido_order_parameter_hessian_tiers.json.

N	Rust	Julia	Python
4	1.02 µs	9.17 µs	19.37 µs
64	3.20 µs	24.49 µs	32.75 µs
256	27.02 µs	201.31 µs	149.48 µs
1024	937.20 µs	5106.65 µs	4113.88 µs
2048	9636.63 µs	38547.01 µs	30616.84 µs

Read-offs:

• Rust wins at every measured N (19.0× faster than Python at N = 4, 3.2× at N = 2048); the dispatcher places it first, matching measurement. The Daido Hessian carries the same rank-one structure as the order-parameter Hessian plus the extra m·θ multiply per element, so the per-call cost is quadratic in N for every tier and Rust's lead is widest here.

Re-run with:

python scripts/bench_daido_order_parameter_hessian_tiers.py

If the measured ordering changes, update the chain comment at _DAIDO_ORDER_PARAMETER_HESSIAN_CHAIN in src/scpn_quantum_control/accel/daido_observables.py to match.

mean_field_force(theta, K) and mean_field_jacobian(theta, K)

The Kuramoto mean-field force F_j = K (S cos θ_j − C sin θ_j) is the phase-coupling term of the dynamics θ̇_j = ω_j + F_j; its stability Jacobian J_jk = (K/N) cos(θ_j − θ_k) − K δ_jk (C cos θ_j + S sin θ_j) is symmetric with a zero row sum (the global-phase Goldstone mode). The force is O(N); the Jacobian is an N × N matrix, so the benchmark stops at N = 2048. Measured 2026-06-23 at the fixed coupling K = 1 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/mean_field_tiers.json.

mean_field_force (K = 1):

N	Rust	Julia	Python
4	1.11 µs	14.66 µs	11.36 µs
64	3.52 µs	19.99 µs	18.25 µs
1024	35.34 µs	51.55 µs	50.73 µs
2048	72.27 µs	110.81 µs	97.58 µs

mean_field_jacobian (K = 1):

N	Rust	Julia	Python
4	2.38 µs	25.26 µs	18.88 µs
64	41.02 µs	59.42 µs	64.66 µs
256	814.42 µs	963.25 µs	999.19 µs
2048	65468.84 µs	88843.01 µs	128051.96 µs

Read-offs:

• Rust is first for both routes at every N, matching the chain ordering. The force lead is the FFI-vs-vectorised-NumPy gap (~10× at small N narrowing to ~1.3× at N = 2048 where the work is bandwidth-bound); the Jacobian builds the full pairwise cos(θ_j − θ_k) matrix, where Rust's elementwise loop avoids NumPy's temporary broadcast and stays ~2× ahead at N = 2048.

Re-run with:

python scripts/bench_mean_field_tiers.py

If the measured ordering changes, update the chain comments at _MEAN_FIELD_FORCE_CHAIN and _MEAN_FIELD_JACOBIAN_CHAIN in src/scpn_quantum_control/accel/kuramoto_mean_field.py to match.

daido_mean_field_force(theta, K, m) and daido_mean_field_jacobian(theta, K, m)

The Daido m-th-harmonic mean field couples each oscillator to the m-th Fourier mode of the ensemble, F_j = K (S_m cos m θ_j − C_m sin m θ_j) = K r_m sin(ψ_m − m θ_j) (the all-to-all mean field at m = 1), driving m-cluster synchronisation; its stability Jacobian J_jl = K m [(1/N) cos(m(θ_j − θ_l)) − δ_jl r_m cos(ψ_m − m θ_j)] is symmetric with a zero row sum (the global-phase Goldstone mode). The force is O(N); the Jacobian is an N × N matrix, so the benchmark stops at N = 2048. Measured 2026-06-24 at the fixed coupling K = 1 and harmonic m = 2 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/daido_mean_field_tiers.json.

daido_mean_field_force (K = 1, m = 2):

N	Rust	Julia	Python
4	1.01 µs	10.96 µs	10.51 µs
64	2.24 µs	10.25 µs	12.81 µs
1024	21.51 µs	35.02 µs	40.99 µs
2048	43.99 µs	125.11 µs	69.25 µs

daido_mean_field_jacobian (K = 1, m = 2):

N	Rust	Julia	Python
4	1.23 µs	8.76 µs	14.80 µs
64	31.70 µs	140.87 µs	51.28 µs
256	823.90 µs	1067.48 µs	960.32 µs
2048	69508.77 µs	85901.52 µs	94583.11 µs

Read-offs:

• Rust is first for both routes at every N, matching the chain ordering. The force is the m-scaled analogue of the mean-field force, with the same FFI-vs-vectorised-NumPy lead; the Jacobian builds the full pairwise cos(m(θ_j − θ_l)) matrix in a fused loop and stays ahead of both the Julia tier and NumPy's broadcast temporary at N = 2048.

Re-run with:

python scripts/bench_daido_mean_field_tiers.py

If the measured ordering changes, update the chain comments at _DAIDO_MEAN_FIELD_FORCE_CHAIN and _DAIDO_MEAN_FIELD_JACOBIAN_CHAIN in src/scpn_quantum_control/accel/daido_mean_field.py to match.

sakaguchi_mean_field_force(theta, K, α) and sakaguchi_mean_field_jacobian(theta, K, α)

The Sakaguchi–Kuramoto mean field adds a phase frustration α to the all-to-all coupling, F_j = K r sin(ψ − θ_j − α) (the frustrated travelling-wave drive); its stability Jacobian J_jl = (K/N) cos(θ_j − θ_l + α) − δ_jl K r cos(ψ − θ_j − α) is non-symmetric for α ≠ 0 (the dynamics are non-variational) yet keeps a zero row sum (the global-phase Goldstone mode). The force is O(N); the Jacobian is an N × N matrix, so the benchmark stops at N = 2048. Measured 2026-06-24 at the fixed coupling K = 1 and frustration α = 0.5 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/sakaguchi_mean_field_tiers.json.

sakaguchi_mean_field_force (K = 1, α = 0.5):

N	Rust	Julia	Python
4	1.01 µs	9.65 µs	14.28 µs
64	2.18 µs	13.14 µs	15.81 µs
1024	20.42 µs	36.27 µs	44.56 µs
2048	38.96 µs	75.27 µs	92.30 µs

sakaguchi_mean_field_jacobian (K = 1, α = 0.5):

N	Rust	Julia	Python
4	1.29 µs	11.08 µs	16.79 µs
64	32.72 µs	54.29 µs	53.25 µs
256	753.43 µs	823.70 µs	867.11 µs
2048	66737.71 µs	71708.45 µs	93044.18 µs

Read-offs:

• Rust is first for both routes at every N, matching the chain ordering. The force is the frustrated analogue of the mean-field force, with the same FFI-vs-vectorised-NumPy lead; the Jacobian builds the full cos(θ_j − θ_l + α) matrix in a fused loop and stays ahead of both the Julia tier and NumPy's broadcast temporary at N = 2048.

Re-run with:

python scripts/bench_sakaguchi_mean_field_tiers.py

If the measured ordering changes, update the chain comments at _SAKAGUCHI_MEAN_FIELD_FORCE_CHAIN and _SAKAGUCHI_MEAN_FIELD_JACOBIAN_CHAIN in src/scpn_quantum_control/accel/sakaguchi_mean_field.py to match.

triadic_mean_field_force(theta, K) and triadic_mean_field_jacobian(theta, K)

The triadic (2-simplex) higher-order Kuramoto mean field couples each oscillator through the squared first moment, F_j = K r² sin(2ψ − 2θ_j) (the Skardal–Arenas drive behind explosive synchronisation, distinct from the second Daido harmonic, which uses the second moment ⟨e^{2iθ}⟩ rather than ⟨e^{iθ}⟩²); its Jacobian J_jl = (2K/N) r cos(2θ_j − θ_l − ψ) − δ_jl 2K r² cos(2ψ − 2θ_j) is non-symmetric (the higher-order mean field is non-variational) yet keeps a zero row sum (the global-phase Goldstone mode). The force is O(N); the Jacobian is an N × N matrix, so the benchmark stops at N = 2048. Measured 2026-06-24 at the fixed coupling K = 1 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 100 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/triadic_mean_field_tiers.json.

triadic_mean_field_force (K = 1):

N	Rust	Julia	Python
4	1.94 µs	9.44 µs	10.74 µs
64	2.25 µs	10.92 µs	13.40 µs
1024	22.81 µs	41.70 µs	45.45 µs
2048	41.43 µs	80.32 µs	95.08 µs

triadic_mean_field_jacobian (K = 1):

N	Rust	Julia	Python
4	1.31 µs	11.75 µs	19.47 µs
64	51.12 µs	83.50 µs	117.88 µs
256	1219.62 µs	1268.54 µs	1938.60 µs
2048	92926.71 µs	97702.08 µs	179851.22 µs

Read-offs:

• Rust is first for both routes at every N, matching the chain ordering. The force is an O(N) reduction with the usual FFI-vs-vectorised-NumPy lead; the Jacobian builds the full cos(2θ_j − θ_l) / sin(2θ_j − θ_l) matrix in a fused loop and stays ahead of both the Julia tier and NumPy's broadcast temporary (≈1.9× faster than NumPy at N = 2048).

Re-run with:

python scripts/bench_triadic_mean_field_tiers.py

If the measured ordering changes, update the chain comments at _TRIADIC_MEAN_FIELD_FORCE_CHAIN and _TRIADIC_MEAN_FIELD_JACOBIAN_CHAIN in src/scpn_quantum_control/accel/triadic_mean_field.py to match.

networked_kuramoto_force(theta, K) and networked_kuramoto_jacobian(theta, K)

The general (graph) Kuramoto force F_j = Σ_k K_jk sin(θ_k − θ_j) for a coupling matrix K and its stability Jacobian (J_jl = K_jl cos(θ_l − θ_j), l ≠ j; J_jj = −Σ_{k≠j} K_jk cos(θ_k − θ_j)). Both routes are O(N²) — the force sums over the full coupling row — so the sizes stop at N = 2048. Measured 2026-06-23 on a dense symmetric random coupling matrix on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 50 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/networked_kuramoto_tiers.json.

networked_kuramoto_force:

N	Rust	Julia	Python
4	2.56 µs	19.29 µs	4.96 µs
64	31.69 µs	48.40 µs	53.82 µs
256	825.17 µs	739.25 µs	915.51 µs
2048	99005.61 µs	87903.85 µs	100478.17 µs

networked_kuramoto_jacobian:

N	Rust	Julia	Python
4	1.58 µs	276.80 µs	14.58 µs
64	60.96 µs	78.91 µs	54.11 µs
256	968.56 µs	1724.86 µs	1043.65 µs
2048	71741.98 µs	77389.81 µs	110316.76 µs

Read-offs:

• Rust is first at small-to-mid N for both routes and wins the Jacobian at N = 2048. For the force, Julia's BLAS-backed reduction edges ahead from N = 256 (739 vs 825 µs) to N = 2048 (88 vs 99 ms, within ~11%); the chain still places Rust first because it wins every small workload and the whole Jacobian route, and the large-N force gap is small. Callers running the force at N ≥ 256 with Julia installed can call _julia_networked_kuramoto_force directly.
• Julia's Jacobian is dominated by FFI overhead at small N (277 µs at N = 4 to marshal the N × N matrix), so it only becomes competitive past N ≈ 1024; Rust and the vectorised NumPy floor both avoid that fixed cost.

Re-run with:

python scripts/bench_networked_kuramoto_tiers.py

If the measured ordering changes, update the chain comments at _NETWORKED_KURAMOTO_FORCE_CHAIN and _NETWORKED_KURAMOTO_JACOBIAN_CHAIN in src/scpn_quantum_control/accel/networked_kuramoto.py to match.

kuramoto_interaction_energy(theta, K), its gradient and its Hessian

The Kuramoto interaction energy E = −½ Σ_jk K_jk cos(θ_j − θ_k) (the Lyapunov function whose gradient flow is the dynamics for symmetric K), its gradient ∂E/∂θ_j = ½ Σ_k (K_jk + K_kj) sin(θ_j − θ_k) (the symmetrised force, components sum to zero) and its Hessian ∂²E/∂θ_i∂θ_l = −½(K_il + K_li) cos(θ_i − θ_l) for l ≠ i (symmetric, zero-row-sum; for symmetric K it equals the negated networked Jacobian). All three are O(N²), so the sizes stop at N = 2048. Measured 2026-06-24 on a dense symmetric random coupling matrix on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 50 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/kuramoto_energy_tiers.json.

kuramoto_interaction_energy:

N	Rust	Julia	Python
4	1.37 µs	23.90 µs	6.23 µs
64	46.84 µs	65.07 µs	62.29 µs
256	888.33 µs	777.90 µs	1253.15 µs
2048	83679.45 µs	65816.60 µs	102431.09 µs

kuramoto_interaction_energy_gradient:

N	Rust	Julia	Python
4	1.43 µs	169.25 µs	11.98 µs
64	33.26 µs	143.04 µs	47.38 µs
256	959.01 µs	1212.74 µs	1143.13 µs
2048	130928.77 µs	118233.64 µs	125975.93 µs

kuramoto_interaction_energy_hessian:

N	Rust	Julia	Python
4	1.97 µs	27.85 µs	11.02 µs
64	37.76 µs	59.99 µs	47.50 µs
256	950.62 µs	956.24 µs	1035.37 µs
2048	129571.12 µs	144461.84 µs	160374.82 µs

Read-offs:

• Rust is first at small-to-mid N for all three routes and wins the Hessian at N = 2048 (130 vs 144 ms vs Julia). Julia's BLAS reduction edges the scalar energy ahead from N = 256 (778 vs 888 µs) to N = 2048 (66 vs 84 ms), and edges the gradient at N = 2048 (118 vs 131 ms); the chain keeps Rust first because it wins every small workload and the Hessian, and the large-N crossovers are small.
• Julia's matrix-marshalling cost dominates the gradient and Hessian at small N (169/28 µs at N = 4), so Rust and the vectorised NumPy floor lead until N ≈ 1024.

Re-run with:

python scripts/bench_kuramoto_energy_tiers.py

If the measured ordering changes, update the chain comments at _KURAMOTO_INTERACTION_ENERGY_CHAIN, _KURAMOTO_INTERACTION_ENERGY_GRADIENT_CHAIN and _KURAMOTO_INTERACTION_ENERGY_HESSIAN_CHAIN in src/scpn_quantum_control/accel/kuramoto_energy.py to match.

sakaguchi_force(theta, K, α) and sakaguchi_jacobian(theta, K, α)

The Kuramoto–Sakaguchi frustrated force F_j = Σ_{k≠j} K_jk sin(θ_k − θ_j − α) and its stability Jacobian (J_jl = K_jl cos(θ_l − θ_j − α), l ≠ j; J_jj = −Σ_{k≠j} K_jk cos(θ_k − θ_j − α)). The frustration α breaks reciprocity, so the Jacobian is asymmetric for α ≠ 0. Both routes are O(N²), so the sizes stop at N = 2048. Measured 2026-06-23 on a dense symmetric random coupling matrix at α = 0.3 on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 50 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/sakaguchi_kuramoto_tiers.json.

sakaguchi_force (α = 0.3):

N	Rust	Julia	Python
4	1.35 µs	16.01 µs	10.07 µs
64	38.04 µs	47.10 µs	46.91 µs
256	907.02 µs	647.75 µs	1105.99 µs
2048	66216.15 µs	59439.56 µs	92195.17 µs

sakaguchi_jacobian (α = 0.3):

N	Rust	Julia	Python
4	2.45 µs	17.18 µs	6.67 µs
64	35.03 µs	56.54 µs	42.25 µs
256	977.16 µs	900.20 µs	945.40 µs
2048	85922.89 µs	85880.37 µs	119125.31 µs

Read-offs:

• Rust is first at small-to-mid N for both routes. As with the networked force, Julia's BLAS-backed reduction edges the force ahead from N = 256 (648 vs 907 µs) to N = 2048 (59 vs 66 ms); the Jacobian is a near-tie at large N (Rust and Julia within 0.05% at N = 2048). The chain keeps Rust first because it wins every small workload, and the large-N gaps are small. The frustration angle adds no measurable cost over the networked kernels — it is a constant offset inside the existing trigonometric argument.

Re-run with:

python scripts/bench_sakaguchi_kuramoto_tiers.py

If the measured ordering changes, update the chain comments at _SAKAGUCHI_FORCE_CHAIN and _SAKAGUCHI_JACOBIAN_CHAIN in src/scpn_quantum_control/accel/sakaguchi_kuramoto.py to match.

local_order_parameter(theta, A) and local_order_parameter_jacobian(theta, A)

The network-local order parameter r_j = |Σ_k A_jk e^{iθ_k}| / Σ_k A_jk (the degree- normalised local coherence, the standard chimera diagnostic) and its Jacobian ∂r_j/∂θ_l = (A_jl/d_j) sin(ψ_j − θ_l). Both routes are O(N²) — each node sums over its full adjacency row — so the sizes stop at N = 2048. Measured 2026-06-23 on a dense random adjacency matrix on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 50 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/local_order_parameter_tiers.json.

local_order_parameter:

N	Rust	Julia	Python
4	2.05 µs	49.34 µs	17.66 µs
64	11.50 µs	43.40 µs	15.17 µs
1024	1563.27 µs	1701.95 µs	15253.64 µs
2048	6625.78 µs	6680.90 µs	22416.56 µs

local_order_parameter_jacobian:

N	Rust	Julia	Python
4	2.96 µs	132.39 µs	28.97 µs
64	14.16 µs	80.50 µs	62.46 µs
1024	4139.15 µs	10770.38 µs	35449.13 µs
2048	25476.57 µs	74026.61 µs	153172.53 µs

Read-offs:

• Rust is first at every N for both routes. The value is a near-tie with Julia at large N (6.63 vs 6.68 ms at N = 2048), but Rust wins the Jacobian decisively (2.9× faster than Julia and 6× faster than NumPy at N = 2048): the floor builds the full N × N sin(ψ_j − θ_l) basis as a temporary, whereas Rust accumulates each node's degree, local mean phase and Jacobian row in a single fused pass. The Python value path stays competitive to N ≈ 256 (its A @ cos θ reduction is BLAS) but the per-node division and the materialised basis cost it the large-N range.

Re-run with:

python scripts/bench_local_order_parameter_tiers.py

If the measured ordering changes, update the chain comments at _LOCAL_ORDER_PARAMETER_CHAIN and _LOCAL_ORDER_PARAMETER_JACOBIAN_CHAIN in src/scpn_quantum_control/accel/local_order.py to match.

local_mean_phase(theta, A) and local_mean_phase_jacobian(theta, A)

The network-local mean phase ψ_j = atan2(Σ_k A_jk sin θ_k, Σ_k A_jk cos θ_k) — the argument of the local complex order Z_j = Σ_k A_jk e^{iθ_k} and the phase partner of the local order parameter r_j (together Z_j/d_j = r_j e^{iψ_j}) — and its Jacobian ∂ψ_j/∂θ_l = A_jl cos(ψ_j − θ_l) / |Z_j|. Both routes are O(N²) — each node sums over its full adjacency row — so the sizes stop at N = 2048. Measured 2026-06-24 on a dense random adjacency matrix on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 50 calls per sample × 7 samples; per-call median. Raw JSON: docs/benchmarks/local_mean_phase_tiers.json.

local_mean_phase:

N	Rust	Julia	Python
4	3.24 µs	29.22 µs	13.15 µs
64	11.57 µs	36.57 µs	17.29 µs
256	96.91 µs	129.55 µs	50.50 µs
2048	5850.38 µs	6637.53 µs	13492.91 µs

local_mean_phase_jacobian:

N	Rust	Julia	Python
4	2.94 µs	31.40 µs	30.99 µs
64	12.36 µs	39.98 µs	32.22 µs
256	160.76 µs	427.28 µs	440.01 µs
2048	26838.08 µs	62790.56 µs	125376.61 µs

Read-offs:

• Rust is first at every N for the Jacobian (2.3× faster than Julia and 4.7× faster than NumPy at N = 2048): Rust accumulates each node's local complex order and Jacobian row in a single fused pass, where the floor materialises the full cos(ψ_j − θ_l) basis. For the value the Python floor's A @ cos θ reduction is BLAS and edges ahead in the mid-N band (≈50 vs 97 µs at N = 256), but the per-node atan2 and the masking cost it the large-N range, so the chain keeps Rust first — fastest at the small and large N that dominate use.

Re-run with:

python scripts/bench_local_mean_phase_tiers.py

If the measured ordering changes, update the chain comments at _LOCAL_MEAN_PHASE_CHAIN and _LOCAL_MEAN_PHASE_JACOBIAN_CHAIN in src/scpn_quantum_control/accel/local_phase.py to match.

kuramoto_euler_trajectory(...) and kuramoto_euler_vjp(...)

The differentiable networked-Kuramoto Euler integrator: the forward θ_{n+1} = θ_n + dt(ω + F(θ_n)) with F_j = Σ_k K_jk sin(θ_k − θ_j) recording the full trajectory, and its reverse-mode adjoint λ_n = λ_{n+1} + dt J(θ_n)ᵀ λ_{n+1} returning ∂L/∂{θ₀, ω, K} for a terminal objective L. Both routes are O(n_steps · N²); measured 2026-06-24 at a fixed step budget n_steps = 50 on a dense random coupling matrix on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 20 calls per sample × 7 samples; per-call median; sizes stop at N = 512. Raw JSON: docs/benchmarks/diff_kuramoto_euler_tiers.json.

kuramoto_euler_trajectory (n_steps = 50):

N	Rust	Julia	Python
4	7.49 µs	25.11 µs	295.91 µs
64	1736.83 µs	1600.64 µs	2280.07 µs
256	32278.06 µs	32304.33 µs	39970.61 µs
512	180532.94 µs	239200.05 µs	279290.74 µs

kuramoto_euler_vjp (n_steps = 50):

N	Rust	Julia	Python
4	34.03 µs	125.09 µs	688.89 µs
64	2932.06 µs	3041.90 µs	5610.73 µs
256	73920.12 µs	71299.27 µs	77958.85 µs
512	391885.99 µs	388681.18 µs	427163.48 µs

Read-offs:

• Rust is decisively first at small N (3–4× faster than Julia at N = 4, where its allocation-free flat-slice loop pays off against per-call overhead) and is the chain's first tier. Both Rust and Julia compile to tight O(n_steps · N²) loops, so at large N they converge to within a few percent and trade the lead run-to-run (e.g. the forward at N = 512 favours Rust, the adjoint at N = 256/512 favours Julia by ≲ 4%, inside the host-noise band of a non-isolated workstation). Both stay well ahead of the NumPy floor, which pays Python-loop overhead per step.

Re-run with:

python scripts/bench_diff_kuramoto_euler_tiers.py

If the measured ordering changes, update the chain comments at _KURAMOTO_EULER_TRAJECTORY_CHAIN and _KURAMOTO_EULER_VJP_CHAIN in src/scpn_quantum_control/accel/diff_kuramoto_euler.py to match.

kuramoto_rk4_trajectory(...) and kuramoto_rk4_vjp(...)

The differentiable networked-Kuramoto fourth-order Runge–Kutta integrator: the forward θ_{n+1} = θ_n + (dt/6)(k1 + 2k2 + 2k3 + k4) (four stages of f = ω + F) recording the full trajectory, and its reverse-mode adjoint backpropagating a terminal cotangent through the four stages to return ∂L/∂{θ₀, ω, K}. Fourth-order accurate where the Euler integrator is first-order, at four force evaluations per step. Both routes are O(n_steps · N²); measured 2026-06-24 at a fixed step budget n_steps = 50 on a dense random coupling matrix on the local Linux runner (Intel i5-11600K @ 3.90 GHz, Python 3.12.3, juliacall, scpn-quantum-engine local build). Inner loop: 20 calls per sample × 7 samples; per-call median; sizes stop at N = 512. Raw JSON: docs/benchmarks/diff_kuramoto_rk4_tiers.json.

kuramoto_rk4_trajectory (n_steps = 50):

N	Rust	Julia	Python
4	25.48 µs	38.67 µs	1173.61 µs
64	4990.44 µs	4961.62 µs	6561.66 µs
256	62255.51 µs	72447.83 µs	85494.35 µs
512	577108.01 µs	363119.40 µs	682349.46 µs

kuramoto_rk4_vjp (n_steps = 50):

N	Rust	Julia	Python
4	66.30 µs	645.83 µs	3290.90 µs
64	13126.00 µs	55235.22 µs	21583.15 µs
256	190324.76 µs	845643.18 µs	245770.48 µs
512	1462893.30 µs	4408806.10 µs	1663594.63 µs

Read-offs:

• Rust is first for the adjoint at every N — decisively so (3× faster than the Julia tier and ahead of NumPy at N = 512): its buffer-reused four-stage backward pass stays allocation-free, where the Julia VJP allocates the stage cotangents each step. The forward trades the lead with Julia (both compile to four O(N²) stage loops; the forward favours Julia at the largest N, inside the host-noise band), and the chain keeps Rust first — fastest on the adjoint, which is the gradient path that dominates optimisation. The NumPy floor stays the slowest tier.

Re-run with:

python scripts/bench_diff_kuramoto_rk4_tiers.py

If the measured ordering changes, update the chain comments at _KURAMOTO_RK4_TRAJECTORY_CHAIN and _KURAMOTO_RK4_VJP_CHAIN in src/scpn_quantum_control/accel/diff_kuramoto_rk4.py to match.

S4 multi-hardware readiness

Command:

scpn-bench s4-multi-hardware-ready

Regenerates data/s4_multi_hardware_control/s4_multi_hardware_readiness_2026-05-06.json and docs/campaigns/s4_multi_hardware_readiness_2026-05-06.md. The harness compiles the Kuramoto-XY programme into neutral-atom and circuit-QED analogue forms, exports Pulser, Bloqade, and IBM pulse-level provider payloads, and wraps each route in an approval-gated execution plan. It performs no provider contact, no emulator execution, and no QPU submission.

S4 IBM pulse-level preregistration

Command:

scpn-bench s4-provider-preregistration

Regenerates data/s4_multi_hardware_control/s4_ibm_pulse_preregistration_2026-05-06.json and docs/campaigns/s4_ibm_pulse_preregistration_2026-05-06.md. The harness reads the S4 readiness payload and packages the IBM pulse-level route as a calibration-review dossier. It performs no provider contact, creates no calibrated pulse schedule, and requests no QPU time.

S4 neutral-atom preregistration

Command:

scpn-bench s4-neutral-atom-preregistration

Regenerates data/s4_multi_hardware_control/s4_neutral_atom_preregistration_2026-05-06.json and docs/campaigns/s4_neutral_atom_preregistration_2026-05-06.md. The harness reads the S4 readiness payload and packages the Pulser/Bloqade neutral-atom routes as provider-object construction review dossiers. It performs no provider SDK construction, no emulator run, no cloud contact, and no QPU submission.

S5 benchmark suite

Command:

scpn-bench s5-benchmark-suite

Regenerates data/s5_benchmark_harness/phase1_benchmark_harness_2026-05-06.json and docs/campaigns/benchmark_harness_phase1_2026-05-06.md. The harness exposes the Phase 1 DLA-parity raw-count dataset through scpn_quantum_control.benchmark_harness, recomputes the published statistics, checks the noiseless classical parity-conservation baseline, and performs no QPU submission.

S5 benchmark registry

Command:

scpn-bench s5-benchmark-registry

Regenerates data/s5_benchmark_harness/benchmark_registry_2026-05-06.json and docs/campaigns/benchmark_harness_registry_2026-05-06.md. The registry lists implemented and planned benchmark families separately so CHSH, BKT, OTOC, and DLA-dimension plans are visible without being presented as available benchmark results.

S6 quantum-kuramoto split audit

Command:

scpn-bench s6-split-audit

Regenerates data/s6_quantum_kuramoto_split/quantum_kuramoto_split_audit_2026-05-07.json and an internal split-audit note. The audit classifies candidate phase, bridge, hardware, and accel modules as reusable, needs-review, or SCPN-specific. It does not create or publish a second package.

S6 quantum-kuramoto boundary review

Command:

scpn-bench s6-boundary-review

Regenerates data/s6_quantum_kuramoto_split/quantum_kuramoto_boundary_review_2026-05-07.json and docs/campaigns/quantum_kuramoto_boundary_review_2026-05-07.md. The review proposes a stable public API surface, documents decisions for needs-review rows, and keeps package-skeleton creation blocked until config/provenance/analysis-dependent refactors are closed.