Parameter derivatives

Analytic derivatives of the waveform with respect to mass_ratio, spin1z and spin2z. Written 2026-09-02, when they first covered the same span as the waveform.

What you get

times, modes, derivatives = model.get_td_waveform_modes_geometric(
    mass_ratio, spin1z, spin2z, with_derivatives=True)

over t in [-3819.1, +124.4] M at mass_ratio = 4, spins (0.3, -0.2), which is the SAME span the waveform covers, 15775 samples. include_merger_ringdown defaults to "auto", matching the waveform call, so the two sibling calls agree by default; with include_merger_ringdown=False the derivatives are the inspiral alone, 15276 samples from -3819.1 M to -0.3 M (measured 2026-09-02).

Since checklist A7 the derivatives ride a flag on the geometric mode method rather than a method of their own, and behind it _modes_and_derivatives goes straight to the fused gradient kernel; see waveform-evaluation.md for why it is not built on top of the waveform call. model.modes_grad(...) is the superseded name and is retained, forwarding to the same implementation.

A mode subset and a time_step are REFUSED with derivatives rather than silently dropped: the fused kernel produces all seven modes on its own grid, and honouring neither argument while accepting it would answer a different question from the one asked. Before 2026-09-02 the derivatives stopped at the amplitude peak while the waveform ran through ringdown, which meant a Fisher integral crossing the join saw the derivative fall to zero.

The waveform comes free. aa_gradients.py:1161-1165 builds the derivative from cP[p] and rot, and the product of those two IS the waveform, so forming it is one complex multiply per mode per sample. You cannot obtain the derivative without obtaining the waveform's ingredients, which is why the combined return exists rather than a derivatives-only path being cheaper.

Batched

modes_grad_batch(model, thetas, times, n_procs=..., with_mr="auto")

Bitwise against the scalar route, not merely close: 0.0 across 5 parameter sets, 7 modes, 3 parameters and 12935 samples, on both the waveform and the derivative, with the fork pool bitwise against a serial loop. So the scalar route's validation transfers exactly.

Process guidance, measured: about 32 processes below B = 256, about 60 above B = 512, never 128. The turnover is NOT compute. The fork pool is created inside every call, costing 0.16 s at 32 processes, 1.43 s at 60 and 2.51 s at 120, and Pool.map pipes 3.7 MB per row.

Two costs that were measured rather than assumed. The recorded v3 batched anti-scaling defect does NOT reach this route, shown by weak scaling holding flat from 8 to 60 workers, because that defect is in the compiled batched VALUE evaluator while this is NumPy processes at one thread each. And hoisting the already-batched merger-ringdown chain was measured and rejected: its theta-only part is 1.8 percent of a row, so the hoist would add 11.1 s to remove 0.14 s at B = 1024.

The seam

The derivative is as continuous at the inspiral to merger-ringdown join as the waveform is. The jump estimator falls as dt cubed, ratio 7.99 per halving, 1.17e-2 down to 3.05e-6 over five refinements, tracking the value blend's own floor to within a few percent.

Underneath that, and this is model physics rather than code: a HARD join would expose a 0.43x to 6.94x arm-to-arm derivative mismatch, worst on (2,2) at 5.42x and almost independent of theta. The ramp spreads that over 15 M. It does not remove it, and the delivered envelope carries 20 to 60 times the curvature of the waveform there as a result.

Accuracy

Finite differences against an independent value path, on the merger-ringdown window specifically: worst 2.54e-3 on (3,3) d/dchi2z at t = +18.4 M, q = 6. chi1z is the best channel at 3.7e-5. chi2z is worst everywhere and the worst times cluster at +13 to +31 M, not in the late tail.

The ramp first read 0.40, and that was attributed rather than accepted: the same comparison BELOW the ramp, on untouched inspiral-only code, reads 0.58, which is the documented blend-weight pathology, and with weights linearized the inspiral reads 2.9e-5, the ramp 2.2e-5 and the pure merger-ringdown 2.5e-3. The join contributes no error of its own.

The clock term matters more than its size suggests. dtau_end/dtheta is only 2e-3 to 3.9e-2 M, but carrying dh_MR/dt . dtau_end and w'(t) . dtau_end is worth an accuracy ratio of 817 in the ramp, 1.81e-2 down to 2.22e-5, for 3.3 ms.

Limits, stated where you will need them

• The device route now covers the merger and ringdown. cpp/mr.cpp gained mr_segment_grad and imr_join_grad, and the host theta-only stage that feeds them is compiled too. 6.894 ms per waveform at B = 1024 on 60 threads, 1.697x over the host-parallel route, at 1.00x accuracy. Where the time now is: the merger-ringdown arm is 2.2 percent of the route and the inspiral arm is 97.8 percent, so further work belongs there. On agreement: the compiled contraction differs from the NumPy reference by 1.9e-10 and CANNOT do better, because the contraction cancels by 8.3e6, which is cond(alpha) = 6.5e16 arriving at evaluation time. The reference is not accurate to 1e-12 either. Arbitrated against a float128 sum of the same terms the two are indistinguishable, 0.59x to 1.42x over eight channels with the compiled route ahead in five.
• Hybridised batch derivatives are inspiral-only, now by explicit pin. Extending them needs the scalar hybrid chain extended first.
• modes_grad_batch does not open _artifact_arm(), so it serves under the environment's PN normaliser rather than the artifact's. Inert on the shipped v3 artifact, since stamp and live flag agree, which is exactly why every bitwise gate is blind to it. It would bite on a pre-adoption artifact.
• evaluate_dt was not delegated alongside evaluate_grad, so two expressions for dh/dt exist with nothing gating them against each other.
• The gradient chain does not implement the correlation odd kernel. nrhjsurrogate/gradients.py:291 and nrhjsurrogate/driver/aa_gradients.py:665 build the parity kernel inline as 0.5 * (k1 -/+ k2) and never apply the variance normalisation that gpr_kernel_basis._sym_rbf applies when the correlation form is on. A model fitted with odd_corr_kernel true would therefore get values from one kernel and derivatives from another. Inert on every released artifact, all of which store odd_corr_kernel false (measured on every group of AA/AdA/MR v1, v2 and v3, 2026-09-02), and the flag is now persisted so the disagreement can be read off the file rather than inferred.
• The AdA arm carries full-IMR gradients, but only through the compiled module.

Where the code is