madspace: stable 2->3 Jacobian when a particle is soft

[oliviermattelaer]

Problem

ColorOrderedMapping occasionally returns weights up to ~1e20 times the typical one, which can throw an integration off by orders of magnitude. The source is TwoToThreeParticleScattering, the 2→3 peel used for any color-ordered set of three or more particles. It weights each point with the Byckling–Kajantie Jacobian 1/(8 sqrt(-G4)), forming the 4×4 Gram determinant G4 (and cos φ) directly from invariants of order s ~ 1e6 GeV². On the s23 range, G4 factorises:

-G4 = λ(s12, ma², t2)/16 · (s23 − s23_min)(s23_max − s23)

When a particle is soft the range is narrow. Terms of order s⁴ then have to cancel to a result smaller than their rounding error. Beyond ~1e-8 of relative softness the block returns garbage or NaN: 1e14 instead of O(1) at block level, in both the "peeled particle soft" and "previous particle soft" limits.

How often it hits

Massless 4- and 5-body phase-space volume at √s = 1 TeV, 8 batches × 2M points, estimate ÷ analytic (π/2)^(n−1) s^(n−2)/((n−1)!(n−2)!):

color order	origin/main, per batch	this PR, per batch	this PR, combined
single chain [0,2,3,4,5,1]	1039, 0.99, 1.00, 1.01, 1.00, 1.01, 0.99, 1e12	0.986 – 1.022	1.003
1+3 split [0,2,1,3,4,5]	1.02, 105, 0.99, 1.00, 1.01, 24, 1.00, 0.99	0.993 – 1.007	0.998
5-parton chain	1.9, 5455, 1.0, 1.7, 706, 194, 19, 0.99	0.980 – 1.094	1.009

The largest weight in the 4-parton chain falls from 2e18× to 1.6e4× the exact volume.

Fix

Both cos φ and the Jacobian are now computed from u = (s23 − s23_min)/(s23_max − s23_min), the position of s23 in the same range the invariant was sampled from:

• cos φ = 2u − 1
• 8 sqrt(−G4) = sqrt(λ) · (s23_max − s23_min) · |sin φ|, with |sin φ| = 2 sqrt(u(1−u))

Any rounding in the range cancels between the invariant's Jacobian and this one. The inverse reads |sin φ| off the momenta. bk_gram4, get_phi_from_s23 and lup_det4 have no remaining users and are removed.

At regular points old and new agree: median relative weight difference 1e-13, 99.9% of points within 1e-6. Forward/inverse round trips improve where the old determinant was ill-conditioned; for example, the worst |det·det_inv − 1| for a 5-body chain with t and W masses goes from 1.0 to 5e-4.

Tests

New madspace/tests/test_two_to_three_soft_limit.py:

• test_jacobian_with_soft_peeled_particle, test_jacobian_with_soft_third_particle: drive the block into both soft limits down to δ = 1e-14 and require the Jacobian to follow its smooth limit (within a factor of 2).
• test_color_ordered_massless_volume: the analytic 4-body volume for the single chain and the 1+3 split, averaged over 8 × 2M points, within 5%.

All 4 fail on origin/main and pass here. Full madspace suite: 1571 passed. The only local failures are 4 test_double_t tests needing scipy, the same on origin/main, plus 2 torch-only modules I skipped. Also built with ENABLE_SIMD=ON (arm64); the 2→3, t-channel, cut and new tests (288) pass there. The CUDA/HIP builds could not be tried locally; the kernel only uses functions already used on those backends (min, max, where, sqrt, acos, fabs).

Not addressed

• 1/sin φ edge peak. It sits at the edges of the s23 range and comes from sampling s23 rather than φ. The mean is correct but the variance is log-divergent; it is the remaining tail in the table above. Sampling u through e.g. u = sin²(πr/2) would cancel it, but that changes the s23 importance sampling, so it is left as a design decision.
• Soft-limit NaN points. A few points deep in the soft limit (≲ 1e-8) still come out as NaN. They get zero weight, and the region's contribution is negligible.
• Inverse kernel. The inverse is still unreliable at those points, and forward/inverse round trips with pT cuts fail for ~1% of points on both origin/main and this branch.
• Relation to madspace: hand every child's pair cuts to ColorOrderedMapping #173. Independent of that PR: different file, and this reproduces with no cuts on origin/main.

🤖 Generated with Claude Code

[ramonpeter]

@oliviermattelaer , thats nice. I will also check how to update the inverse then.
I need to fully check the math, but if this generally holds this might also be worth to add in the madspace paper.
This would then improve over Rikkerts implementation.

[ramonpeter]

@oliviermattelaer I pushed five commits on top of yours. Your two commits are unchanged. They address the "Not addressed" list (the 1/sin φ edge peak and the unreliable inverse) plus a few precision issues found while checking it. The full madspace suite passes on this branch and also merged with current main (incl. #173), and with ENABLE_SIMD=ON.

1. Arcsine map for s23 (removes the 1/sin φ tail)

s23 is linear in cos φ and the 2→3 measure is flat in φ. Any s23 density that is finite at the kinematic limits therefore leaves an integrable 1/|sin φ| in the weight, with a log-divergent variance. The random number of s23 is now first mapped flat in φ/2 and then handed to the existing s23 importance sampling, which is kept. The weight is bounded at the edges. With flat s23 sampling the block weight is exactly independent of r_s23. New option arcsine_s23 (default on) on TwoToThreeParticleScattering and ColorOrderedMapping.

Massless phase space at 1 TeV, 8 × 1M points, flat → arcsine:

	σ/μ	max w/μ
4-parton chain	11.3 → 4.5	8.1e3 → 340
1+3 split	7.9 → 4.6	5.4e3 → 350
5-parton chain	25.1 → 6.0	3.1e4 → 920
4-chain, pT/m_ij/ΔR cuts	8.8 → 3.4	9.6e3 → 190

Volumes are unchanged: the analytic value without cuts, and the fiducial volume with cuts.

2. Inverse / round trip

• With cuts: about 3% of points failed the forward/inverse round trip, carrying about 4% of the fiducial volume. Now none fail in 200k points. There were two causes:
  • rotate (frame of the 2→2 LAB block) switched between two conventions on the exact test q_t == 0. Forward and inverse round the beam-side momentum differently, so r_φ came back rotated, typically by exactly π/2. It is now the minimal rotation z → q, smooth across the beam axis and identical to the old one on it.
  • Block B applied its smin bound only when m2_2 < EPS. The inverse reads m² off the momentum (~1e-11 for a massless particle), so the bound switched off and the sampled s23 range changed. The massless test is now relative to the energy scale. The forward is unchanged.
• Near the edges of the s23 range: forward and inverse went through s23 itself, which only knows s23 − s23_min to ~1e-16 s. The block now carries u = (s23 − s23_min)/Δ and 1 − u in both directions. The new s23_arcsine_sample does the arcsine map and the importance sampling in one step via expm1/log1p, and the inverse reads u, 1 − u off the momenta. r_s23 now comes back to ~1e-17 absolute (was ~1e-8) and det·det_inv to ~1e-11 (was ~1e-5), also within 1e-12 of the edges.
• Inverse chain systems: these are now summed from the outgoing momenta, soft end first, instead of from beam-sized differences.

3. On-shell momenta, momentum still conserved by construction

Light-like momenta came out off shell by up to √|p²|/E ≈ 4e-5. This is why test_t_channel_masses[4 particles, massless-color_ordered] is flaky on main (fails for 2 of 40 seeds). Three causes:

• rotate_two_ref did a single Gram–Schmidt pass. With p3 close to pa the frame was not orthonormal, the rotation changed |p1|, and both daughters came out off shell. A second pass fixes it.
• Soft daughters formed as p12 − p1 (or boosted against a fast p12) lost precision. Now the softer daughter is boosted onto its mass shell in light-cone components and the harder one is p12 minus it, so p1 + p2 = p12 still holds by construction. This is used in the 2→3 and 2→2 LAB kernels.
• Each 2→3 block got its system as pa + pb − p3, built from beam-sized momenta. The kernels now take p12 directly (pb → p12, formed as before unless p12_condition is set), and the color-ordered walk hands each block the chain output of the previous one.

Max √|p²|/E is now ~4e-7 in the partonic frame (was 4e-5). Momentum conservation is unchanged, and the existing conservation tests pass with their original tolerances.

4. Boost to the hadronic frame

boost_beam lost a factor e^{2|y|} of precision for particles moving against the boost. It now takes the external masses and boosts those particles in light-cone components. It only does so for momenta that are already on shell to rounding; everything else keeps the linear boost. The resulting shifts are absorbed by the hardest outgoing particle, so conservation is exactly as before.

Lab-frame max |p² − m²|/E²:

• color-ordered: 2.8e-10 → 1e-12
• propagator: 5.9e-11 → 7.8e-12
• RAMBO: 5e-12 → 1.9e-13

What remains in the other modes comes from their partonic momenta.

Behaviour changes

• arcsine_s23 (on by default) and the new rotate change which event a given random number produces. The rotate change also affects the LAB 2→2 block in other modes. Trained flows or VEGAS grids need retraining; arcsine_s23=False restores the old s23 sampling.
• Kernel signatures in instruction_set.yaml changed: pb → p12 in the 2→3 kernels, the scattering kernels take u, u_c instead of s23, and boost_beam takes masses. New instructions: s23_arcsine, s23_arcsine_sample, s23_position (+ inverses).

6. Soft and collinear limits

This fixes the NaNs deep in the soft limit. I probed the block down to 1e-16 softness at the corners of the random numbers, and about a quarter of those points gave a non-finite inverse. There were two causes:

• rotate_two_ref floored |p3⊥| at 1e-12. A p3 below ~1e-12 GeV therefore gave a non-unit axis, which put the momenta off shell by ~1e4 GeV². The new light-cone boost had the same floor for a p12 almost at rest. Both frames are now normalised without a fixed scale.
• When p1 or p3 lies along pa, the s23 range has zero width and the weight is zero, but the inverse divided 0/0. Empty ranges were guarded with an absolute EPS, which is smaller than one ulp of s23 or t, so the guard did nothing. The sampling ranges now get a relative minimum width. The Jacobian still uses the true width, so the weight stays zero.

All probed points are now finite in both directions, and on shell to 1e-9 wherever the weight is nonzero. Round trips below ~1e-8 softness are still limited by the precision of the float momenta.

arcsine_s23 no longer supports a Breit–Wigner s23: it now raises an error, since no resonance is expected in this block. Without the arcsine map, Breit–Wigner sampling still works.

New tests: test_two_to_three_arcsine.py, test_color_ordered_precision.py and test_two_to_three_soft_limit.py, plus a strict inverse test in test_2to3_scattering.py. Numbers and details are in the commit messages.