Paper deep dive
Infrared Subtraction with Artificial Intelligence
Wenjie He, Xiaohui Liu, Yandong Liu, Zhan Wang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 9/30/2026, 7:42:14 AM
Summary
This paper presents an AI-developed framework for local infrared subtraction in perturbative QCD calculations, utilizing Large Language Models (LLMs) and neural networks. The method, based on Projection to Born (P2B) and Effective Field Theory (EFT) matching, separates integrable radiation terms from a finite 'Born contact' contribution. Two implementations are described: one using neural networks for phase space projection and EFT cumulant matching, and another using analytic constructions to calculate the contact directly from EFT singular distributions and finite radiation integrals. The approach eliminates the need for slicing parameters, reuses lower-order calculations, and demonstrates feasibility on modest hardware (Apple M1 CPU). Applications include NLO 3- and 4-jet production in electron-positron annihilation and NNLO dijet production.
Entities (11)
Relation Signals (8)
LLM → develops → Infrared Subtraction
confidence 97% · Under human physics guidance, an LLM develops two implementations.
AI-developed local infrared subtraction → uses → Projection to Born
confidence 95% · We present AI-developed local infrared subtraction, building on projection to Born and EFT matching.
AI-developed local infrared subtraction → uses → Effective Field Theory
confidence 95% · building on projection to Born and EFT matching.
Born Contact → determinedby → EFT singular distribution
confidence 94% · The contact is determined using the EFT singular distribution in a resolution observable such as N-jettiness
Neural Network → usedfor → phase space projection
confidence 92% · One uses a neural network for phase space projection
AI-developed local infrared subtraction → appliedto → NLO 3-jet production
confidence 90% · As a demonstration, we reconstruct the full NLO correction for massless 3- and 4-jet production
AI-developed local infrared subtraction → appliedto → NNLO dijet production
confidence 90% · The attempt to the NNLO dijet production is also made
AI-developed local infrared subtraction → comparedwith → EERAD3
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We present AI-developed local infrared subtraction, building on projection to Born and EFT matching. The framework separates an integrable radiation term from a finite contribution at Born kinematics, referred to as the Born contact. The contact is determined using the EFT singular distribution in a resolution observable such as N-jettiness $\tau_N$. Under human physics guidance, an LLM develops two implementations. One uses a neural network for phase space projection and fits the contact by matching to EFT cumulants. The other uses an analytic construction that keeps the Born momenta fixed while integrating over radiation. It combines the EFT $\delta(\tau_N)$ coefficient with finite 4-dimensional radiation integrals to calculate the contact term directly. This gives a local subtraction formula without a slicing parameter, while reusing existing lower-order radiation calculations and EFT singular predictions. As a demonstration, we reconstruct the full NLO correction for massless 3- and 4-jet production in electron-positron annihilation. The attempt to the NNLO dijet production is also made by recursively using the NLO P2B construction with the LLM designing machine-learning controls to reduce the variance of the contact integral. The tested predictions are in good agreement with EERAD3. The numerical calculation and projection-network training use a 2020 Apple M1 MacBook, without GPU acceleration, illustrating the feasibility of the construction with modest computing resources. The appendices develop an extension of the local subtraction to 3-jet NNLO, giving explicit radiation maps and a proposed contact formula. We also show how to integrate over NNLO radiation while keeping the Born momenta fixed, for any number of massless final-state jets. Our results demonstrate how AI can help higher-order calculations by constructing infrared subtraction and improving its numerical integration.
Tags
Links
- Source: https://arxiv.org/abs/2609.36007v1
- Canonical: https://arxiv.org/abs/2609.36007v1
Trouble viewing inline? Open PDF directly →
Full Text
167,853 characters extracted from source content.
Expand or collapse full text
Infrared Subtraction with Artificial Intelligence Wenjie He Email: 202531101011@mail.bnu.edu.cn Affiliation: School of Physics and Astronomy, Beijing Normal University, and Key Laboratory of Multiscale Spin Physics (Beijing Normal University), Ministry of Education, Beijing 100875, China Xiaohui Liu Email: xiliu@bnu.edu.cn Affiliation: School of Physics and Astronomy, Beijing Normal University, and Key Laboratory of Multiscale Spin Physics (Beijing Normal University), Ministry of Education, Beijing 100875, China Affiliation: Southern Center for Nuclear Science Theory (SCNT), Institute of Modern Physics, Chinese Academy of Science, Huizhou 516000, China Yandong Liu Email: ydliu@bnu.edu.cn Affiliation: School of Physics and Astronomy, Beijing Normal University, and Key Laboratory of Multiscale Spin Physics (Beijing Normal University), Ministry of Education, Beijing 100875, China Zhan Wang Email: zhanwang@mail.bnu.edu.cn Affiliation: School of Physics and Astronomy, Beijing Normal University, and Key Laboratory of Multiscale Spin Physics (Beijing Normal University), Ministry of Education, Beijing 100875, China September 27, 2026 Abstract We present AI-developed local infrared subtraction, building on projection to Born and effective-field-theory (EFT) matching. The framework separates an integrable radiation term from a finite contribution at Born kinematics, referred to as the Born contact. The contact is determined using the EFT singular distribution in a resolution observable such as N-jettiness τN _N. Under human physics guidance, a large language model (LLM) develops two implementations. One uses a neural network for phase space projection and fits the contact by matching to EFT cumulants. The other uses an analytic construction that keeps the Born momenta fixed while integrating over radiation. It combines the EFT δ(τN)δ( _N) coefficient with finite four-dimensional radiation integrals to calculate the contact term directly. This gives a local subtraction formula without a slicing parameter, while reusing existing lower-order radiation calculations and EFT singular predictions. As a demonstration, we reconstruct the full next-to-leading order (NLO) correction for massless three- and four-jet production in electron–positron annihilation. The attempt to the next-to-next-to-leading order (NNLO) di-jet production is also made by recursively using the NLO P2B construction with machine-learning controls to reduce the variance of the contact integral. The tested differential predictions are in good agreement with EERAD3. The numerical calculation and projection-network training use the CPU of a 2020 Apple M1 MacBook, without GPU acceleration, illustrating the feasibility of the construction with modest local computing resources. The appendices develop an extension of the local subtraction to three-jet NNLO, giving explicit radiation maps and a proposed contact formula. We also show how to integrate over NNLO radiation while keeping the Born momenta fixed, for any number of massless final-state jets. Our results show how AI could help higher-order calculations by constructing infrared subtraction and improving its numerical integration. I Introduction Real-radiation corrections are essential to precision QCD predictions. Additional partons generate resolved jets, redistribute recoil, and change event shapes and the acceptance of experimental cuts. Their contribution must therefore be included with its measurement dependence when extracting the strong coupling, constraining parton distributions, or comparing collider data with the Standard Model. The same radiation becomes singular when a parton is soft or two partons are collinear. Infrared safety ensures cancellation with virtual corrections in the physical cross section, but does not make the separate phase-space integrals finite [78, 80, 60]. Turning this cancellation into a differential numerical prediction remains a central problem in higher-order calculations. The treatment of unresolved radiation has developed alongside perturbative jet physics, from infrared-safe jet cross sections in the 1970s [93] and calculations of jet structure in electron–positron annihilation in the early 1980s [59] to general prescriptions based on universal soft and collinear factorization. The Frixione–Kunszt–Signer and Catani–Seymour methods established such prescriptions at next-to-leading order (NLO) in the 1990s [63, 36], opening the way to automated real-radiation subtraction [62]. At next-to-next-to-leading order (NNLO), simultaneous unresolved emissions introduce overlapping singularities. Sector decomposition provided a way to disentangle these singularities in real-radiation integrals [4]; antenna, sector-improved residue, nested soft-collinear, and local analytic sector subtraction developed complementary organizations of their cancellation [71, 51, 35, 85]. Phase-space slicing developed in parallel, with invariant-mass cutoffs and separate soft and collinear cutoffs used to isolate unresolved radiation [72, 73]. At NNLO, the qTq_T method exploited the singular small-recoil distribution for colour-singlet production [38]. Fully differential top-quark decay was subsequently calculated at NNLO using the invariant mass of the hadronic final state as a slicing variable [66]. A factorization theorem in soft-collinear effective theory (SCET) [12, 13, 15, 14, 23] determines the contribution below the mass cut, while the region above it is obtained from t→Wb+jett→ Wb+jet at NLO. This demonstrated how integrated singular information could replace the construction of counterterms for double-unresolved radiation. The definition and factorization of N-jettiness [94, 76] provided a resolution variable for processes with several hard directions; its use in NNLO calculations extended this strategy to jet production [27, 67]. Recently, recoil-free winner-take-all jet axes have also extended qTq_T-based slicing to jet final states [64]. These developments have made NLO calculations broadly automated [3, 26] and NNLO predictions available for a growing range of collider observables. The public NNLOJET framework covers jet production in lepton and hadron collisions [74], while recent antenna and local-sector developments pursue greater generality and automation [87, 24]. At N3LO, several complementary constructions have produced differential predictions. Projection to Born has been applied to jet production in deep-inelastic scattering [49, 68], Higgs decay to bottom quarks in combination with N-jettiness slicing [89], and gluon-fusion Higgs production [40]. Independently, qTq_T subtraction combined with NNLO antenna subtraction yielded the Higgs rapidity distribution, with an approximation for then-unknown third-order collinear terms [47], followed by neutral- and charged-current Drell–Yan distributions using the complete third-order small-qTq_T ingredients [44, 45]. In these calculations, the small-qTq_T factorization formula supplies the unresolved N3LO contribution, while antenna subtraction treats colour-singlet-plus-jet production at NNLO in the resolved region. Predictions with fiducial cuts have also been obtained through transverse-momentum subtraction and resummation in several implementations [25, 33, 43, 90]. More recently, antenna subtraction has been extended to the N3LO unresolved limits themselves and applied directly to e+e−→2e^+e^-→ 2 jets, yielding the two-jet rate and leading-jet energy distribution with singlet contributions neglected [46]. Fully differential Higgs-pair production has also reached N3LO in the heavy-top limit using qTq_T slicing [42]. Further progress includes the three-loop zero-jettiness soft function [39, 8, 9] and a transverse-momentum subtraction construction with a winner-take-all axis [79] for semi-inclusive deep-inelastic scattering at N3LO [53]. Despite this progress, extending this reach to more external legs and higher orders still requires substantial analytic work and numerical integration. Moreover, achieving the broad automation established at NLO remains a challenge at higher orders, where the subtraction terms, momentum maps, and integrated contributions often require dedicated analytic construction even when the amplitudes are available. To illustrate these difficulties, we compare subtraction and slicing at NLO for an N-parton Born process with no incoming coloured partons. We denote an m-parton configuration by Φm _m, its phase-space measure by dΦmd _m, and the value of an infrared-safe measurement on that configuration by JmJ_m. We write R for the real contribution and V for the ultraviolet-renormalized virtual contribution. With the Born term omitted, local subtraction expresses the NLO correction as [63, 36] δσNLO[J] δ _NLO[J] =∫dΦN+1[RJN+1−SJ~N] = _N+1 [RJ_N+1-S J_N ] (1) +∫dΦN[V+∫dΦradS]JN. + _N [V+ _radS ]J_N. Here the subtraction map m^:ΦN+1→ΦN m: _N+1→ _N takes a real-emission configuration to an N-parton Born configuration while preserving on-shell kinematics and total four-momentum. It approaches the corresponding Born limit when radiation becomes unresolved. We evaluate the projected measurement as J~N=JN(m^ΦN+1) J_N=J_N( m _N+1). The counterterm S reproduces R in each single-unresolved limit, making the first line finite. We integrate S over the unresolved radiation at fixed mapped Born momenta, including the phase-space Jacobian; sums over counterterms and maps are implicit. The resulting poles cancel those of V in dimensional regularization, so the second line is also finite in four dimensions. Slicing instead introduces a dimensionless resolution variable τ that vanishes on Born kinematics and resolves all relevant unresolved limits. With a cut tc>0t_c>0, the same correction takes the form δσNLO[J] δ _NLO[J] =∫dΦNΣN,sing(1)(ΦN,tc)JN = _N\, _N,sing^(1)( _N;t_c)J_N (2) +∫τ>tcdΦN+1RJN+1+Δpow[J;tc]. + _τ>t_cd _N+1\,RJ_N+1+ _pow[J;t_c]. The quantity ΣN,sing(1) _N,sing^(1) is the NLO singular cumulant below the cut, differential in Born kinematics. It combines virtual and unresolved real contributions and can be obtained from a factorization theorem [66, 67]. The second term is a finite real-emission integral above the cut. The remainder Δpow _pow contains the power-suppressed contributions, including the difference between the real and Born measurements below the cut, and vanishes as tc→0t_c→ 0. Dropping it defines the finite-cut slicing approximation. The two formulas show how the difficulties noted above arise in different parts of the calculation. Subtraction cancels the real singularity near each unresolved configuration and preserves the full measurement dependence, avoiding a neglected below-cut contribution. This can improve numerical convergence, but constructing the counterterms, their momentum maps, and their integrated forms becomes demanding as the multiplicity and perturbative order increase. Slicing based on factorization uses an integrated singular prediction calculated in effective field theory, reducing the need to construct local counterterms. The resolved contribution reuses a calculation with one additional jet at one lower perturbative order. However, its numerical cost grows when a small cut is needed to suppress power corrections. This is because lowering the cut exposes more singular regions of the resolved integral and produces large logarithmic contributions of opposite sign above and below the cut, whose cancellation requires increasingly precise numerical integration. In this work, we use a large language model (LLM) to advance multijet infrared subtraction, building on the established connection between projection to Born and EFT singular predictions [31, 57, 34]. Human guidance sets the framework and physical constraints, while the LLM devises momentum maps and contact formulas and implements the numerical calculations. In one implementation, it constructs the phase-space projection with a neural network and fits the Born contact through EFT matching. This tests whether the expressive power of neural networks can accommodate the kinematic and infrared constraints of subtraction. Furthermore, the LLM extends the same framework to local subtraction. It devises analytic maps that keep the Born momenta fixed during radiation integration and derives a formula that calculates the contact directly from the EFT δ(τ)δ(τ) coefficient and finite four-dimensional radiation integrals. The resulting construction evaluates both radiation and contact in four spacetime dimensions over the full phase space, without a small slicing parameter, while recycling existing lower-order radiation calculations and EFT singular predictions. We demonstrate both implementations for three-jet production at NLO and extend the local subtraction to four jets. At NNLO, we calculate the two-jet contact recursively from our NLO P2B construction. The LLM also designs a machine-learning strategy that reduces the variance of the contact integral. The appendices also develop the three-jet NNLO extension, with explicit radiation maps at fixed Born momenta and a proposed local contact formula. Human guidanceP2B + EFT matchingLLM−-SubtractionAnalytic / neural projectionδ singular inputNN-jettiness, qTq_T, …Complete fixed-orderresult Learned control functions Reduce integration variance Figure 1: Roles of human physics guidance, the LLM, machine learning and EFT input in the subtraction construction. Humans choose the P2B + EFT matching framework and specify the physical constraints. The LLM develops neural and analytic projections, contact formulas and their numerical implementation. EFT supplies the Born-dependent singular distribution. The contact is fitted through matching or calculated from finite radiation integrals and the EFT δ(τN)δ( _N) coefficient. In the LLM-developed local subtraction, learned control functions reduce the integration variance, with their known integrals added back to preserve the result. The radiation term and Born contact together give the fixed-order prediction. I P2B + EFT matching for multijet production Projection to Born (P2B) separates an integrable radiation contribution from a finite Born contribution by evaluating the higher-multiplicity input on the difference between real and projected measurements [31]. Its connection to differential N-jettiness subtraction provides an EFT determination of the Born contribution [67, 57]. We refer to this organization as P2B + EFT matching. The P2B-improved τ0 _0 construction in Ref. [34] treats colour-singlet production. It evaluates the difference between real and Born-projected measurements locally. The remaining contribution, differential in the projected Born variables, is the Born contact in our notation. This contact is still obtained by slicing, so a small cutoff is needed to suppress the residual power corrections. Building on these developments, we use AI to construct multijet infrared subtraction in two ways. First, we follow the established P2B + EFT matching strategy, using the expressive power of neural networks to construct the phase-space projection and determining the Born contact by direct fitting through EFT matching. This tests the ability of neural networks to handle the kinematic and infrared constraints of subtraction. Furthermore, we extend the framework to local subtraction. An analytic construction allows radiation integration at fixed Born momenta, so that the contact can be calculated directly from the EFT δ(τN)δ( _N) coefficient and finite four-dimensional radiation integrals. Both radiation and contact are then evaluated over the full phase space without a small slicing parameter. We demonstrate both implementations for three jets at NLO, and local subtraction for four jets at NLO and two jets at NNLO. Appendices B–D develop the extension to three-jet NNLO. Both constructions share the following P2B + EFT matching framework. For a process whose Born term starts at αsp _s^p, let ΣN(k) _N^(k) denote the coefficient at order αsp+k _s^p+k, with k≥1k≥ 1, and ℐN+1(k−1)I_N+1^(k-1) the coefficient at the same order in the calculation with one additional jet. Let M M denote the P2B projection. It leaves a Born configuration unchanged, M^ΦN=ΦN M _N= _N. We define JM^(Φ)=JN(M^Φ)J_ M( )=J_N( M ) and ΔJ=J−JM J=J-J_ M. Input multiplicities and branch or sector labels are displayed when needed. At any such order, the construction has the form ΣN(k)[J] _N^(k)[J] =ℐN+1(k−1)[ΔJ] =I_N+1^(k-1)[ J] (3) +∫dΦNcN(k)(ΦN)JN(ΦN). + _N\,c_N^(k)( _N)J_N( _N). This organization recycles a calculation with one additional jet at one lower perturbative order and an EFT singular prediction at the target order, preserving an advantage of slicing. The radiation calculation uses its existing internal subtraction, and M M must preserve its unresolved limits so that ΔJ J admits local cancellation and an integrable radiation term. The finite density cN(k)c_N^(k) depends on the projection prescription and contains the remaining Born contribution in the same normalization. To determine this contact, choose a dimensionless resolution observable τN _N that vanishes on Born configurations, and weight functions wa(ΦN)w_a( _N) on Born phase space. For t>0t>0, define Wa,t(Φ)=wa(M^Φ)Θ[t−τN(Φ)]W_a,t( )=w_a( M ) [t- _N( )] and its projected measurement Wa,t,M^(Φ)=wa(M^Φ)W_a,t, M( )=w_a( M ). The matching relation is ∫dΦNcN(k)wa=limt→0 _N\,c_N^(k)w_a= _t→ 0 \ ΣN,sing(k)[Wa,t] _N,sing^(k)[W_a,t] (4) −ℐN+1(k−1)[Wa,t−Wa,t,M^]. -I_N+1^(k-1)[W_a,t-W_a,t, M] \. Here ΣN,sing(k) _N,sing^(k) is the EFT singular cumulant at the same order, including its Born dependence and the finite coefficient of δ(τN)δ( _N) in the singular distribution. The subtracted term is the radiation contribution in Eq. (3), evaluated with J=Wa,tJ=W_a,t. Equation (4) thus determines a weighted moment of the contact for each choice of waw_a. Weights localized in different Born regions, or chosen from a function basis, probe its Born dependence. A neural-network or other parametrization can therefore be fitted by requiring its weighted integrals to reproduce these moments, as illustrated by the three-jet NLO fit in Sec. VII. For the analytic implementation, we parametrize the higher-multiplicity phase space by Born coordinates and radiation variables ξ. The map G(ΦN,ξ)G( _N,ξ) generates radiation that projects back to the prescribed Born configuration, M^G(ΦN,ξ)=ΦN MG( _N,ξ)= _N. Holding the Born momenta fixed, we integrate only over ξ, so the matching relation determines the contact at each Born point. Appendix A gives the analytic NLO projection and its inverse for any number of massless final-state jets. Appendix D derives a direct contact calculation that avoids this small-cut limit. At NLO, we add the finite part of the radiation integral at fixed Born momenta to the EFT δ(τN)δ( _N) coefficient. Denote this coefficient by H1(b)H_1(b) at a Born point b, and let 1(ζ,b)Z_1(ζ;b) be the tree-level radiation integral weighted by τNζ _N^ζ, with the same projection prescription. The contact is then cN(1)(b)=H1(b)+1(b),1(b)=[ζ0]1(ζ,b).c_N^(1)(b)=H_1(b)+A_1(b), _1(b)=[ζ^0]Z_1(ζ;b). (5) Here [ζ0][ζ^0] selects the constant term in the Laurent expansion of 1Z_1 around ζ=0ζ=0; this constant is 1A_1. To calculate it, we combine radiation contributions evaluated at coordinates related by fixed rescalings, so that their common soft and collinear terms cancel before integration. We also include the finite contributions from the changes in integration range. The result is a sum of convergent integrals in four spacetime dimensions, retaining the nonsingular radiation over the full phase space. No small numerical value of ζ or slicing cutoff is required. Together with Eq. (3), this gives a local subtraction formula at NLO while reusing existing radiation and EFT results. Once the projection and contact are fixed, Eq. (3) predicts different observables without further matching. The numerical examples compare the two implementations for three jets at NLO and demonstrate local subtraction for four jets. They also demonstrate recursive local contact extraction at NNLO in the simpler two-jet case, where the observables require only one contact coefficient. Appendices B–D propose the extension of the local subtraction to three-jet NNLO, combining radiation integration at fixed Born momenta with a finite-difference representation of the contact. I The role of artificial intelligence Machine learning supports event reconstruction, simulation, and analysis in high-energy physics [2], while LLM agents have begun to develop analyses and coordinate simulation workflows [52, 88, 65]. In perturbative QCD, applications include machine learning for loop-integral differential equations [32, 82], AI-assisted integration-by-parts reduction [95, 92], and factorization-aware matrix-element emulation [86]. LLM-driven bootstrapping has also reached the nine-loop six-gluon maximally helicity-violating amplitude in planar =4N=4 super-Yang–Mills theory [48]. Work related to infrared subtraction includes GAN generation of subtracted event samples using prescribed collinear counterterms [30] and neural importance sampling within NNLO subtraction [75]. These developments motivate using AI to devise the subtraction itself. The projection must satisfy kinematic constraints and reproduce the soft and collinear limits, but its form away from these limits is not uniquely fixed. Recent LLM-guided searches for mathematical constructions illustrate the ability to devise algorithms under explicit constraints [58, 1, 91]. Here the LLM devises momentum maps that preserve on-shell kinematics, four-momentum conservation, particle symmetries, and unresolved limits. The analytic construction gives explicit radiation maps at fixed Born momenta for massless final-state jets of any multiplicity, while the neural network parametrizes the projection away from unresolved limits. The LLM also implements both constructions, including training for the neural network. Human physics guidance specifies the constraints, while the LLM develops and refines the construction. Within this construction, neural networks provide the capacity to represent functions of many variables. Universal approximation on compact domains [50] motivates using them to parametrize the projection away from soft and collinear limits. We test whether this flexibility can be combined with exact kinematic and infrared constraints to produce an integrable radiation term. A neural parametrization could also represent the signed Born contact; the numerical fit presented here uses a spline basis. Machine learning has a further role after local subtraction has made the radiation integrand integrable. A function learned from integration samples can be subtracted from the exact weight and its known integral added back, reducing fluctuations without changing the prediction. We realize this strategy with fitted Gaussian and polynomial control functions in the two-jet NNLO calculation of Sec. VII.3. As illustrated in Fig. 1, humans choose the framework and specify the physical constraints. The LLM devises the projections and contact formulas and implements the numerical calculation, including the learned controls for integration. For the local NLO calculation, it also assembles the EFT δ(τN)δ( _N) coefficient from established hard, jet and soft contributions. We used GPT-5.6 and GPT-6 for this construction and implementation. The neural implementation tests a flexible representation of the projection. The LLM derives radiation maps at fixed Born momenta and finite contact integrals, yielding the local subtraction detailed in Appendices A and D. The two-jet NNLO example combines recursive contact extraction with learned variance reduction. Appendices B–D develop the extension to three-jet NNLO. Table 1: Physical inputs and construction restrictions for AI-developed P2B subtraction. The exclusions concern the new subtraction being constructed. The lower-order radiation input may use its existing subtraction, as indicated by the dash. Ingredient Input or permitted use Excluded from the new construction Projection Born information, parton kinematics and identities, universal infrared limits and geometric constraints. Reuse of local counterterms or momentum maps from existing subtraction schemes. Radiation The complete calculation with one additional jet at one lower perturbative order, including its internal subtraction. — Contact EFT singular predictions at fixed Born kinematics, including the coefficient of δ(τN)δ( _N), after fixing the projection. Separate virtual or integrated-counterterm contributions as contact-fitting targets. Prediction A common fixed projection and contact; independent reference comparisons. Observable-dependent changes to the fixed projection or contact. IV NLO workflow In both NLO implementations, we first construct a projection satisfying the physical constraints. We then fix the projection, determine its contact from EFT input and predict several observables using the same map and contact. To test the LLM’s ability to carry out this construction, we required it to develop the maps and contact calculation from the given constraints, without importing local counterterms or momentum maps from existing subtraction schemes. Table 1 summarizes the inputs and restrictions, and Fig. 2 shows the sequence. Human guidance P2B + EFT matching Physics constraints LLM designs map Analytic / neural Check constraints Train if neural Projection M^4 M_4 refine parametrization P2B integral FM^[J]F_ M[J] Determine Born contact NLO observable distributions Comparison with EERAD3 EFT singular predictionat fixed Born kinematics Figure 2: Common workflow for the two NLO implementations. Humans specify the framework and physical constraints. The LLM develops a projection and checks its kinematic and infrared properties, including training for the neural option. With the projection fixed, the contact is fitted through EFT matching or calculated directly using finite radiation integrals and the EFT δ(τ3)δ( _3) coefficient. Each implementation uses its own common projection and contact for all final distributions. (a) Projection module Four momenta, labels Ei/QE_i/Q, sij/Q2s_ij/Q^2 Branch selection 10 kinematic inputs + 6 label entries Mapper MLP16→32→32→1116→ 32→ 32→ 11 Shape + orientation Constrained decoder One state Φ~3 _3 Analytic limiting behaviour; symmetry constraints; geometric training criterion(b) Real-emission integral Integration sample Φ4 _4, exact R4R_4 Projection (Φ4,Φ~3)( _4, _3) Measure both states with specified support Integrate weighted measurement difference (c) Contact module EFT singular cumulants weighted by wa(b)w_a(b) EFT minus P2B t→0t→ 0 limit Contact as a function of Born coordinates Born contact cM^(b)c_ M(b) with uncertainty The same projection and contact enter all observable distributions(d) Learned integration controls for the local NNLO contact Radiation integrand for the contact Learn a control from separate samples Integrate exact integrand minus control Add the known control integral Control held fixed during independent integration; unchanged physical prediction Figure 3: Modules for projection, radiation integration and contact determination. Panels (a)–(c) describe the neural implementation at NLO. The projection network has two hidden layers with sigmoid linear unit (SiLU) activations. Its 11 outputs comprise three shape parameters and eight coefficients for two orientation vectors. The real-emission integral uses the exact matrix element. Contact fitting uses EFT singular cumulants weighted by functions of the Born variables, in the limit t→0t→ 0. Here Φ~3=M^4Φ4 _3= M_4 _4. Panel (d) shows the additional use of machine learning in local two-jet NNLO contact integration. A control function is trained on separate samples, then held fixed while the exact residual is integrated and the known control integral restored. Panels (a)–(c) of Fig. 3 give the three modules of the neural implementation. The projection module selects a branch labelled by a soft parton or collinear pair. Here Q is the centre-of-mass energy, pip_i is a parton four-momentum, Ei=pi0E_i=p_i^0, and sij=2pi⋅pjs_ij=2p_i· p_j, with p2=(p0)2−||2p^2=(p^0)^2-|p|^2. Four energy fractions and six pair invariants, together with six learned components encoding this label, enter a multilayer perceptron (MLP). Its outputs determine a single projected state through a constrained parametrization. The network is trained on kinematic information; pointwise real weights and reference NLO distributions do not enter its training or selection. The second module integrates the exact real coefficient against the difference between real and projected measurements. The contact module fits a function of Born coordinates to the t→0t→ 0 difference between EFT singular cumulants and the P2B integral. Panel (d) illustrates the learned variance reduction used for the local NNLO contact below. In the analytic implementation, the map allows radiation to be generated while the Born momenta remain fixed. The radiation term uses the real and projected measurements as before, while the contact is calculated from the finite radiation contribution 1A_1 and the EFT δ(τ3)δ( _3) coefficient H1H_1. The next section describes the two constructions. All numerical calculations presented here, including projection-network training and contact extraction, use the CPU of a 2020 Apple M1 MacBook, without GPU acceleration. V NLO construction and contact extraction We consider a pure virtual-photon current, nf=5n_f=5 massless quark flavours, and μR=Q _R=Q with the MS¯ MS coupling. We work in the centre-of-mass frame with flavour-blind, rotationally invariant measurements. With as=αs(Q)/(2π)a_s= _s(Q)/(2π) and the photon-mediated two-parton Born cross section σ0 _0, the three-jet expansion is σ[J]σ0=asA[J]+as2B[J]+(as3). σ[J] _0=a_sA[J]+a_s^2B[J]+O(a_s^3). (6) The displayed three-jet NLO results give the coefficient B in this expansion. All coefficient densities below have the corresponding powers of asa_s and the common factor σ0 _0 removed and use Lorentz-invariant phase-space measures. For the rotationally invariant measurements used here, we average the kernels over the overall event orientation. The real channels are qq¯ggq qg and four-quark production, with the flavour sums, interference terms, and correlated azimuthal samples supplied by EERAD3 [70]. These samples pair configurations differing by a π/2π/2 azimuthal rotation of a daughter pair in its rest frame. V.1 Neural projection and contact fitting Let Γm0 _m^0 select JADE jets [10] with E-scheme recombination and strict y23>0.02y_23>0.02; no y34y_34 cut is applied. Here y23y_23 and y34y_34 denote the 3→23→ 2 and 4→34→ 3 jet transition values. The neural projection M^4:Φ4→Φ3 M_4: _4→ _3 is evaluated only for Γ40=1 _4^0=1. Its projected counterevent contributes only when Γ30(M^4Φ4)=1 _3^0( M_4 _4)=1. For a measurement J supported within this three-jet region, define FM^[J]=∫Γ40=1dΦ4R4(Φ4)[J4(Φ4) F_ M[J]= _ _4^0=1d _4\,R_4( _4) [J_4( _4) (7) −Γ30(M^4Φ4)J3(M^4Φ4)]. - _3^0( M_4 _4)J_3( M_4 _4) ]. Here R4dΦ4R_4d _4 includes the normalization of the real contribution to B. All additional physical cuts are evaluated on each state separately. In a single-unresolved limit with a resolved three-parton limiting configuration, the map approaches the corresponding energy flow. The real and projected weights therefore cancel their leading singular behaviour before integration, giving an integrable FM^[J]F_ M[J] for a fixed infrared-safe measurement without a slicing boundary. At finite resolution, the map changes the finite subtraction prescription. In this normalization, Eq. (3) at k=1k=1 becomes B[J]=FM^[J]+∫dΦ3(b)cM^(b)J3(b).B[J]=F_ M[J]+ _3(b)\,c_ M(b)J_3(b). (8) Here b denotes the Born momenta, ordered by decreasing energy, and cM^c_ M is a density relative to dΦ3(b)d _3(b), with the parton-label assignments summed. The contact cM^c_ M contains the finite combination of virtual and projected real contributions. It depends on the map and support prescription, while their sum in Eq. (8) does not. The fit parametrizes its ratio to the Born density, as specified in Sec. VII.1. For the projection entering Eq. (8), we parametrize the shape and orientation of a massless three-parton state. An analytic reference is obtained by removing the selected soft parton or combining the selected collinear pair. We subtract the mean of the three spatial momenta, assign massless energies, and rescale to total energy Q. This recovers the Born state in the corresponding unresolved limit; an equal-energy configuration is used if the reference degenerates away from that limit. The network deforms this reference through Li L_i =lnuianchor+2d1+dtanhℓi, = u_i^anchor+ 2d1+d _i, (9) ui u_i =eLi∑jeLj,Ei=Q2(1−ui). = e^L_i _je^L_j, E_i= Q2(1-u_i). Here uianchor=1−2Eianchor/Qu_i^anchor=1-2E_i^anchor/Q specifies the reference shape, and ℓi _i are network outputs. We use d=2Es/Qd=2E_s/Q for a soft parton of energy EsE_s and d=(1−cosθij)/2d=(1- _ij)/2 for a collinear pair of opening angle θij _ij, so that the deformation vanishes in the corresponding unresolved limit. The orientation is constructed covariantly, with the required particle-exchange symmetries, and all three momenta remain on shell and conserve total four-momentum. The geometric training criterion contains neither the four-parton matrix element nor the NLO cumulants. Since Eq. (7) is integrated over the original four-parton phase space, no Jacobian for the projection is required. With the projection fixed, we determine the contact through its moments in Born coordinates. We use geometric N-jettiness [94, 76] with reference directions chosen by global minimization, τN(Φm) _N( _m) =1Qmin^a∑k=1mmina=1,…,Nna⋅pk, = 1Q _\ n_a\ _k=1^m _a=1,…,Nn_a· p_k, (10) τ3(Φ4) _3( _4) =mini<jEi+Ej−|i+j|Q. = _i<j E_i+E_j-|p_i+p_j|Q. Here na=(1,^a)n_a=(1, n_a) with |^a|=1| n_a|=1, so na⋅pk=Ek−^a⋅kn_a· p_k=E_k- n_a·p_k. For four massless partons, the minimizing partition consists of one pair and two singletons. Their axes align with the respective momentum sums, yielding the second line. On three massless partons, τ3=0 _3=0. For N=3N=3, the measurement Wa,tW_a,t of Eq. (4) has Born value wa(b)w_a(b). Its contact moment in the present normalization is CM^,a=∫dΦ3(b)cM^(b)wa(b).C_ M,a= _3(b)\,c_ M(b)w_a(b). (11) Since three-jettiness vanishes on every projected Born state, the contact contributes CM^,aδ(τ3)C_ M,aδ( _3) to the Born-weighted distribution. The matching relation in Eq. (4) reads CM^,a=limt→0Bsing[Wa,t]−FM^[Wa,t].C_ M,a= _t→ 0 \B_sing[W_a,t]-F_ M[W_a,t] \. (12) Here BsingB_sing is the leading-power EFT prediction for the same measurement Wa,tW_a,t and phase-space cuts, including its δ(τ3)δ( _3) term and plus distributions [67]. The logarithms of t cancel between Bsing[Wa,t]B_sing[W_a,t] and FM^[Wa,t]F_ M[W_a,t], while the nonsingular cumulant vanishes as t→0t→ 0. The remaining constant determines the contact moment. This requires the finite EFT δ(τ3)δ( _3) coefficient, including the hard contribution, as well as the logarithmic terms. We determine the Born-dependent contact by fitting its weighted integrals to these moments. Section VII.1 gives the parametrization and numerical results. V.2 Analytic projection and local contact The second implementation calculates the contact directly at each Born point. The analytic map in Appendix A merges a selected pair and rescales the remaining momenta, preserving the mass shell, total four-momentum and single-unresolved limits. Its inverse generates the four-parton state from the prescribed Born momenta and three radiation variables. This defines the common projection and phase-space measure needed for the radiation term and the contact. At fixed Born kinematics, we evaluate 1A_1 using the finite four-dimensional integrals in Eq. (55) and add the EFT δ(τ3)δ( _3) coefficient H1H_1 to obtain cM^=H1+1c_ M=H_1+A_1. Appendix D derives these weights and the compensations that restore the full integration domain. The resulting local subtraction formula evaluates both radiation and contact without a small slicing parameter and includes the exact nonsingular contributions. Section VII.2 compares its three-jet predictions with the same distributions used for the neural implementation and presents the extension to four jets. The extension of the local subtraction to three-jet NNLO is proposed in Appendices B–D. VI Two-jet NNLO example The two-jet case provides a simple test of the P2B organization at NNLO. The radiation term in Eq. (3) is obtained from the NLO correction to three-jet production. We reuse this complete NLO calculation, including its internal subtraction, and apply P2B to remove the remaining singularity associated with the approach to two-parton Born kinematics [31, 49]. The radiation term preserves the full infrared cancellation within the NLO input, including its integrated contributions. The measurement ΔJ J is evaluated on every real configuration and every subtraction configuration. In a single-unresolved limit it must approach the corresponding lower-multiplicity measurement difference; on reaching the two-parton Born configuration it vanishes. These limits must be preserved, with sufficient suppression for the remaining integration to be finite. No extra-jet generation cut may remove a nonzero contribution to ΔJ J. A completed NLO subtraction, including one constructed with AI, can therefore be reused at the measurement level. We use the same massless, pure-photon theory as at NLO. For measurements of thrust [61], the parton energy–energy correlation (EEC) [11], and the two highest jet energies, the Born value jBj_B is independent of the orientation of the back-to-back pair. The Born contact therefore reduces to one scalar. Writing the coefficient of σ0as2 _0a_s^2 as ℬ2B_2, Eq. (3) becomes ℬ2[J]=ℛ3[J−jB]+c2jB,B_2[J]=R_3[J-j_B]+c_2j_B, (13) where ℛ3R_3 is the three-jet NLO correction in this normalization. For these rotationally invariant measurements the three-to-two projection is trivial, since only the Born value jBj_B is needed. The LLM therefore implemented this step without a neural network. For this two-direction Born configuration, geometric two-jettiness with back-to-back thrust axes is τ2=1−T _2=1-T. A Born-dependent or axis-sensitive measurement requires a nontrivial projection. The scalar contact c2c_2 contains the finite combination of double-virtual and Born-projected radiation contributions. We calculate it directly as c2=H2+2c_2=H_2+A_2, recursively using the NLO P2B construction as described in Appendix D.4. Here H2H_2 is the coefficient of δ(τ2)δ( _2) in the NNLO EFT thrust distribution, with the hard, jet and soft contributions combined [18, 77], and 2A_2 is the finite radiation contribution. Both the radiation term and the contact are evaluated by local subtraction without a slicing cut. Section VII.3 presents the numerical results and the learned controls used to stabilize the contact integral. Appendices B and D develop the extension to a Born-dependent three-jet NNLO contact. VII Numerical results VII.1 Three-jet NLO with a fitted contact Figure 4: Complete three-jet NLO coefficient distributions with strict JADE/E y23>0.04y_23>0.04 and no y34y_34 cut. The projection and the 127-parameter contact are common to all panels. Blue bands combine P2B Monte Carlo and contact bootstrap errors; black points show EERAD3 with Monte Carlo errors. In all distribution plots, O denotes the observable on the horizontal axis. Lower panels show Δ/σΔ / _ , where Δ is the difference from the reference and σΔ _ its combined statistical error. The same event sample contributes to different bins, and the fitted contact introduces additional correlations. The matching prescription in Eq. (12) uses EFT singular cumulants in the limit t→0t→ 0. For simplicity in this numerical test, we use full reference cumulants from EERAD3 at t=0.004,0.025,0.09t=0.004,0.025,0.09 to constrain the contact fit. At these finite thresholds, the exact identity CM^,a=B[Wa,t]−FM^[Wa,t]C_ M,a=B[W_a,t]-F_ M[W_a,t] includes the nonsingular contribution and requires no small-t approximation. These Monte Carlo cumulants carry statistical errors and require weighted integrations to resolve the Born dependence, making contact extraction more demanding than with an analytic EFT cumulant evaluated at chosen Born points [76, 67]. For the fit, we use the ordered energy fractions xi=2Ei/Qx_i=2E_i/Q, with x1≥x2≥x3x_1≥ x_2≥ x_3 and x3=2−x1−x2x_3=2-x_1-x_2. After integrating over the overall orientation, dΦ3=Q2dx1dx2/(128π3)d _3=Q^2dx_1dx_2/(128π^3). Let B3(b)B_3(b) be the three-parton Born coefficient and B¯3(b)=∑π∈S3B3(bπ) B_3(b)= _π∈ S_3B_3(b_π) its sum over the six assignments of q,q¯,gq, q,g to these momenta. We expand the dimensionless ratio cM^(b)/B¯3(b)c_ M(b)/ B_3(b) in bicubic splines with 127 independent coefficients and C2C^2 continuity at x1=x2=2/3x_1=x_2=2/3. The fit integrals therefore use the Born weight dΦ3B¯3d _3\, B_3. The coefficients may have either sign, as required for the finite virtual contribution. The statistical uncertainty is estimated from 128 bootstrap replicas of the cumulant fit. The resulting contact is common to all distributions in Fig. 4. We reconstruct the full three-jet NLO coefficient B by combining the subtracted real-emission integral with the fitted Born contact, including the finite virtual contribution. The contact is determined in the region y23>0.02y_23>0.02. Figure 4 presents the resulting distributions with the tighter cut y23>0.04y_23>0.04, using the same projection and contact throughout. Total broadening BTB_T [37] is measured on the original partons; the ordered jet energies follow exclusive three-jet recombination, with total visible energy Evis=QE_vis=Q. For this comparison, the EERAD3 reference is evaluated with its own subtraction prescription, applying the cuts separately to real events and counterevents. The comparison therefore probes both the observable dependence and the change in the three-jet cut. Figure 4 shows good agreement with EERAD3, reproducing the full NLO coefficient with a single fitted contact and no observable-dependent refit. For broadening and leading energy, the integrated NLO coefficients are 286.43±1.00286.43± 1.00 and 287.75±1.25287.75± 1.25, respectively, compared with 287.09±0.28287.09± 0.28 and 287.72±1.21287.72± 1.21 from EERAD3. For second-leading energy we obtain 288.33±0.92288.33± 0.92, compared with 285.50±0.76285.50± 0.76, a difference of 2.36 combined standard deviations. These normalizations agree within two combined standard errors after accounting for the common contact. The EERAD3 second-leading-jet reference lies slightly below its other determinations, also by less than two combined standard errors, consistent with fluctuations between independent reference runs. The largest binwise differences are 1.71, 3.04 and 1.91 combined standard deviations for broadening, leading and second-leading energy, respectively. For this numerical study, we deliberately restricted the computational budget to a single M1 laptop with limited running time, keeping the integration statistics modest and the contact parametrization to 127 coefficients. At fixed sampling, narrow bins have fewer contributing events, while away from unresolved limits real and projected configurations can populate neighbouring bins. Such fluctuations are common to histogrammed subtraction calculations. EERAD3 reduces the sampling cost with adaptive VEGAS grids [81] and phase-space parametrizations adapted to unresolved limits [70]. Our calculation also uses adaptive integration, but neither the geometric projection training nor the sampling was optimized for individual distributions. Sampling adapted to the bins and a projection optimized for variance can improve precision per event, while more precise cumulants weighted by functions of the Born variables can constrain a richer contact representation. These improvements offer a direct route to higher precision for the reconstructed NLO distributions. VII.2 Three- and four-jet NLO with a local contact We next calculate the same three-jet NLO distributions using the analytic local subtraction of Sec. V.2. The map in Appendix A generates radiation at fixed Born momenta and defines a common projection and measure for the radiation term and the contact. At each Born point, Eq. (57) gives cM^=H1+1c_ M=H_1+A_1. The LLM obtains H1H_1 by expanding the three-jettiness factorization formula [76] to one loop and collecting the coefficient of δ(τ3)δ( _3). It derives the hard contribution from the finite three-parton virtual interference [71] supplied by EERAD3, converting it to the MS¯ MS hard-function convention. It then combines this with the one-loop quark and gluon jet functions [17, 16] and the three-jettiness soft function, including its finite geometry-dependent terms [76]. All terms use μR=Q _R=Q, the geometric jettiness normalization with reference scales Qi=2EiQ_i=2E_i, and the dimensionless τ3 _3 with the plus-distribution convention of Eq. (56). The LLM also includes the finite δ(τ3)δ( _3) terms generated when expressing the jet and soft distributions in this common variable. The finite radiation contribution 1A_1 is calculated directly with the weights of Eq. (55), including all finite-interval compensations. Adding H1H_1 completes the contact at each Born point. Both the radiation term and contact are integrated in four dimensions, including all nonsingular contributions. This gives a local subtraction calculation without a small slicing parameter; its inputs are the tree matrix elements and the EFT singular prediction. No contact fit or reference cumulant enters this prediction. Figure 5: Three-jet NLO coefficient distributions from the analytic local subtraction with a directly calculated contact, using the same observables and physical cuts as Fig. 4. Green, orange and purple curves show the P2B radiation term FM^F_ M, the Born-weighted EFT δ(τ3)δ( _3) coefficient H1H_1, and the Born-weighted finite radiation integral 1A_1, respectively. Their sum is the blue result; its band includes the correlation between radiation and contact. Black points reproduce the EERAD3 reference of Fig. 4. Lower panels compare the total with the reference in combined standard errors, with grey bands at ±2± 2. The calculation uses the maps of Appendix A and the local contact weights of Appendix D, over the full radiation phase space. Figure 5 separates the P2B radiation contribution, labelled FM^F_ M, from the Born-weighted H1H_1 and 1A_1 terms. Their sum agrees well with EERAD3. The calculation uses 3.36×1083.36× 10^8 samples and gives a common integrated coefficient B=286.680±0.277B=286.680± 0.277, with a relative statistical error of 0.10%0.10\%. The largest binwise differences from EERAD3 are 1.561.56, 2.962.96 and 1.901.90 combined standard errors. These results demonstrate complete multijet NLO predictions from the AI-constructed analytic projection and direct local contact calculation. Figure 7 compares the contact 1+H1A_1+H_1 obtained with finite rescalings ρ=0.1ρ=0.1 and 0.010.01 and by direct evaluation of the P2B integral 1A_1 in Eq. (52). The same EFT δ(τ3)δ( _3) coefficient H1H_1 is added in all three cases. Here ρ∈(0,1)ρ∈(0,1) is a freely chosen rescaling parameter, not a cutoff. The exact compensations in Appendix D restore the full phase space for any such choice. The two rescalings use independent samples, while Eq. (52) is evaluated on those same samples. The differences shown below the distributions account for their correlations. All three determinations agree, with the largest binwise difference from the direct integral below 2.42.4 standard errors. This supports both the equivalence of the two contact formulas and the independence of the extracted contact from the free parameter ρ. We further show four-jet NLO predictions in Fig. 7, using the analytic projection and local contact with four-parton Born kinematics. The contact combines the EFT δ(τ4)δ( _4) coefficient, calculated directly by the LLM, with finite radiation integrals, without a small slicing parameter. All distributions agree well with EERAD3.11 1 The LLM identified 7**s12 as a typo in the four-parton one-loop nfn_f contribution of EERAD3 2.0.0. Its symbolic checks show that replacing it by 7*s12 restores the expected scaling with the invariants and the symmetry under simultaneous quark–antiquark and gluon exchange. The correction also makes the leading qgqg and q¯g qg collinear residues vanish, as required for this nfn_f term. We use this LLM-suggested matrix-element modification here and for K4K_4 in Fig. 11. Figure 6: Comparison of three contact determinations for the analytic projection. Each curve shows the Born integral in Eq. (8), denoted BcontactB_contact, for the analytic projection and the cuts of Fig. 4. The black histogram uses the finite radiation integral 1A_1 in Eq. (52) plus the EFT δ(τ3)δ( _3) coefficient H1H_1. Blue circles and orange squares use Eq. (55) for 1A_1, with ρ=0.1ρ=0.1 and 0.010.01. Error bars and the black band denote one standard error; markers are displaced within each bin for visibility. Lower panels show the differences from the black result in standard errors, including shared-sample correlations; grey bands indicate ±2± 2. As defined in Appendix D, ρ is a free rescaling parameter, not a cutoff. Figure 7: Four-jet NLO coefficient distributions from the analytic local subtraction with a directly calculated contact, with strict JADE/E y34>0.04y_34>0.04. Here 4C_4 is the coefficient of as3a_s^3 in σ[J]/σ0σ[J]/ _0. Total broadening is measured on the partons, and the jet energies follow exclusive four-jet recombination with Evis=QE_vis=Q. Blue curves and bands show the P2B + EFT prediction and one standard error; black points show the EERAD3 antenna reference with Monte Carlo errors. Lower panels show the differences in standard errors, including shared-sample correlations, with grey bands at ±2± 2. VII.3 Two-jet NNLO We now turn to the two-jet NNLO results, comparing a directly calculated contact with one determined by cumulant matching. The local construction recursively uses our NLO P2B contact, as described in Appendix D.4, and requires no slicing cut. As in the NLO examples, the LLM constructs the EFT input H2H_2 by combining the known hard and quark-jet contributions with the analytic two-loop hemisphere-soft constant [18, 77]. It collects the coefficient of δ(τ2)δ( _2) in the NNLO singular thrust distribution, using the plus-distribution convention of Appendix D. Figure 8: Two-jet NNLO coefficient distributions ℬ2B_2. Red dashed curves use local subtraction with the directly calculated contact. The red bands combine P2B and contact errors, retaining the common contact correlation between Born-supported bins. Green solid curves show the cumulant-matched calculation, and black points the EERAD3 reference. Lower panels compare the red prediction with the reference in combined statistical errors; grey bands indicate ±2± 2. Shaded upper-panel bins contain Born support. There are no event-selection cuts; EEC includes all ordered parton pairs and self correlations, and JADE/E clustering at ycut=0.04y_ cut=0.04 includes two-, three-, and four-jet events. Local cancellation makes the contact integral finite, but signed contributions can still produce substantial Monte Carlo fluctuations. We reduce these fluctuations with control functions learned from previous integration samples. Write the radiation contribution to the contact as 2=∫dμ2A_2= μ\,W_2, where dμdμ includes all radiation variables and channel sums. Let q be a normalized sampling density on this domain, so that 2/qW_2/q is the Monte Carlo weight. For a learned function hθh_θ with fitted parameters θ and known integral IθI_θ, we use the identity 2=Iθ+⟨2−hθq⟩q,Iθ=∫dμhθ.A_2=I_θ+ W_2-h_θq _q, I_θ= μ\,h_θ. (14) The brackets denote an average over qdμq\,dμ. Gaussian and polynomial control functions are fitted to the signed integrand 2W_2 using separate samples, then held fixed during integration. Their known integrals are added back, so the learned approximation affects the variance without changing the contact. No inclusive cross section is used to train or normalize the prediction. The contact shown in Fig. 8 is calculated directly with about 22 CPU-hours on the M1 laptop, excluding training. The machine-learning controls designed by the LLM are used to reduce the estimated variance of the integration by about 80%80\%. Further optimization to stabilize the integration is needed and is left to future work. The P2B contribution ℛ3[J−jB]R_3[J-j_B] is obtained following Appendix D.4. For comparison, the green curves use the same contact-matching strategy as the neural NLO implementation. As a proof of concept, that calculation uses a full reference thrust cumulant at τ2≤0.3 _2≤ 0.3, giving c2=5.538±0.329c_2=5.538± 0.329. Both implementations use one contact for thrust, EEC, and the two jet energies, with no observable-dependent adjustment. For an EEC bin ℓ , the measurement is Jℓ=∑i,j(EiEj/Q2)ℓ(zij)J_ = _i,j(E_iE_j/Q^2)1_ (z_ij), where ℓ1_ is its indicator, zij=(1−cosχij)/2z_ij=(1- _ij)/2, and χij _ij is the opening angle. The sum includes i=ji=j. At Born level, this measurement gives weight 1/21/2 in each of the bins containing z=0z=0 and z=1z=1. Thrust has unit Born weight in the first bin, while both jet energies have unit weight in the bin containing E/Q=1/2E/Q=1/2. Jet energies are ranked after JADE/E clustering at ycut=0.04y_ cut=0.04, including all two-, three-, and four-jet events. These distributions give the NNLO coefficient, including the bins with Born support. The reference uses EERAD3 and restores the Born contribution with the known NNLO correction to the total cross section [69]. Figure 8 shows good agreement with the local-contact prediction, including the bins with Born support. The largest binwise differences are 2.362.36, 1.261.26, 1.831.83 and 1.281.28 combined standard errors for thrust, EEC, leading and second-leading jet energy, respectively. The agreement across these observables demonstrates the feasibility of recursive local contact extraction at NNLO. Similar to the NLO examples, larger samples and improved variance reduction provide a direct route to higher precision. Together, these results demonstrate two complementary uses of AI. Neural projections and fitted contacts reconstruct the full three-jet NLO coefficient, while analytic maps and direct contact integrals give local predictions for three and four jets at NLO and two jets at NNLO. The NNLO calculation also shows how machine learning can improve the integration of the constructed subtraction. Appendices B–D develop the three-jet NNLO extension and test its radiation contribution. VIII Conclusions and outlook We have developed infrared subtraction for multijet production with AI within the P2B + EFT matching framework. Human guidance specifies the framework and physical constraints, while the LLM devises the projections and contact formulas and implements the numerical calculation. The first implementation uses the expressive power of a neural network to construct a phase-space projection and determines its Born contact by fitting. Furthermore, the analytic implementation returns a local subtraction formula. Radiation is integrated at fixed Born momenta, and the contact is calculated directly from the EFT δ(τN)δ( _N) coefficient and finite four-dimensional radiation integrals. Both radiation and contact are evaluated over the full phase space without a small slicing parameter. The construction thus combines local subtraction with slicing’s advantage of recycling lower-order radiation calculations and EFT singular predictions. The numerical examples reconstruct the complete three-jet NLO coefficient in both implementations and the four-jet NLO coefficient with local subtraction, giving good agreement with EERAD3. The neural network implementation determines the contact through cumulant matching, while the analytic local calculations use the EFT δ(τN)δ( _N) coefficient and the directly evaluated radiation contribution 1A_1. At NNLO, recursive use of our NLO P2B construction gives the two-jet contact directly. Gaussian and polynomial machine-learning controls designed by the LLM reduce the variance of its finite radiation integral while preserving its value through exact integral compensation. The resulting predictions agree well with EERAD3. Machine learning therefore serves both to represent the subtraction ingredients and to improve their numerical integration. The calculations use the CPU of a 2020 M1 MacBook, without GPU acceleration, demonstrating that the construction is practical with modest local resources. For three jets at NNLO, Appendices B–D give explicit inverse maps and paired radiation weights that allow all Born momenta to remain fixed, together with a proposed finite-difference representation of the Born contact. The existing NLO matrix elements and subtraction terms are reused, while their phase-space integrations are reorganized at fixed Born momenta. Numerical checks test the inverse maps, radiation cancellation and integration at fixed Born kinematics. Appendix C extends this construction to any number of massless final-state jets and gives explicit maps and measures for radiation integration at fixed Born momenta. A complete study of the NNLO extension, including the numerical determination of 2A_2, its combination with the EFT coefficient of δ(τ3)δ( _3), and complete three-jet NNLO predictions, is left for future work. Further progress should focus on efficient evaluation of the local contact integrals at NNLO and higher perturbative orders. This calls for phase-space maps that preserve local cancellation and importance sampling adapted to the complete subtracted weight. LLM-designed machine-learning controls provide a complementary way to reduce variance while preserving the integral. Together with larger samples and parallel evaluation, these developments could make local subtraction practical for more demanding multijet calculations. Our results show how AI could help higher-order calculations by constructing infrared subtraction and making its numerical evaluation more efficient. Acknowledgements. We thank Xuan Chen for helpful discussions. This work is supported by the National Natural Science Foundation of China under Grant Nos. 12625505 and 12547109 and by the Fundamental Research Funds for the Central Universities, Beijing Normal University. X. L. would like to thank the Erwin-Schrödinger International Institute for Mathematics and Physics at the University of Vienna for partial support during the Programme ‘New Paradigms for Harnessing Quantum Field Theory at Colliders’, July 27–August 28, 2026. References [1] (2026) A counterexample to the Jacobian conjecture. Note: Technical note External Links: Link Cited by: §I. [2] K. Albertsson et al. (2018) Machine Learning in High Energy Physics Community White Paper. J. Phys. Conf. Ser. 1085 (2), p. 022008. External Links: 1807.02876, Document Cited by: §I. [3] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. -S. Shao, T. Stelzer, P. Torrielli, and M. Zaro (2014) The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations. JHEP 07, p. 079. External Links: 1405.0301, Document Cited by: §I. [4] C. Anastasiou, K. Melnikov, and F. Petriello (2004) A new method for real radiation at NNLO. Phys. Rev. D 69, p. 076010. External Links: hep-ph/0311311, Document Cited by: §I. [5] B. C. Aveleira, A. Gehrmann-De Ridder, T. Gehrmann, N. Glover, G. Heinrich, and C. T. Preuss (2025) EERAD3 version 2: QCD corrections in hadronic colour-singlet decays. SciPost Phys. Codeb. 59, p. 1. External Links: 2503.20610, Document Cited by: §B.4. [6] P. Banerjee, P. K. Dhani, and V. Ravindran (2018) Gluon jet function at three loops in QCD. Phys. Rev. D 98 (9), p. 094016. External Links: 1805.02637, Document Cited by: §D.2. [7] D. Baranowski, A. Behring, K. Melnikov, L. Tancredi, and C. Wever (2023) Beam functions for N-jettiness at N3^3LO in perturbative QCD. JHEP 2023 (02), p. 073. External Links: 2211.05722, Document Cited by: §D.2. [8] D. Baranowski, M. Delto, K. Melnikov, A. Pikelner, and C. Wang (2024) One-loop corrections to the double-real emission contribution to the zero-jettiness soft function at N3LO in QCD. JHEP 2024 (04), p. 114. External Links: 2401.05245, Document Cited by: §I. [9] D. Baranowski, M. Delto, K. Melnikov, A. Pikelner, and C. Wang (2025) Zero-Jettiness Soft Function to Third Order in Perturbative QCD. Phys. Rev. Lett. 134 (19), p. 191902. External Links: 2409.11042, Document Cited by: §I. [10] W. Bartel et al. (1986) Experimental Studies on Multi-Jet Production in e+ e- Annihilation at PETRA Energies. Z. Phys. C 33, p. 23. External Links: Document Cited by: §V.1. [11] C. L. Basham, L. S. Brown, S. D. Ellis, and S. T. Love (1978) Energy Correlations in electron - Positron Annihilation: Testing QCD. Phys. Rev. Lett. 41, p. 1585. External Links: Document Cited by: §VI. [12] C. W. Bauer, S. Fleming, and M. E. Luke (2000) Summing Sudakov logarithms in B→XsγB→ X_s effective field theory.. Phys. Rev. D 63, p. 014006. External Links: hep-ph/0005275, Document Cited by: §I. [13] C. W. Bauer, S. Fleming, D. Pirjol, and I. W. Stewart (2001) An Effective field theory for collinear and soft gluons: Heavy to light decays. Phys. Rev. D 63, p. 114020. External Links: hep-ph/0011336, Document Cited by: §I. [14] C. W. Bauer, D. Pirjol, and I. W. Stewart (2002) Soft collinear factorization in effective field theory. Phys. Rev. D 65, p. 054022. External Links: hep-ph/0109045, Document Cited by: §I. [15] C. W. Bauer and I. W. Stewart (2001) Invariant operators in collinear effective theory. Phys. Lett. B 516, p. 134–142. External Links: hep-ph/0107001, Document Cited by: §I. [16] T. Becher and G. Bell (2011) The gluon jet function at two-loop order. Phys. Lett. B 695, p. 252–258. External Links: 1008.1936, Document Cited by: §VII.2. [17] T. Becher and M. Neubert (2006) Toward a NNLO calculation of the anti-B —>> X(s) gamma decay rate with a cut on photon energy. I. Two-loop result for the jet function. Phys. Lett. B 637, p. 251–259. External Links: hep-ph/0603140, Document Cited by: §VII.2. [18] T. Becher and M. D. Schwartz (2008) A precise determination of αs _s from LEP thrust data using effective field theory. JHEP 07, p. 034. External Links: 0803.0342, Document Cited by: §VI, §VII.3. [19] G. Bell, A. Broggio, S. Edelmann, M. A. Lim, and R. Rahn (2026) NNLO soft functions for heavy-quark pair production at hadron colliders. External Links: 2608.09621 Cited by: §D.2. [20] G. Bell, K. Brune, G. Das, and M. Wald (2022) Automated Calculation of Beam Functions at NNLO. PoS L2022, p. 026. External Links: 2208.04847, Document Cited by: §D.2. [21] G. Bell, B. Dehnadi, T. Mohrmann, and R. Rahn (2024) The NNLO soft function for N-jettiness in hadronic collisions. JHEP 2024 (07), p. 077. External Links: 2312.11626, Document Cited by: §D.2. [22] G. Bell, R. Rahn, and J. Talbert (2019) Generic dijet soft functions at two-loop order: correlated emissions. JHEP 2019 (07), p. 101. External Links: 1812.08690, Document Cited by: §D.2. [23] M. Beneke, A. P. Chapovsky, M. Diehl, and T. Feldmann (2002) Soft collinear effective theory and heavy to light currents beyond leading power. Nucl. Phys. B 643, p. 431–476. External Links: hep-ph/0206152, Document Cited by: §I. [24] G. Bertolotti, G. Limatola, P. Torrielli, and S. Uccirati (2025) Advances in Local Analytic Sector Subtraction: massive NLO and elements of NNLO automation. External Links: 2503.14629 Cited by: §I. [25] G. Billis, B. Dehnadi, M. A. Ebert, J. K. L. Michel, and F. J. Tackmann (2021) Higgs pT Spectrum and Total Cross Section with Fiducial Cuts at Third Resummed and Fixed Order in QCD. Phys. Rev. Lett. 127 (7), p. 072001. External Links: 2102.08039, Document Cited by: §I. [26] E. Bothmann et al. (2019) Event Generation with Sherpa 2.2. SciPost Phys. 7 (3), p. 034. External Links: 1905.09127, Document Cited by: §I. [27] R. Boughezal, C. Focke, X. Liu, and F. Petriello (2015) W-boson production in association with a jet at next-to-next-to-leading order in perturbative QCD. Phys. Rev. Lett. 115 (6), p. 062002. External Links: 1504.02131, Document Cited by: §I. [28] R. Boughezal, X. Liu, and F. Petriello (2015) N-jettiness soft function at next-to-next-to-leading order. Phys. Rev. D 91 (9), p. 094035. External Links: 1504.02540, Document Cited by: §D.2. [29] R. Brüser, Z. L. Liu, and M. Stahlhofen (2018) Three-Loop Quark Jet Function. Phys. Rev. Lett. 121 (7), p. 072003. External Links: 1804.09722, Document Cited by: §D.2. [30] A. Butter, T. Plehn, and R. Winterhalder (2020) How to GAN Event Subtraction. SciPost Phys. Core 3, p. 009. External Links: 1912.08824, Document Cited by: §I. [31] M. Cacciari, F. A. Dreyer, A. Karlberg, G. P. Salam, and G. Zanderighi (2015) Fully Differential Vector-Boson-Fusion Higgs Production at Next-to-Next-to-Leading Order. Phys. Rev. Lett. 115 (8), p. 082002. Note: [Erratum: Phys.Rev.Lett. 120, 139901 (2018)] External Links: 1506.02660, Document Cited by: Appendix B, §I, §I, §VI. [32] F. Calisto, R. Moodie, and S. Zoia (2024) Learning Feynman integrals from differential equations with neural networks. JHEP 2024 (07), p. 124. External Links: 2312.02067, Document Cited by: §I. [33] S. Camarda, L. Cieri, and G. Ferrera (2021) Drell–Yan lepton-pair production: qT resummation at N3L accuracy and fiducial cross sections at N3LO. Phys. Rev. D 104 (11), p. L111503. External Links: 2103.04974, Document Cited by: §I. [34] J. Campbell, T. Neumann, and G. Vita (2025) Projection-to-Born-improved subtractions at NNLO. JHEP 05, p. 172. External Links: 2408.05265, Document Cited by: §I, §I. [35] F. Caola, K. Melnikov, and R. Röntsch (2017) Nested soft-collinear subtractions in NNLO QCD computations. Eur. Phys. J. C 77 (4), p. 248. External Links: 1702.01352, Document Cited by: §I. [36] S. Catani and M. H. Seymour (1997) A General algorithm for calculating jet cross-sections in NLO QCD. Nucl. Phys. B 485, p. 291–419. Note: [Erratum: Nucl.Phys.B 510, 503–504 (1998)] External Links: hep-ph/9605323, Document Cited by: §I, §I. [37] S. Catani, G. Turnock, and B. R. Webber (1992) Jet broadening measures in e+e−e^+e^- annihilation. Phys. Lett. B 295, p. 269–276. External Links: Document Cited by: §VII.1. [38] S. Catani and M. Grazzini (2007) An NNLO subtraction formalism in hadron collisions and its application to Higgs boson production at the LHC. Phys. Rev. Lett. 98, p. 222002. External Links: hep-ph/0703012, Document Cited by: §I. [39] W. Chen, F. Feng, Y. Jia, and X. Liu (2022) Double-real-virtual and double-virtual-real corrections to the three-loop thrust soft function. JHEP 2022 (12), p. 094. External Links: 2206.12323, Document Cited by: §I. [40] X. Chen, T. Gehrmann, E. W. N. Glover, A. Huss, B. Mistlberger, and A. Pelloni (2021) Fully Differential Higgs Boson Production to Third Order in QCD. Phys. Rev. Lett. 127 (7), p. 072002. External Links: 2102.07607, Document Cited by: §I. [41] X. Chen, D. Chicherin, E. Fox, N. Glover, M. Marcoli, V. Sotnikov, H. Sun, H. Zhang, and S. Zoia (2026) Four-Jet Rate in Electron-Positron Annihilation at Order αs4 _s^4. Phys. Rev. Lett. 136 (25), p. 251901. External Links: 2602.18185, Document Cited by: §C.3. [42] X. Chen, Y. Dai, H. T. Li, S. Li, H. Shao, and J. Wang (2026) Fully differential Higgs boson pair production at N3^3LO with top quark mass effects. JHEP 2026 (06), p. 005. External Links: 2601.19990, Document Cited by: §I. [43] X. Chen, T. Gehrmann, E. W. N. Glover, A. Huss, P. F. Monni, E. Re, L. Rottoli, and P. Torrielli (2022) Third-Order Fiducial Predictions for Drell-Yan Production at the LHC. Phys. Rev. Lett. 128 (25), p. 252001. External Links: 2203.01565, Document Cited by: §I. [44] X. Chen, T. Gehrmann, N. Glover, A. Huss, T. Yang, and H. X. Zhu (2022) Dilepton Rapidity Distribution in Drell-Yan Production to Third Order in QCD. Phys. Rev. Lett. 128 (5), p. 052001. External Links: 2107.09085, Document Cited by: §I. [45] X. Chen, T. Gehrmann, N. Glover, A. Huss, T. Yang, and H. X. Zhu (2023) Transverse mass distribution and charge asymmetry in W boson production to third order in QCD. Phys. Lett. B 840, p. 137876. External Links: 2205.11426, Document Cited by: §I. [46] X. Chen, P. Jakubčík, M. Marcoli, and G. Stagnitto (2025) Jet production at electron-positron colliders at next-to-next-to-next-to-leading order in QCD. Phys. Lett. B 869, p. 139804. External Links: 2505.10618, Document Cited by: §I. [47] L. Cieri, X. Chen, T. Gehrmann, E. W. N. Glover, and A. Huss (2019) Higgs boson production at the LHC using the qTq_T subtraction formalism at N3^3LO QCD. JHEP 2019 (02), p. 096. External Links: 1807.11501, Document Cited by: §I. [48] Claude (Anthropic) (2026) Nine-loop six-gluon MHV amplitude in planar =4N=4 super-Yang–Mills theory: large data files. Note: Zenodo datasetResults and validation External Links: Document, Link Cited by: §I. [49] J. Currie, T. Gehrmann, E. W. N. Glover, A. Huss, J. Niehues, and A. Vogt (2018) N3^3lo corrections to jet production in deep inelastic scattering using the projection-to-born method. JHEP 2018 (05), p. 209. External Links: 1803.09973, Document Cited by: §I, §VI. [50] G. Cybenko (1989) Approximation by superpositions of a sigmoidal function. Math. Control Signals Syst. 2 (4), p. 303–314. External Links: Document Cited by: §I. [51] M. Czakon (2010) A novel subtraction scheme for double-real radiation at NNLO. Phys. Lett. B 693, p. 259–268. External Links: 1005.0274, Document Cited by: §I. [52] S. Diefenbacher, A. Hallin, G. Kasieczka, M. Krämer, A. Lauscher, and T. Lukas (2025) Agents of Discovery. External Links: 2509.08535 Cited by: §I. [53] L. Dong, S. Fang, J. Gao, H. T. Li, D. Y. Shao, H. X. Zhu, and Y. J. Zhu (2026) Two-Dimensional Transverse-Momentum Subtraction and Semi-Inclusive Deep-Inelastic Scattering at N3^3LO in QCD. External Links: 2603.29673 Cited by: §I. [54] M. A. Ebert, B. Mistlberger, and G. Vita (2020) N-jettiness beam functions at N3^3LO. JHEP 2020 (09), p. 143. External Links: 2006.03056, Document Cited by: §D.2. [55] M. A. Ebert, B. Mistlberger, and G. Vita (2020) Transverse momentum dependent PDFs at N3^3LO. JHEP 2020 (09), p. 146. External Links: 2006.05329, Document Cited by: §D.2. [56] M. A. Ebert, B. Mistlberger, and G. Vita (2021) TMD fragmentation functions at N3^3LO. JHEP 2021 (07), p. 121. External Links: 2012.07853, Document Cited by: §D.2. [57] M. A. Ebert and F. J. Tackmann (2020) Impact of isolation and fiducial cuts on qT_T and N-jettiness subtractions. JHEP 03, p. 158. External Links: 1911.08486, Document Cited by: §I, §I. [58] J. S. Ellenberg, C. S. Fraser-Taliente, T. R. Harvey, K. Srivastava, and A. V. Sutherland (2025) Generative Modeling for Mathematical Discovery. External Links: 2503.11061 Cited by: §I. [59] R. K. Ellis, D. A. Ross, and A. E. Terrano (1981) The perturbative calculation of jet structure in e+e−e^+e^- annihilation. Nucl. Phys. B 178, p. 421–456. External Links: Document Cited by: §I. [60] R. K. Ellis, W. J. Stirling, and B. R. Webber (1996) QCD and collider physics. Cambridge University Press. Note: Chap. 3, p. 58 and 72 External Links: Document Cited by: §I. [61] E. Farhi (1977) A QCD Test for Jets. Phys. Rev. Lett. 39, p. 1587–1588. External Links: Document Cited by: §VI. [62] R. Frederix, S. Frixione, F. Maltoni, and T. Stelzer (2009) Automation of next-to-leading order computations in QCD: The FKS subtraction. JHEP 2009 (10), p. 003. External Links: 0908.4272, Document Cited by: §I. [63] S. Frixione, Z. Kunszt, and A. Signer (1996) Three jet cross-sections to next-to-leading order. Nucl. Phys. B 467, p. 399–442. External Links: hep-ph/9512328, Document Cited by: §I, §I. [64] R. Fu, R. Rahn, D. Y. Shao, W. J. Waalewijn, and B. Wu (2025) qT Slicing with Multiple Jets. Phys. Rev. Lett. 135 (17), p. 171903. External Links: 2412.05358, Document Cited by: §I. [65] H. Gao et al. (2026) LQCDMaster: Agentic Scientific Computing for Lattice Quantum Chromodynamics Research. External Links: 2607.15001 Cited by: §I. [66] J. Gao, C. S. Li, and H. X. Zhu (2013) Top Quark Decay at Next-to-Next-to Leading Order in QCD. Phys. Rev. Lett. 110 (4), p. 042001. External Links: 1210.2808, Document Cited by: §I, §I. [67] J. Gaunt, M. Stahlhofen, F. J. Tackmann, and J. R. Walsh (2015) N-jettiness Subtractions for NNLO QCD Calculations. JHEP 2015 (09), p. 058. External Links: 1505.04794, Document Cited by: Appendix A, §I, §I, §I, §V.1, §VII.1. [68] T. Gehrmann, A. Huss, J. Niehues, A. Vogt, and D. M. Walker (2019) Jet production in charged-current deep-inelastic scattering to third order in QCD. Phys. Lett. B 792, p. 182–186. External Links: 1812.06104, Document Cited by: §I. [69] A. Gehrmann-De Ridder, T. Gehrmann, E. W. N. Glover, and G. Heinrich (2009) NNLO moments of event shapes in e+e- annihilation. JHEP 2009 (05), p. 106. External Links: 0903.4658, Document Cited by: §VII.3. [70] A. Gehrmann-De Ridder, T. Gehrmann, E. W. N. Glover, and G. Heinrich (2014) EERAD3: Event shapes and jet rates in electron-positron annihilation at order αs3 _s^3. Comput. Phys. Commun. 185, p. 3331. External Links: 1402.4140, Document Cited by: Appendix A, §V, §VII.1. [71] A. Gehrmann-De Ridder, T. Gehrmann, and E. W. N. Glover (2005) Antenna subtraction at NNLO. JHEP 2005 (09), p. 056. External Links: hep-ph/0505111, Document Cited by: §B.1, §B.2, §C.1, Appendix C, §I, §VII.2. [72] W. T. Giele, E. W. N. Glover, and D. A. Kosower (1993) Higher order corrections to jet cross-sections in hadron colliders. Nucl. Phys. B 403, p. 633–670. External Links: hep-ph/9302225, Document Cited by: §I. [73] B. W. Harris and J. F. Owens (2002) The Two cutoff phase space slicing method. Phys. Rev. D 65, p. 094032. External Links: hep-ph/0102128, Document Cited by: §I. [74] A. Huss et al. (2026) NNLOJET: A parton-level event generator for jet cross sections at NNLO QCD accuracy. SciPost Phys. Codeb. 69, p. 1. External Links: 2503.22804, Document Cited by: §I. [75] T. Janßen, R. Poncelet, and S. Schumann (2025) Sampling NNLO QCD phase space with normalizing flows. JHEP 2025 (09), p. 194. External Links: 2505.13608, Document Cited by: §I. [76] T. T. Jouttenus, I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn (2011) The Soft Function for Exclusive N-Jet Production at Hadron Colliders. Phys. Rev. D 83, p. 114030. External Links: 1102.4344, Document Cited by: §I, §V.1, §VII.1, §VII.2. [77] R. Kelley, M. D. Schwartz, R. M. Schabinger, and H. X. Zhu (2011) The two-loop hemisphere soft function. Phys. Rev. D 84, p. 045022. External Links: 1105.3676, Document Cited by: §VI, §VII.3. [78] T. Kinoshita (1962) Mass singularities of Feynman amplitudes. J. Math. Phys. 3, p. 650–677. External Links: Document Cited by: §I. [79] A. J. Larkoski, D. Neill, and J. Thaler (2014) Jet Shapes with the Broadening Axis. JHEP 04, p. 017. External Links: 1401.2158, Document Cited by: §I. [80] T. D. Lee and M. Nauenberg (1964) Degenerate Systems and Mass Singularities. Phys. Rev. 133, p. B1549–B1562. External Links: Document Cited by: §I. [81] G. P. Lepage (1978) A New Algorithm for Adaptive Multidimensional Integration. J. Comput. Phys. 27, p. 192. External Links: Document Cited by: §VII.1. [82] Y. Liu, Y. Xu, and Y. Zhang (2025) Uncovering Singularities in Feynman Integrals via Machine Learning. External Links: 2510.10099 Cited by: §I. [83] M. Luo, T. Yang, H. X. Zhu, and Y. J. Zhu (2020) Quark Transverse Parton Distribution at the Next-to-Next-to-Next-to-Leading Order. Phys. Rev. Lett. 124 (9), p. 092001. External Links: 1912.05778, Document Cited by: §D.2. [84] M. Luo, T. Yang, H. X. Zhu, and Y. J. Zhu (2021) Unpolarized quark and gluon TMD PDFs and FFs at N3^3LO. JHEP 2021 (06), p. 115. External Links: 2012.03256, Document Cited by: §D.2. [85] L. Magnea, E. Maina, G. Pelliccioli, C. Signorile-Signorile, P. Torrielli, and S. Uccirati (2018) Local analytic sector subtraction at NNLO. JHEP 2018 (12), p. 107. Note: [Erratum: JHEP 06, 013 (2019)] External Links: 1806.09570, Document Cited by: §I. [86] D. Maître and H. Truong (2021) A factorisation-aware matrix element emulator. JHEP 2021 (11), p. 066. External Links: 2107.06625, Document Cited by: §I. [87] M. Marcoli (2026) Recent progress in antenna subtraction at NNLO and N3^3LO. External Links: 2603.05747 Cited by: §I. [88] T. Menzo, A. Roman, S. Gleyzer, K. Matchev, G. T. Fleming, S. Höche, S. Mrenna, and P. Shyamsundar (2025) HEPTAPOD: Orchestrating High Energy Physics Workflows Towards Autonomous Agency. External Links: 2512.15867 Cited by: §I. [89] R. Mondini, M. Schiavi, and C. Williams (2019) N3^3LO predictions for the decay of the Higgs boson to bottom quarks. JHEP 2019 (06), p. 079. External Links: 1904.08960, Document Cited by: §I. [90] T. Neumann and J. Campbell (2023) Fiducial Drell-Yan production at the LHC improved by transverse-momentum resummation at N4LLp+N3LO. Phys. Rev. D 107 (1), p. L011506. External Links: 2207.07056, Document Cited by: §I. [91] OpenAI (2026) Finite time blowup for Navier–Stokes. External Links: Link Cited by: §I. [92] Z. Song, T. Yang, Q. Cao, M. Luo, and H. X. Zhu (2026) Explainable AI-assisted optimization for Feynman integral reduction. JHEP 2026 (06), p. 225. External Links: 2502.09544, Document Cited by: §I. [93] G. F. Sterman and S. Weinberg (1977) Jets from Quantum Chromodynamics. Phys. Rev. Lett. 39, p. 1436. External Links: Document Cited by: §I. [94] I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn (2010) N-Jettiness: An Inclusive Event Shape to Veto Jets. Phys. Rev. Lett. 105, p. 092002. External Links: 1004.2489, Document Cited by: §I, §V.1. [95] M. von Hippel and M. Wilhelm (2025) Refining Integration-by-Parts Reduction of Feynman Integrals with Machine Learning. JHEP 2025 (05), p. 185. External Links: 2502.05121, Document Cited by: §I. Appendix A An analytic NLO projection This appendix gives the analytic projection and its inverse used in the local NLO calculation. The projection maps a real configuration to Born kinematics; the inverse generates radiation while keeping all Born momenta fixed. We denote the projection by M^:ΦN+1→ΦN M: _N+1→ _N. It must preserve on-shell kinematics and four-momentum, and approach the underlying Born state in every single-unresolved limit. The tree-level real coefficient RN+1R_N+1 must give an integrable product with ΔJ=JN+1−JN∘M J=J_N+1-J_N M, with physical cuts evaluated on each state separately. For massless final states with P=(Q,)P=(Q,0), choose the pair h=(i,j)h=(i,j) minimizing dijd_ij and define Kij K_ij =pi+pj,dij=Ei+Ej−|ij|, =p_i+p_j, d_ij=E_i+E_j-|K_ij|, (15) λ λ =1−dij/Q, =1-d_ij/Q, p¯1 p_1 =(|ij|,ij)λ,p¯a=psaλ(a=2,…,N), = (|K_ij|,K_ij)λ, p_a= p_s_aλ (a=2,…,N), where sas_a labels the spectators in increasing parton-label order; the merged pair is assigned Born label 1. We denote this projection for a specified pair h=(i,j)h=(i,j) by M^h M_h; the global map M M selects the pair minimizing dijd_ij, with a fixed ordering used to break ties. The projected momenta obey p¯a2=0 p_a^2=0 and ∑ap¯a=P _a p_a=P. In a soft or collinear limit, the selected pair contains the unresolved radiation and λ→1λ→ 1, recovering the Born momenta. For a measurement Lipschitz-continuous in energy flow, |J(X)−J(Y)|≤LJdE(X,Y)|J(X)-J(Y)|≤ L_Jd_E(X,Y) for a finite constant LJL_J, where dE(X,Y)=Q−1min∑abfabfab|a−b′|d_E(X,Y)=Q^-1 _\f_ab\ _abf_ab|n_a-n _b| is the minimum energy-transport cost between states of total energy Q. Here fab≥0f_ab≥ 0 transfers energy from direction an_a in X to b′n _b in Y, with ∑bfab=Ea _bf_ab=E_a and ∑afab=Eb′ _af_ab=E _b, the corresponding particle energies. For thrust and the C-parameter, this gives |ΔJ|=O(dij/Q)| J|=O( d_ij/Q). In a soft-collinear limit, dij/Q=O(yθij2)d_ij/Q=O(y _ij^2), suppressing dy/ydθij2/θij2dy/y\,d _ij^2/ _ij^2; here y is the soft energy divided by Q and θij _ij is the pair’s opening angle. Conversely, fix b=(p¯1,…,p¯N)b=( p_1,…, p_N) with e1=p¯10e_1= p_1^0 and choose 0<λ<10<λ<1. Only the N−1N-1 spectator momenta are rescaled. The Born leg p¯1 p_1 is replaced by two massless particles with total momentum KijK_ij, fixed by four-momentum conservation, psa p_s_a =λp¯a,Kij=(1−λ)P+λp¯1, =λ p_a, K_ij=(1-λ)P+λ p_1, (16a) Kij2 K_ij^2 =(1−λ)[Q2(1−λ)+2λQe1]>0. =(1-λ) [Q^2(1-λ)+2λ Qe_1 ]>0. Write p¯1=e1(1,) p_1=e_1(1,n), and k=λe1k=λ e_1. Choose orthonormal vectors x,ye_x,e_y transverse to n and set ⟂(ϕ)=cosϕx+sinϕye_ (φ)= φ\,e_x+ φ\,e_y. The individual momenta are pi0 p_i^0 =Kij0+kcosθ2,pj=Kij−pi, = K_ij^0+k θ2, p_j=K_ij-p_i, (16b) i _i =12[(k+Kij0cosθ)+Kij2sinθ⟂(ϕ)]. = 12 [(k+K_ij^0 θ)n+ K_ij^2 θ\,e_ (φ) ]. Here θ,ϕθ,φ specify the direction of pip_i in the KijK_ij rest frame, with polar axis n, −1<cosθ<1-1< θ<1 and 0≤ϕ<2π0≤φ<2π. Equation (16b) boosts the two back-to-back massless daughters to the total centre-of-mass frame, ensuring pi2=pj2=0p_i^2=p_j^2=0. Define Gh(b,ξ)G_h(b,ξ) as the complete event pi,pj,ps2,…,psN\p_i,p_j,p_s_2,…,p_s_N\, with ξ=(λ,cosθ,ϕ)ξ=(λ, θ,φ). Reapplying the specified-pair projection gives M^h(Gh(b,ξ))=b M_h(G_h(b,ξ))=b. The selector χh _h is one when h is the pair chosen by M M and zero otherwise. Selecting charts with χh=1 _h=1 therefore ensures M^(Gh(b,ξ))=b M(G_h(b,ξ))=b. Before this selection, the labelled Lorentz-invariant measure on each chart is dΦN+1=dΦN(b)Qe1λ2N−332π3⏟h(b,ξ)dλdcosθdϕ.d _N+1=d _N(b)\, Qe_1λ^2N-332π^3_J_h(b,ξ)dλ\,d θ\,dφ. (17) Summing branches with χh _h covers phase space without double counting. Adding and subtracting the projected real measurement gives the NLO correction ΣN(1)[J]=∫dΦN+1RN+1ΔJ+∫dΦN(b)cM^(b)JN(b), _N^(1)[J]= _N+1R_N+1 J+ _N(b)c_ M(b)J_N(b), (18) where cM^c_ M combines the virtual and projected real contributions into a finite density relative to dΦNd _N. For the pointwise matching below, M M in Eq. (18) must be the projection in Eq. (15), with inverse and measure given by Eqs. (16)–(17). Direct sampling in ΦN+1 _N+1 needs no additional Jacobian, but does not permit changing this projection while reusing the contact obtained from Eq. (19). The inverse map allows cM^(b)c_ M(b) to be determined separately at every resolved N-body Born phase-space point b. The required radiation input is the tree-level jettiness cumulant above t, evaluated at this fixed Born point, cM^(b) c_ M(b) =limt→0[ΣN,sing(1)(b,t)+ℛ>t(b)], = _t→ 0 [ _N,sing^(1)(b,t)+R_>t(b) ], (19) ℛ>t(b) _>t(b) =∫dΦN+1RN+1Θ(τN−t)|b. = . _N+1\,R_N+1 ( _N-t) |_b. The notation |b|_b denotes a density relative to dΦN(b)d _N(b). Generate the real momenta with Gh(b,ξ)G_h(b,ξ), select the branch with χh=1 _h=1, and integrate only the radiation variables with the Jacobian in Eq. (17). One then evaluates τN _N on each generated state and includes its weight when τN>t _N>t. For the geometric measure in Eq. (10), this is simply λ<1−tλ<1-t because τN=1−λ _N=1-λ. The Born point b is held fixed throughout. Here ΣN,sing(1)(b,t) _N,sing^(1)(b,t) is the NLO correction to the EFT cumulant for τN≤t _N≤ t, differential in the same Born variables at leading power and including the coefficient of δ(τN)δ( _N) [67]. Its addition to ℛ>tR_>t compensates the negative radiation contribution to the below-t cumulant. Equation (19) then determines cM^(b)c_ M(b) point by point. This limiting relation fixes the contact convention. Appendix D evaluates the same contact by finite four-dimensional integrals, so neither its extraction nor Eq. (18) requires a small slicing parameter. For N=3N=3, Table A confirms the kinematic identities and phase-space normalization using momenta and phase-space weights alone. Table 2: Phase-space checks of the analytic 4→34→ 3 map. Residuals are maximum componentwise errors, with momenta normalized to Q and squared masses to Q2Q^2. The scans use 4096 points per branch and 4096 independent four-body points. The volume test uses 2282^28 points in 32 independent batches; its uncertainty is one standard error. Check Result On-shell and momentum conservation <4×10−14<4× 10^-14 Complete four-body round trip <5×10−13<5× 10^-13 Fixed-Born reconstruction <3×10−16<3× 10^-16 Volume / [Q4/(24576π5)][Q^4/(24576π^5)] 1.00006±0.000081.00006± 0.00008 We now test local cancellation with the four-parton tree matrix elements from EERAD3 [70]. Fix a resolved Born point with 2Ei/Q=(0.8,0.7,0.5)2E_i/Q=(0.8,0.7,0.5), for which 1−T=0.21-T=0.2 and C=9/14C=9/14. Choose ¯1 p_1 along +z+z and ¯2 p_2 in the xzxz plane with positive x, fixing the overall orientation. For either measurement J=1−T,CJ=1-T,C, the radiation integral normalized to the three-parton Born density B3(b)B_3(b) is ℱJ(b,η) _J(b;η) =1B3(b)∑h∫dξhχhR4(Gh) = 1B_3(b) _h ξ\,J_h _hR_4(G_h) (20) ×[J4(Gh)−J3(b)]Θ(ymin−η), ×[J_4(G_h)-J_3(b)]\, (y_ -η), where ymin=mini<jsij/Q2y_ = _i<js_ij/Q^2. Replacing the measurement difference by J4(Gh)J_4(G_h) or J3(b)J_3(b) defines the separate real and projected integrals, ℱJRF_J^R and ℱJPF_J^P, so that ℱJ=ℱJR−ℱJPF_J=F_J^R-F_J^P. The powers of αs/(2π) _s/(2π) are stripped from the tree coefficients R4R_4 and B3B_3. All radiation variables are integrated using Eqs. (16)–(17), while b remains fixed. Figure 9: Local P2B cancellation at a fixed three-parton Born point. The upper panels show the qq¯ggq qg channel with the qgqg pair (1,3)(1,3) generated by Eq. (16), at ϕ=0.7φ=0.7. Write δ=1−λδ=1-λ and w=1−cosθw=1- θ. The collinear, soft and soft-collinear paths use (δ,w)=(ϵ,0.8),(ϵ,20ϵ),(ϵ2,2ϵ)(δ,w)=(ε,0.8),(ε,20ε),(ε^2,2ε), respectively. The ordinate JW_J denotes the logarithmic radiation density for J=1−TJ=1-T. Blue and orange show ωhR4J4/B3 _hR_4J_4/B_3 and ωhR4J3(b)/B3 _hR_4J_3(b)/B_3; green shows their difference. Here ω=δω=δ in (a) and ω=δwω=δ w in (b,c), corresponding to logarithmic radiation coordinates. The lower panels integrate all six branches and both partonic channels. Their lower strips show ℱJ(η)/ℱJrefF_J(η)/F_J^ref, with ℱJref=ℱJ(10−10)F_J^ref=F_J(10^-10) and correlated statistical errors. Grey bands give the relative one-standard-error uncertainty of the reference integral. The technical cut η on pair invariants probes the numerical stability of the locally subtracted integral. Figure 9 shows the local cancellation and the convergence of this integral. The separate real and projected logarithmic radiation densities approach nonzero limits, while their difference vanishes in soft, collinear and soft-collinear scans. The integral includes both qq¯ggq qg and four-quark channels with nf=5n_f=5, using 3232 independent batches of 6×2166× 2^16 radiation points. Every accepted configuration projects back to the same b, with maximum componentwise error below 2×10−16Q2× 10^-16Q. At η=10−10η=10^-10 we obtain ℱ1−T=1.1372±0.0013F_1-T=1.1372± 0.0013 and ℱC=3.0338±0.0035F_C=3.0338± 0.0035. Reducing η from 10−610^-6 to 10−1010^-10 changes either result by less than 0.011%0.011\%, while the separate real and projected integrals continue to grow. These tests demonstrate local infrared cancellation with the Born phase-space point held fixed throughout. Appendix B A proposal for local three-jet NNLO subtraction at fixed Born kinematics This appendix develops the radiation part of the proposed extension of the local subtraction to three-jet NNLO, using a four-jet NLO calculation as input [31]. To preserve local cancellation at a specified three-parton Born point, the real contribution and its counterterm must be generated on a common phase space with those Born momenta fixed. The LLM constructs it by composing the inverse in Appendix A with the inverse of the NLO subtraction map. We first give the phase space and paired radiation weights, then discuss their integrability and numerical checks. Appendix D gives the complementary proposal for direct contact calculation. Let R5R_5 denote the five-parton tree contribution, SaS_a the counterterms grouped by their maps m^a:Φ5→Φ4 m_a: _5→ _4, and K4K_4 the sum of the renormalized four-parton one-loop contribution and integrated counterterms, finite at resolved four-parton kinematics. For an infrared-safe four-jet measurement fmf_m on an m-parton state, the NLO correction is ℐ4(1)[f] _4^(1)[f] =∫dΦ5[R5f5−∑aSaf4(m^aΦ5)] = _5 [R_5f_5- _aS_af_4( m_a _5) ] (21) +∫dΦ4K4f4. + _4\,K_4f_4. Only this correction, at the order of the three-jet NNLO coefficient, is needed. The matrix elements, counterterms, maps and finite K4K_4 are reused without modification, while their four- and five-parton integrations are replaced by the common phase space below. No new integrated counterterms or d-dimensional calculation are required. Flavour sums, symmetry factors and normalization are inherited from the NLO input, whose generation cuts must cover the new measurement. B.1 A common phase space We illustrate the pairing with antenna subtraction [71]; other NLO schemes can follow the same strategy with compatible inverse maps, measures and infrared bounds. Use the P2B projection M^4:Φ4→Φ3 M_4: _4→ _3 of Eq. (15), ordering its Born output by decreasing energy. The inverse sums six pair choices and six assignments of the Born momenta to the merged leg and spectators. Here h includes both choices. Write G4,hG_4,h for the N=3N=3 inverse in Eqs. (16)–(17), generating x=G4,h(b,ξ)x=G_4,h(b,ξ). At fixed b, its radiation measure is dμ4b=∑hχhhdξd _4^b= _h _hJ_hdξ. The Born measure dΦ3(b)d _3(b) is restricted to this energy ordering. Denote the two angles in ξ by θ4,ϕ4 _4, _4, to distinguish them from the antenna azimuth ϕφ. For each antenna channel, write x=(I,K,q1,q2)x=(I,K,q_1,q_2), where I,KI,K are the massless parent four-momenta replacing partons i,j,ki,j,k under m^a m_a, and q1,q2q_1,q_2 are unchanged spectators. In the I+KI+K rest frame, divide momenta by sIK s_IK, where sIK=(I+K)2s_IK=(I+K)^2, giving a dimensionless total (1,)(1,0). For r,s>0r,s>0 and r+s<1r+s<1, choose the massless daughters p^i p_i =12(1−s)(1,0,0,1), = 12(1-s)(1,0,0,1), (22) p^j p_j =(r+s2,rs(1−r−s)1−s,0,r+s2−r1−s), = ( r+s2, rs(1-r-s)1-s,0, r+s2- r1-s ), p^k p_k =(1,)−p^i−p^j. =(1,0)- p_i- p_j. The hatted momenta satisfy 2p^i⋅p^j=r2 p_i· p_j=r and 2p^j⋅p^k=s2 p_j· p_k=s. The antenna subtraction map gives the back-to-back parents c c =1−r−s,y=sr+s, =1-r-s, y= sr+s, (23) ρant _ant =1+4y(1−y)rsc, = 1+ 4y(1-y)rsc, ai a_i =1+ρant−2ys2(1−s),ak=1−ρant−2yr2(1−r), = 1+ _ant-2ys2(1-s), a_k= 1- _ant-2yr2(1-r), I I =aip^i+yp^j+akp^k, =a_i p_i+y p_j+a_k p_k, K K =(1,)−I^. =(1,0)- I. Rotate I I onto the direction of I in the I+KI+K rest frame, rotate the decay plane about that direction by ϕφ, multiply all daughters by sIK s_IK, and boost from the I+KI+K rest frame to the total centre-of-mass frame. The two spectators are unchanged. This defines the inverse p=Gaant(x,r,s,ϕ)p=G_a^ant(x;r,s,φ) explicitly. For 0≤ϕ<2π0≤φ<2π, define the azimuthal origin using the spectator with the larger component transverse to I in the antenna rest frame. If both components vanish, use any nonparallel Cartesian reference axis. The resulting phase space obeys m^aGaant(x,r,s,ϕ) m_aG_a^ant(x;r,s,φ) =x,M^5,a=M^4∘m^a, =x, M_5,a= M_4 m_a, (24) M^5,aGaant(G4,h(b,ξ),r,s,ϕ) M_5,aG_a^ant(G_4,h(b,ξ);r,s,φ) =b, =b, dΦ5 _5 =dΦ4(x)sIK32π3drdsdϕ⏟dΦant,a. =d _4(x)\, s_IK32π^3dr\,ds\,dφ_d _ant,a. Thus all three Born momenta are fixed exactly in both members of a subtraction pair. The phase-space Jacobian is the product of the measure in Appendix A and the antenna measure. Table 3 checks this construction against the NLO subtraction maps. Table 3: Kinematic checks of the common fixed-Born phase space at 2Ei/Q=(0.8,0.7,0.5)2E_i/Q=(0.8,0.7,0.5). Of 600 attempted four-parton configurations, 308 satisfy the selected-chart condition. Each is tested at azimuths ϕφ and ϕ+π/2φ+π/2. Entries are maximum componentwise residuals normalized to Q, or Q2Q^2 for the mass shell, using extended-precision geometry. Check Maximum residual Mass shell <7×10−17<7× 10^-17 Momentum conservation <2×10−16<2× 10^-16 Subtraction inverse, m^aGaant=x m_aG_a^ant=x <3×10−13<3× 10^-13 Fixed Born point, M^5,aGaantG4,h=b M_5,aG_a^antG_4,h=b <3×10−16<3× 10^-16 B.2 Local subtraction on the common phase space Different counterterms have different maps. We therefore partition the real kernel into contributions RaR_a and assign each the projection M^5,a=M^4∘m^a M_5,a= M_4 m_a. An additional smooth partition suppresses counterterms whose reduced four-parton state becomes singular while the real state remains resolved. Let va=τ3(m^ap)v_a= _3( m_ap), with the minimized geometric measure of Eq. (10), and define the dimensionless transport distance between the three daughters and their two parents by ℓa=1Qmin∑l=i,j,k0≤zl≤El∑lzl=EI[zl|l−I|+(El−zl)|l−K|]. _a= 1Q _ subarrayc0≤ z_l≤ E_l\\ _lz_l=E_I subarray _l=i,j,k [z_l|n_l-n_I|+(E_l-z_l)|n_l-n_K| ]. (25) All energies and unit directions in Eq. (25) refer to the total centre-of-mass frame. The variable zlz_l is the energy transferred from daughter l to parent I; the remainder goes to K. Sorting |l−I|−|l−K||n_l-n_I|-|n_l-n_K| and assigning energy to I in that order until EI=I0E_I=I^0 is reached solves the minimization. It bounds the change of an energy-flow Lipschitz measurement, |J5−J4|≤LJℓa|J_5-J_4|≤ L_J _a, and also |τ3(p)−va|≤ℓa| _3(p)-v_a|≤ _a. Choose ha h_a =h(vava+ℓa), =h\! ( v_av_a+ _a ), (26) h(z) h(z) =0,z≤12,w2(3−2w),w=4z−2,12<z<34,1,z≥34. = cases0,&z≤ 12,\\ w^2(3-2w), w=4z-2,& 12<z< 34,\\ 1,&z≥ 34. cases Near a single-unresolved limit of channel a with resolved x, ha=1h_a=1. It vanishes when va→0v_a→ 0 while τ3(p) _3(p) stays nonzero, suppressing a counterterm singularity absent from the real configuration. The complementary contribution is retained, so this partition removes no phase-space region. A suitable signed partition of the real contribution is d^ij d_ij =sijQ(Ei+Ej+|i+j|),ua=min(d^ij,d^jk), = s_ijQ(E_i+E_j+|p_i+p_j|), u_a= ( d_ij, d_jk), (27) ωa _a =ua−2∑cuc−2,D=R5−∑chcSc, = u_a^-2 _cu_c^-2, D=R_5- _ch_cS_c, Ra R_a =haSa+ωaD,∑aRa=R5. =h_aS_a+ _aD, _aR_a=R_5. Here d^ij=dij/Q d_ij=d_ij/Q is the dimensionless deficit corresponding to Eq. (15). The RaR_a need not be positive. The 30 channels have three distinct labels (i,j,k)(i,j,k) with i<ki<k; j is the parton removed by the subtraction map. We average R5R_5 and SaS_a over the same 5!5! label permutations, grouping counterterms with identical maps, and average K4K_4 over 4!4! permutations. These averages relabel integration variables and introduce no extra multiplicity. On an accepted chart, x=G4,h(b,ξ)x=G_4,h(b,ξ) and p=Gaant(x,r,s,ϕ)p=G_a^ant(x;r,s,φ), set J5=J(p)J_5=J(p), J4=J(x)J_4=J(x) and ΔJm=Jm−J3(b) J_m=J_m-J_3(b). The common local weight is WaJ W_a^J =RaΔJ5−SaΔJ4 =R_a J_5-S_a J_4 (28) =ωaDΔJ5+haSa(J5−J4)−(1−ha)SaΔJ4. = _aD J_5+h_aS_a(J_5-J_4)-(1-h_a)S_a J_4. Every difference is formed before integration, at the same five-parton point and its mapped four-parton state. Scalar antennae require an additional treatment of collinear spin correlations [71]. At fixed x,r,sx,r,s, define ⟨WaJ⟩ϕ=[WaJ(ϕ)+WaJ(ϕ+π/2)]/2 W_a^J _φ=[W_a^J(φ)+W_a^J(φ+π/2)]/2. This average preserves the full azimuthal integral. The leading cos2ϕ 2φ and sin2ϕ 2φ harmonics cancel within this pair, which is also the unit for estimating statistical errors. At fixed b, the radiation contribution and its combination with the required NNLO contact take the form FJ(b) F_J^P(b) =∑a∫dμ4bdΦant,a⟨WaJ⟩ϕ+∫dμ4bK4ΔJ4, = _a _4^b\,d _ant,a\, W_a^J _φ+ _4^b\,K_4 J_4, (29) Σ3(2)[J] _3^(2)[J] =∫dΦ3(b)[FJ(b)+c(2)(b)J3(b)]. = _3(b) [F_J^P(b)+c_P^(2)(b)J_3(b) ]. Here P denotes the combined choice of M^4 M_4, m^a m_a and RaR_a. Integrating over b and using ∑aRa=R5 _aR_a=R_5 recovers the real and four-parton terms of Eq. (21); the remaining projected contribution is absorbed into c(2)c_P^(2). For the constant measurement J=1J=1, F1=0F_1^P=0 sample by sample. The four-parton weight K4ΔJ4K_4 J_4 is separately integrable. B.3 Integrability and pointwise contact matching For a fixed Born point b away from soft and collinear boundaries, the maps give the bounds dE(x,b)≤vad_E(x,b) v_a and ℓa≤min(r,s) _a (r,s), with dEd_E defined in Appendix A. The positive constants C, ′C , JC_J and c0c_0 below remain finite on compact resolved Born regions. The first bound follows by transporting the merged pair onto its total direction and rescaling; the second follows from the explicit antenna inverse. Combining the latter bound with ℓa≤sIK/Q _a s_IK/Q and min(sij,sjk)/Q2≤2ua (s_ij,s_jk)/Q^2≤ 2u_a gives ℓa≤ua1/4 _a _a^1/4. With umin=minaua≤v=τ3(p)u_ = _au_a≤ v= _3(p), these imply ωa|ΔJ5| _a| J_5| ≤Jv1/8, _Jv^1/8, (30) ha|J5−J4| h_a|J_5-J_4| ≤Jmin(va,r,s). _J (v_a, r, s). Each antenna counterterm has a logarithmic unresolved measure, bounded by drds/(rs)dr\,ds/(rs) times a reduced tree density. Integrating the second line therefore yields vav_a times at most two logarithms. On the support of 1−ha1-h_a, the transport bound instead implies r,s>c0va2r,s>c_0v_a^2, so the last term of Eq. (28) is bounded after antenna integration by a reduced tree density times va(1+|logva|2) v_a(1+| v_a|^2). Both terms are integrable at the remaining four-parton boundary. For the first term, NLO factorization gives cancellation of the leading single-unresolved residues in D at fixed resolved four-parton kinematics, after the correlated azimuthal average. A further bound is needed when the reduced state also becomes unresolved. For two ordered scales κ and κzκ z, the following inequality gives a useful estimate of the integral over their ratio z. ∫01dzmin(1,κ−Azγ)≤′(1+|logκ|). _0^1 dzz (1,Cκ^-Az^γ) (1+| κ|). (31) Here A≥0A≥ 0 and γ>0γ>0. If the remaining integral over z is bounded by Eq. (31), the positive power in Eq. (30) makes the outer logarithms integrable. Establishing this bound for the complete D contribution in all overlapping soft and collinear limits remains necessary. If p approaches two jets while M^5,ap=b M_5,ap=b remains resolved, uau_a stays nonzero and ωa=O(umin2) _a=O(u_ ^2) ensures the required suppression. Finally, K4K_4 has only logarithmic enhancements of its single-unresolved tree behaviour, rendered integrable by ΔJ4=O(va) J_4=O( v_a). These statements use the established factorization properties of the NLO input and apply after this azimuthal average. The Monte Carlo variance still depends on the sampling density. Contact matching uses the same radiation integral at fixed b, with jettiness as the measurement. For t>0t>0, let ℛ>t(2),(b) _>t^(2),P(b) =FΘ(τ3−t)(b), =F_ ( _3-t)^P(b), (32) c(2)(b) c_P^(2)(b) =limt→0[Σ3,sing(2)(b,t)+ℛ>t(2),(b)]. = _t→ 0 [ _3,sing^(2)(b,t)+R_>t^(2),P(b) ]. This relation defines the contact separately at each resolved three-parton point, integrating only radiation variables. Here Σ3,sing(2) _3,sing^(2) is the NNLO counterpart of the below-t EFT cumulant in Eq. (19), including the coefficient of δ(τ3)δ( _3). Both terms use the same jettiness definition and leading-power Born variables. The positive t defines this limiting characterization; Appendix D develops a finite-difference representation intended to evaluate it without numerical extrapolation. A different projection prescription changes the finite contact, which must therefore be determined with the same P used in Eq. (29). On the mapped four-parton state x, τ3(x)=1−λ _3(x)=1-λ. This explicit coordinate dependence also determines the powers of the vanishing radiation coordinates used in the finite contact construction. The real-state measurement is evaluated separately, preserving the NLO input’s local cancellation without requiring the map to preserve jettiness. Equations (29) and (32) specify the ingredients of the proposed local three-jet NNLO subtraction. Appendix D develops the finite contact representation on the same phase space at fixed b. The numerical checks below test the radiation contribution. B.4 Numerical checks of the fixed-Born radiation contribution Figures 10 and 11 use the complete four-jet NLO input from EERAD3 [5], with nf=5n_f=5, μR=Q _R=Q and 2Ei/Q=(0.8,0.7,0.5)2E_i/Q=(0.8,0.7,0.5). All three Born momenta remain fixed throughout. For J=1−T,CJ=1-T,C, define ℱJ(2)(b,η)=FJ(b,η)B¯3(b),F_J^(2)(b;η)= F_J^P(b;η) B_3(b), (33) Here B¯3(b)=∑π∈S3B3(bπ) B_3(b)= _π∈ S_3B_3(b_π) is the density on the energy-ordered Born region, with bπb_π running over all six permutations of the Born momenta; powers of αs/(2π) _s/(2π) are stripped. The five-parton weights include both correlated azimuths before statistical accumulation. For the cut scan, the minimum pair invariant ymin=mini<jsij/Q2y_ = _i<js_ij/Q^2 must exceed η at both azimuths; this common condition multiplies the entire paired weight. We also evaluate η=0η=0 directly. The four-parton K4ΔJ4K_4 J_4 contribution is included for every cut, and no auxiliary jettiness threshold is imposed. For Fig. 10, the radiation, Born-projected and paired densities are log⟨RaJ5−SaJ4⟩ϕ/B¯3J_ R_aJ_5-S_aJ_4 _φ/ B_3, log⟨Ra−Sa⟩ϕJ3(b)/B¯3J_ R_a-S_a _φJ_3(b)/ B_3 and log⟨WaJ⟩ϕ/B¯3J_ W_a^J _φ/ B_3, respectively. The brackets denote the two-azimuth average defined above Eq. (29); the paired density is evaluated directly from the complete weight. The paths are generated by the common inverse in Eq. (24). The upper row takes δ4=1−λ=r=ϵ _4=1-λ=r=ε at fixed s and four-parton decay angle θ4 _4. The lower row takes δ4=ϵ _4=ε, 1−cosθ4=Aϵ1- _4=Aε and r+s=Bϵr+s=Bε at fixed r/(r+s)r/(r+s). The columns use cyclic Born assignments (1,2,3)(1,2,3), (2,3,1)(2,3,1) and (3,1,2)(3,1,2), with outer pair (3,4)(3,4) and antenna labels (i,j,k)=(2,1,3)(i,j,k)=(2,1,3). From left to right, (cosθ4,s)=(−0.35,0.22),(0.20,0.35),(0.55,0.47)( _4,s)=(-0.35,0.22),(0.20,0.35),(0.55,0.47) above, and (A,B,r/(r+s))=(8,1,0.3),(12,1.3,0.4),(16,1.6,0.5)(A,B,r/(r+s))=(8,1,0.3),(12,1.3,0.4),(16,1.6,0.5) below. Both rows use (ϕ4,ϕ)=(0.73,0.419),(1.23,1.119),(1.73,1.819)( _4,φ)=(0.73,0.419),(1.23,1.119),(1.73,1.819). The logarithmic Jacobian logJ_ is the product of the four-parton and antenna Jacobians, multiplied by δ4r _4r in the upper row and δ4(1−cosθ4)(r+s)2 _4(1- _4)(r+s)^2 in the lower row. No technical cut is applied. The C-parameter gives the same limiting behaviour, and separate single-unresolved scans at resolved four-parton kinematics test the NLO input’s real–counterterm cancellation. Figure 10 shows three double-collinear and three double-soft paths. The locally subtracted radiation density and its Born-projected counterpart approach the same limit, while their paired difference decreases linearly with ϵε down to 10−910^-9. These tests support local cancellation along the sampled paths at this nontrivial fixed Born point. Figure 10: Local cancellation in the three-jet NNLO radiation contribution at fixed Born kinematics 2Ei/Q=(0.8,0.7,0.5)2E_i/Q=(0.8,0.7,0.5), for thrust J=1−TJ=1-T. The upper and lower rows approach three double-collinear and three double-soft limits, respectively. The curves show the magnitudes of the locally subtracted radiation density (blue), its Born-projected counterpart (orange), and their paired difference (green), including the four-jet NLO counterterms. All densities include the logarithmic phase-space Jacobian and are normalized to the tree-level Born density. The paired difference is evaluated before integration and vanishes along all six paths. No technical cut is applied. The paths, normalization and correlated azimuthal average are specified in the text. These checks concern the radiation contribution and do not include the NNLO contact. Figure 11: NNLO radiation integrals at the fixed three-parton Born point of Fig. 10. The ordinate ℱJ(2)F_J^(2) is the radiation integral normalized to the tree-level Born density, with the coupling stripped as in Eq. (33). The upper panels show the paired five-parton integral (blue), the four-parton term containing the finite one-loop and integrated-counterterm weight K4K_4 (orange), and their sum (green). In each sector a, the paired integrand is Ra[J5(p)−J3(b)]−Sa[J4(m^ap)−J3(b)]R_a[J_5(p)-J_3(b)]-S_a[J_4( m_ap)-J_3(b)], where RaR_a is the real sector weight and SaS_a its NLO counterterm, evaluated on the same five-parton state p; b is the common Born point. The middle panels divide the sum by its mean over 10−11≤η≤10−710^-11≤η≤ 10^-7, accounting for correlations. The first two rows share aligned axes, with the directly evaluated η=0η=0 point shown separately. The lower panels show the cumulative mean at η=0η=0 over 160 independent batches. Bands and error bars denote one standard error, including the independent four-parton uncertainty in the sum. The parameter η is a technical cut for numerical-stability checks. These results evaluate the locally subtracted radiation contribution without a small slicing parameter; the NNLO contact is a separate contribution. Figure 11 integrates the same paired construction for thrust and the C-parameter. The upper panels separate the paired five-parton term from the finite four-parton contribution and show their sum. The small-cut results agree with the direct η=0η=0 calculation, and the lower panels show the cumulative estimate settling as independent batches are added. These results demonstrate numerical integration of the NNLO radiation contribution at fixed Born momenta without a small slicing parameter for the measurements tested here. Figure 11 integrates the complete paired weight with importance sampling. To reduce the cost of high-precision evaluation, we use ⟨wf⟩+⟨wh−wf⟩ w_f + w_h-w_f , where wfw_f and whw_h are the importance-weighted local differences evaluated with double- and higher-precision arithmetic. The two averages use independent samples, while each correction evaluates both weights at the same coordinates. In numerically delicate regions the auxiliary wfw_f is set to zero and the correction includes the full whw_h. This rearrangement preserves the locally subtracted radiation integral. A four-hour run on the same M1 laptop produced 160 independently randomized batches, each containing 8192 fast weights and 64 high-precision corrections. Zero-weight unselected charts remain in the sample count. The finite four-parton contribution uses a separate sample of 320 batches of 64 points. The direct zero-cut results are ℱ1−T(2)=39.32±0.89F_1-T^(2)=39.32± 0.89 and ℱC(2)=91.61±2.01F_C^(2)=91.61± 2.01, with relative uncertainties of 2.25%2.25\% and 2.19%2.19\%. The paired five-parton contributions are −0.19±0.66-0.19± 0.66 and −4.78±1.47-4.78± 1.47, respectively; the corresponding four-parton contributions are 39.51±0.5939.51± 0.59 and 96.39±1.3796.39± 1.37. The middle panels use the arithmetic mean over the five cuts 10−11≤η≤10−710^-11≤η≤ 10^-7 as their reference, accounting for all common-sample correlations. The cumulative estimates include the same independent four-parton contribution at every batch count. These integrals contain the measurement differences Jm−J3(b)J_m-J_3(b) and test the locally subtracted radiation contribution. The proposed contact extraction uses the finite-difference weights developed in Appendix D on the same fixed-Born phase space. Algorithm 1 specifies the estimator used for Fig. 11, including the frozen importance density, the independent fast and correction samples, and the separately sampled four-parton contribution. In the code, M4, m[a], g.h and g.omega denote M^4 M_4, m^a m_a, hah_a and ωa _a, respectively. Both precision levels form the complete local difference before integration. The guard acts only on the auxiliary fast weight, whose omission is restored by the high-minus-fast correction. Rejected charts remain in the sample count, and both azimuths form one observation. The ancillary reproduction instructions specify the sampling transforms, normalization conventions and numerical treatment for an independent reproduction of Figs. 10 and 11. Figure 11 uses J=1−T,CJ=1-T,C. Algorithm 1: Integration of locally subtracted NNLO radiation at fixed Born momenta, with the sampling and precision treatment used in Fig. 11. ⬇ 1 # Fixed-Born NNLO radiation. Q=1; nf=5; mu_R=Q; J=(1-T,C). 2 # b=((.4,0,0,.4),(.35,sqrt(3)/8,0,-.275), 3 # (.25,-sqrt(3)/8,0,-.125)); JB=J(b). 4 # R,S: average all 120 labels, then multiply by 5 # 8192*pi**7/Bbar3; K4: all 24, 2048*pi**5/Bbar3. 6 # Bbar3 is the SUM over six Born labels, not a mean. 7 eta = [10.**(-k) for k in range(3,13)] + [0.] 8 proposal = load_frozen("production-proposal.json") 9 # q(u)=.2+.8*q_tree(u); full support, frozen 32 boxes. 10 11 def draw_four(v): # Seven coordinates; Appendix A. 12 h,perm,delta,c4,phi4,pdelta,pc = outer_sample(v) 13 # Six pairs and six Born assignments, uniform. 14 # delta=1-lambda; densities include both mixture 15 # branches; explicit transforms in companion file. 16 x = G4(b,h,perm,1-delta,c4,phi4) 17 if not selected_minimum_pair(x,h): 18 return x,0,delta # Keep zero in sample count. 19 w4 = 36*b[perm[0]].energy*(1-delta)**3 20 w4 /= 16*pi**2*pdelta*pc 21 return x,w4,delta 22 23 def packet(u,precision): # No outer q(u) division yet. 24 a = int(30*u[0]) # All (i,j,k), i<k, distinct. 25 x,w4,delta = draw_four(u[1:8]) 26 if w4 == 0: return zeros((11,2)) 27 r,s,phi,qant = balanced_antenna(u[8:13],delta**2) 28 # qant: full density in dr ds dphi. 20% uniform, 29 # 80% three-corner mixture; no infrared floor. 30 blocks = []; ymin = infinity 31 for angle in (phi,phi+pi/2): 32 p,jac = inverse_antenna(x,a,r,s,angle) 33 # 60-digit kinematics; jac=(I+K)^2/(32*pi^3). 34 ymin = min(ymin,min_invariant(p)) 35 if precision == "fast" and ymin<1e-8: 36 return zeros((11,2)) # Check before fast maps. 37 g = geometry(p,precision) # Cache all 30 maps. 38 active = [c for c in CHANNELS if g.h[c]>0] 39 needed = union(active,[a]) 40 if precision == "fast" and min( 41 g.reduced_min[c] for c in needed)<1e-8: 42 return zeros((11,2)) # L only; H is included in the correction. 43 if precision == "high": assert_close(M4(m[a](p)),b) 44 blocks.append((p,g,jac,active,needed)) 45 values = [] 46 for p,g,jac,active,needed in blocks: 47 R,S = native_kernels(p,needed,precision) 48 # Persistent processes; full label averages; 49 # exact groups sharing a map, sign R-S. 50 D = accurate_sum([R]+[-g.h[c]*S[c] for c in active]) 51 W = accurate_sum([g.omega[a]*D*(J(p)-JB), 52 g.h[a]*S[a]*(J(p)-J(x)), 53 -(1-g.h[a])*S[a]*(J(x)-JB)]) 54 # High: quad kernels and 50-digit sums; if the 55 # packet error bound exceeds 1e-18*abs(W), 56 # recompute the unaggregated packet at SAME p. 57 values.append(w4*30*jac*W/(2*qant)) 58 value = accurate_sum(values) # Average both azimuths. 59 return (ymin>array(eta))[:,None]*value[None,:] 60 61 def design(batch,seed,power): 62 sd = SeedSequence([seed,batch]).generate_state(1)[0] 63 raw = scrambled_Sobol(14,sd,2**power,float64) 64 raw = clip_to_open_unit_cube(raw) 65 return proposal.sample_many(raw[:,:13],raw[:,13]) 66 # Returns (u,q), where q is the FULL mixture density. 67 68 L=[]; C=[] # Independent fast and correction batch means. 69 for batch in range(160): # Complete batches 0,…,159. 70 ul,ql = design(batch,2609179501,13) # 8192 points. 71 uc,qc = design(batch,2609179502,6) # 64 points. 72 low=[]; correction=[] 73 for u,q in zip(ul,ql): 74 low.append(packet(u,"fast")/q) 75 for u,q in zip(uc,qc): 76 l = packet(u,"fast")/q 77 h = packet(u,"high")/q # SAME u as its l. 78 correction.append(h-l) # Includes every guard-zero. 79 # Retry failures at same coordinates; never drop them. 80 L.append(mean(low,axis=0)); C.append(mean(correction,axis=0)) 81 save_complete_batch(L[-1],C[-1]) 82 # The 4h run completed 128 primary + 32 extension batches. 83 # Its incomplete final batch is saved, but not plotted. 84 # Independent finite-four data, already importance weighted. 85 # Each point integrates w4*K4(x)*(J(x)-JB), with its own 86 # proposal weight. Reuse the 320 saved means (64 points each). 87 K = load_independent_K4_means() 88 km = mean(K,axis=0); VK = cov_rows(K)/320 89 paired = mean(L,axis=0)+mean(C,axis=0) # Two means. 90 F = paired+km # Same K4 at EVERY eta, including eta=0. 91 # cov_rows uses rowvar=False, ddof=1; flatten keeps batches. 92 VP = cov_rows(flatten_batches(L))/160 93 VP += cov_rows(flatten_batches(C))/160 94 VF = VP+tile(VK,(11,11)) # Full 22 x 22 covariance. 95 SE = sqrt(diag(VF)).reshape(11,2) 96 # Ratio: A=mean(F[1e-11<=eta<=1e-7]), independently for 97 # each J; use grad(F_i/A) @ VF @ grad(F_i/A) for variance. 98 # Running eta=0: replace 160 by each complete prefix size; 99 # use the SAME km,VK. Last point equals F[-1],SE[-1]. Appendix C Fixed-Born NNLO radiation for arbitrary jet multiplicity The radiation construction in Appendix B extends to any fixed finite N≥2N≥ 2 for massless final-state jets. The outer inverse already holds for general N, while the antenna inverse acts on three partons and leaves all spectators unchanged. We give the phase-space measure at fixed Born momenta and combine each real contribution with its counterterm before integration, assuming an (N+1)(N+1)-jet NLO antenna input with the map in Eq. (23). The integrability argument applies on compact resolved Born regions and to the energy-flow Lipschitz measurements specified in Appendix A. Together with the required Born contact, these ingredients specify a proposal for extending the local subtraction to higher jet multiplicities. Let RN+2R_N+2 be the tree contribution, SaS_a the NLO counterterms grouped by their maps m^a:ΦN+2→ΦN+1 m_a: _N+2→ _N+1, and KN+1K_N+1 the finite one-loop plus integrated-counterterm contribution. For an infrared-safe (N+1)(N+1)-jet measurement fmf_m on m partons, the NLO correction is ℐN+1(1)[f] _N+1^(1)[f] =∫dΦN+2[RN+2fN+2 = _N+2 [R_N+2f_N+2 (34) −∑aSafN+1(m^aΦN+2)] - _aS_af_N+1( m_a _N+2) ] +∫dΦN+1KN+1fN+1. + _N+1\,K_N+1f_N+1. For e+e−→Ne^+e^-→ N jets, these terms are of order αsN _s^N. The NLO kernels, maps and integrated counterterms are recycled [71]; their outer phase-space integrations are replaced below. This requires access to the differential NLO contributions rather than previously integrated histograms. C.1 Radiation phase space at fixed Born momenta Fix b=(p¯1,…,p¯N)b=( p_1,…, p_N), with p¯A2=0 p_A^2=0 and ∑Ap¯A=P=(Q,) _A p_A=P=(Q,0), away from soft and collinear boundaries. Order the Born momenta by decreasing energy, breaking ties by a fixed ordering of direction components. This also removes the label ambiguity for the equal-energy two-parton Born state. A chart h chooses a pair (i,j)(i,j) among N+1N+1 labels and a permutation σ assigning the Born momenta to the merged leg and spectators sAs_A. Define e=p¯σ(1)0e= p_σ(1)^0, =¯σ(1)/en= p_σ(1)/e, k=λek=λ e and ξ=(λ,cosθ,ϕ0)ξ=(λ, θ, _0). The point x=GN+1,h(b,ξ)x=G_N+1,h(b,ξ) is xsA x_s_A =λp¯σ(A),A=2,…,N, =λ p_σ(A), A=2,…,N, (35) Kij K_ij =(1−λ)P+λp¯σ(1), =(1-λ)P+λ p_σ(1), Kij2 K_ij^2 =(1−λ)[Q2(1−λ)+2λQe], =(1-λ) [Q^2(1-λ)+2λ Qe ], xi0 x_i^0 =12(Kij0+kcosθ),xj=Kij−xi, = 12(K_ij^0+k θ), x_j=K_ij-x_i, i _i =12[(k+Kij0cosθ)+Kij2sinθ⟂(ϕ0)]. = 12 [(k+K_ij^0 θ)n+ K_ij^2 θ\,e_ ( _0) ]. Here 0<λ<10<λ<1, and the transverse basis and angular ranges are those of Eq. (16b). These momenta are massless and sum to P for every N. The projection M^N+1 M_N+1 merges the minimum-deficit pair as in Eq. (15) and orders the result. The selector χh _h selects the chart when this pair is (i,j)(i,j), giving M^N+1x=b M_N+1x=b. The (N+12)N! N+12N! charts yield dμN+1b _N+1^b =∑hχhQeλ2N−332π3dλdcosθdϕ0, = _h _h Qeλ^2N-332π^3dλ\,d θ\,d _0, (36) dΦN+1 _N+1 =dΦN(b)dμN+1b. =d _N(b)\,d _N+1^b. The Born measure is restricted to the stated ordering. Summing assignments restores labelled phase space without adding a physical symmetry factor. The minimum deficit is at most the energy of the two least energetic partons, so λ≥(N−1)/(N+1)>0λ≥(N-1)/(N+1)>0 on selected charts. Thus the rescaling introduces no zero denominator for any N≥2N≥ 2. For channel a, write x=(I,K,q1,…,qN−1)x=(I,K,q_1,…,q_N-1), where I,KI,K are the massless antenna parents and the qAq_A are unchanged spectators. Generate their three daughters using Eqs. (22)–(23), with the rotations and boost specified there. This defines p=Gaant(x,r,s,ϕ)∈ΦN+2p=G_a^ant(x;r,s,φ)∈ _N+2. The azimuthal reference uses the spectator with the largest transverse component, with a nonparallel Cartesian axis when all such components vanish. The inverse and factorized measure obey m^aGaant(x,r,s,ϕ) m_aG_a^ant(x;r,s,φ) =x, =x, (37) M^N+2,a M_N+2,a =M^N+1∘m^a, = M_N+1 m_a, M^N+2,ap M_N+2,ap =M^N+1x=b, = M_N+1x=b, dΦN+2 _N+2 =dΦN(b)dμN+1bsIK32π3drdsdϕ. =d _N(b)\,d _N+1^b\, s_IK32π^3dr\,ds\,dφ. The antenna factorization is independent of the spectator multiplicity [71]. On each selected chart, the combined Jacobian is Qeλ2N−3sIK/(1024π6)Qeλ^2N-3s_IK/(1024π^6). There are six radiation coordinates at fixed b, irrespective of N. Reapplying the NLO subtraction map recovers exactly the same x used in its counterterm, and both states project to the same b. C.2 Local pairing and integrability Use minimized geometric N-jettiness as defined in Eq. (10). For N+1N+1 massless partons its minimizing partition can be chosen as one pair and N−1N-1 singletons. Hence va=τN(x)=1−λv_a= _N(x)=1-λ on a selected chart. Define ℓa _a, hah_a, uau_a and ωa _a by Eqs. (25)–(27), using this vav_a and the N+2N+2 real labels. Include all (N+2)(N+12)(N+2) N+12 geometric channels (i,j,k)(i,j,k) with distinct labels and i<ki<k, setting Sa=0S_a=0 for absent counterterms. As in Appendix B, use the same full (N+2)!(N+2)! label average for RN+2,SaR_N+2,S_a and an (N+1)!(N+1)! average for KN+1K_N+1, preserving the input’s flavour sums and symmetry factors. Set D D =RN+2−∑chcSc,Ra=haSa+ωaD, =R_N+2- _ch_cS_c, R_a=h_aS_a+ _aD, (38) ΔJN+2 J_N+2 =J(p)−JN(b),ΔJN+1=J(x)−JN(b), =J(p)-J_N(b), J_N+1=J(x)-J_N(b), WaJ W_a^J =RaΔJN+2−SaΔJN+1 =R_a J_N+2-S_a J_N+1 =ωaDΔJN+2+haSa[J(p)−J(x)] = _aD J_N+2+h_aS_a[J(p)-J(x)] −(1−ha)SaΔJN+1. -(1-h_a)S_a J_N+1. Form this complete weight on the common phase space before integration. For scalar antennae, average it at ϕφ and ϕ+π/2φ+π/2 to cancel the leading collinear spin harmonics as in Appendix B. The two states share Born coordinates but need not share jettiness. With NP_N denoting this projection and partition prescription, the radiation integral at fixed b and the NNLO correction are FJN(b) F_J^P_N(b) =∑a∫dμN+1bdΦant,a⟨WaJ⟩ϕ = _a _N+1^b\,d _ant,a W_a^J _φ (39) +∫dμN+1bKN+1ΔJN+1, + _N+1^b\,K_N+1 J_N+1, ΣN(2)[J] _N^(2)[J] =∫dΦN(b)[FJN(b)+cN(2)(b)JN(b)]. = _N(b) [F_J^P_N(b)+c_P_N^(2)(b)J_N(b) ]. The identity ∑aRa=RN+2 _aR_a=R_N+2 restores the physical real and one-loop terms of Eq. (34) after integration over b. The remaining projected term defines the Born contact. For J=1J=1, the radiation weight vanishes before integration. To establish integrability for general N, first note that the unchanged spectators contribute nothing to the transport distance ℓa _a. The bounds dE(x,b)≤vad_E(x,b) v_a, ℓa≤min(r,s) _a (r,s) and ℓa≤ua1/4 _a _a^1/4 therefore hold as in Appendix B. Set v=τN(p)v= _N(p). Every partition of N+2N+2 momenta into N clusters contains an internal pair. The triangle inequality bounds that pair’s deficit by its cluster’s deficit, giving umin≤vu_ ≤ v for every N. Since va≤v+ℓav_a≤ v+ _a, these facts imply dE(p,b) d_E(p,b) ≤(v+ua1/8), ( v+u_a^1/8), (40) ωa|ΔJN+2| _a| J_N+2| ≤Jv1/8, _Jv^1/8, ha|J(p)−J(x)| h_a|J(p)-J(x)| ≤Jmin(va,r,s). _J (v_a, r, s). For the second line, use ωa≤min[1,(v/ua)2] _a≤ [1,(v/u_a)^2] and split at ua=vu_a=v. The constants may depend on N and the compact Born region, but the positive exponents do not. The last line makes the explicit antenna difference integrable with its logarithmic measure. On the support of 1−ha1-h_a, r,s>c0va2r,s>c_0v_a^2, leaving a reduced tree density times va(1+|logva|2) v_a(1+| v_a|^2) after antenna integration. The remaining (N+1)(N+1)-parton single-unresolved integration is finite. For D, the NLO input cancels single-unresolved residues, including collinear spin correlations after averaging over ϕφ and ϕ+π/2φ+π/2. When the reduced state also becomes unresolved, controlling D requires the uniform bound discussed in Appendix B. If p approaches fewer than N jets while its projection remains at the resolved b, uau_a cannot vanish. Otherwise ℓa→0 _a→ 0 and va→0v_a→ 0 would force b to the same lower-multiplicity boundary. Hence ωa=O(umin2) _a=O(u_ ^2) suppresses such limits. Finally, ΔJN+1=O(va) J_N+1=O( v_a) renders the logarithmically enhanced single-unresolved behaviour of KN+1K_N+1 integrable. The required contact is defined at each Born point by cN(2)(b)=limt→0[ΣN,sing(2)(b,t)+FΘ(τN−t)N(b)].c_P_N^(2)(b)= _t→ 0 [ _N,sing^(2)(b,t)+F_ ( _N-t)^P_N(b) ]. (41) The EFT cumulant uses the same leading-power Born variables and includes the coefficient of δ(τN)δ( _N), with its process-dependent two-loop hard contribution. Only radiation coordinates are integrated in F. The availability of this EFT input and of the (N+1)(N+1)-jet NLO calculation is separate from the explicit kinematic construction above. C.3 Recursion and the four-jet example The inverse maps admit a recursive form. Starting with x1=GN+1,h(b,ξ)x_1=G_N+1,h(b,ξ) and ℳ^1=M^N+1 M_1= M_N+1, let aka_k choose an antenna insertion into xk∈ΦN+kx_k∈ _N+k, with radiation variables ξk=(rk,sk,ϕk) _k=(r_k,s_k, _k). On each chosen sequence of channels, xk+1 x_k+1 =Gakant(xk,ξk), =G_a_k^ant(x_k; _k), (42) ℳ^k+1 M_k+1 =ℳ^k∘m^ak,ℳ^kxk=b, = M_k m_a_k, M_kx_k=b, dμN+kb _N+k^b =dμN+1b∏j=1k−1dΦant,aj. =d _N+1^b _j=1^k-1d _ant,a_j. The inverse identity m^akGakant=id m_a_kG_a_k^ant=id proves preservation of b by induction; phase-space factorization proves the product measure. Each insertion adds three coordinates. At NNLO only x1,x2x_1,x_2 are required, and Eq. (38) gives their locally subtracted weight. Further iteration generates additional radiation while keeping b fixed; higher perturbative orders also require the corresponding radiation subtraction and EFT singular distribution. For N=3N=3, these equations reproduce Appendix B. For N=4N=4, they give 240 outer charts, 60 antenna channels and the Jacobian Qeλ5sIK/(1024π6)Qeλ^5s_IK/(1024π^6), with no planarity restriction on the four-parton Born state. A five-jet NLO input has been implemented and checked in the recent NNLOJET four-jet NNLO calculation [41], using OpenLoops2 for the six-parton tree and five-parton one-loop amplitudes. Such an input could support a future application beyond three jets. The radiation subtraction at fixed Born momenta and contact matching at these higher multiplicities have not been evaluated numerically here. Appendix D Local contact extraction at NLO and NNLO This appendix derives the finite radiation integrals that, together with the EFT δ(τN)δ( _N) coefficient, directly determine the Born contact. The integrals use radiation configurations generated at fixed Born momenta and cover the full phase space without a small slicing parameter. We first illustrate the finite-rescaling construction in one variable, then give the explicit NLO result used in Sec. VII.2. We next describe the recursive two-jet NNLO contact used in Sec. VII.3 and a proposed three-jet NNLO extension using the radiation input of Appendix B. D.1 A one-dimensional example Consider a radiation density g(x)/xg(x)/x on 0<x<10<x<1, where x measures the distance from the singular boundary x=0x=0. Let g be regular away from x=0x=0, with g(x)=g0+O(xα)g(x)=g_0+O(x^α), α>0α>0 and g0=g(0)g_0=g(0). The integral above a small cut t contains −g0lnt-g_0 t. Its singular counterpart below the cut is H+g0lntH+g_0 t, where H is a known constant in the singular cumulant. Contact matching therefore gives c c =limt→0[H+g0lnt+∫t1dxg(x)]=H+A, = _t→ 0 [H+g_0 t+ _t^1 dxxg(x) ]=H+A, (43) A A =∫01dx[g(x)−g0]. = _0^1 dxx[g(x)-g_0]. The logarithms cancel, leaving the finite radiation contribution A. For example, g(x)=1+xg(x)=1+x gives A=1A=1 and c=H+1c=H+1. We can obtain the same A using only interior values of g. Choose the fixed rescaling x↦x/2x x/2 and define Δg(x)=g(x)−g(x/2) g(x)=g(x)-g(x/2) and ℓ=ln2 = 2. Since Δg=O(xα) g=O(x^α), its integrals with 1/x1/x and lnx/x x/x are finite. Changing variables from x to x/2x/2 in the shifted terms gives ∫01dxΔg(x) _0^1 dxx g(x) =∫1/21dxg(x)−g0ℓ, = _1/2^1 dxxg(x)-g_0 , (44) ∫01dxlnxΔg(x) _0^1 dxx x\, g(x) =−ℓA+∫1/21dxln(2x)g(x)−g0ℓ22. =- A+ _1/2^1 dxx (2x)g(x)- g_0 ^22. Eliminating g0g_0 between these two identities yields A= A= ∫01dx(12−lnxℓ)[g(x)−g(x/2)] _0^1 dxx ( 12- x )[g(x)-g(x/2)] (45) +∫1/21dx(12+lnxℓ)g(x). + _1/2^1 dxx ( 12+ x )g(x). Both integrals converge separately. In the first, the common value g0g_0 cancels within each sampled difference. The second restores the finite interval left by the change of variables. It is part of the answer and cannot be dropped. Substituting g(x)=1+xg(x)=1+x again gives A=1A=1. The integration covers the entire interval 0<x<10<x<1; 1/21/2 is a rescaling factor, not an infrared cut. Equation (45) thus determines the contact from finite weights without evaluating g at x=0x=0 or taking a small numerical limit. Alternatively, identify A as the constant term of an auxiliary radiation integral. This also provides a direct derivation of the finite weights in Eq. (45), which is the route we extend to NLO and NNLO below. For Reζ>0Reζ>0, define (ζ)=∫01dxg(x)xζ=g0ζ+A+O(ζ).Z(ζ)= _0^1 dxxg(x)x^ζ= g_0ζ+A+O(ζ). (46a) Thus A=[ζ0]A=[ζ^0]Z. To evaluate this constant using only interior values of g, change variables in the term containing g(x/2)g(x/2), including the interval [1/2,1][1/2,1]. This gives (ζ)=11−2ζ[ (ζ)= 11-2^ζ [ ∫01dxxζΔg(x) _0^1 dxxx^ζ g(x) (46b) −2ζ∫1/21dxζg(x)]. -2^ζ _1/2^1 dxxx^ζg(x) ]. Since (1−2ζ)−1=−1/(ℓζ)+1/2+O(ζ)(1-2^ζ)^-1=-1/( ζ)+1/2+O(ζ), the constant terms of the two weights are [ζ0]xζ1−2ζ [ζ^0] x^ζ1-2^ζ =12−lnxℓ, = 12- x , (46c) [ζ0]−2ζxζ1−2ζ [ζ^0] -2^ζx^ζ1-2^ζ =12+lnxℓ. = 12+ x . Substituting these weights into Eq. (46b) gives exactly Eq. (45). The logarithmically weighted difference is integrable, so this constant term can be taken inside the integrals. The strategy is therefore to identify the finite contribution through the auxiliary integral and then convert it into convergent integrals by finite rescaling. The numerical calculation uses these finite weights, not a small numerical ζ. The auxiliary exponent is not a dimensional regulator. At higher orders, repeated differences remove the additional logarithms of the unresolved radiation coordinates. D.2 Singular-distribution convention and contact For the physical calculation, fix a resolved Born configuration b=(p¯1,p¯2,p¯3)b=( p_1, p_2, p_3). Radiation variables generate a four-parton state Φ4 _4 or a five-parton state Φ5 _5 while these three Born momenta remain fixed. The minimized geometric three-jettiness is τ3(Φm) _3( _m) on either state, as defined in Eq. (10). As in the one-dimensional example, the contact is the sum of the coefficient of δ(τ3)δ( _3) in the EFT singular distribution and a finite radiation contribution. At NLO (k=1k=1) or NNLO (k=2k=2), the singular cumulant differential in b is Σ3,sing(k)(b,t)=Hk(b)+∑j=02k−1ak,j(b)j+1lnj+1t. _3,sing^(k)(b,t)=H_k(b)+ _j=0^2k-1 a_k,j(b)j+1 ^j+1t. (47) Here t is the upper limit on three-jettiness, ak,j(b)a_k,j(b) are the coefficients of its singular logarithms, and Hk(b)H_k(b) is the constant term, including the finite virtual contribution. These are EFT inputs in the same Born variables, resolution definition, scale and normalization as the radiation calculation. Logarithms are referred to the dimensionless interval [0,1][0,1]. As in Eq. (43), HkH_k must still be combined with the finite radiation integral to obtain the P2B contact. The required EFT input includes the finite hard, jet and soft contributions. Numerical methods have been developed for two-loop soft and beam functions [28, 22, 20, 21, 19], while jet [29, 6] and N-jettiness beam functions [54, 7, 83, 55, 84, 56] are known analytically through three loops. These ingredients must be combined with the hard function and conventions for the process under consideration. We use the following four-dimensional radiation inputs. At NLO, R4(Φ4)R_4( _4) is the four-parton tree coefficient. At NNLO, channel a has a five-parton real contribution Ra(Φ5)R_a( _5) and a local counterterm Sa(Φ5)S_a( _5) whose measurement is evaluated on Φ4=m^aΦ5 _4= m_a _5. The assigned real contributions obey ∑aRa=R5 _aR_a=R_5, where R5R_5 is the complete five-parton tree coefficient; their definition is Eq. (27). The remaining input K4(Φ4)K_4( _4) is the renormalized four-parton one-loop coefficient plus its integrated NLO counterterms, finite at resolved four-parton kinematics. The measure dμ4bd _4^b integrates only the radiation coordinates ξ of the inverse map G4,h(b,ξ)G_4,h(b,ξ), with h labelling its branches. It includes the Jacobian and selected branch sum, so dΦ4=dΦ3(b)dμ4bd _4=d _3(b)\,d _4^b. In channel a, dμ5,ab=dμ4bdΦant,ad _5,a^b=d _4^bd _ant,a includes the extra-emission measure of Eq. (24). Both generated states have the same prescribed Born point. The brackets ⟨X⟩ϕ=[X(ϕ)+X(ϕ+π/2)]/2 X _φ=[X(φ)+X(φ+π/2)]/2 use the correlated antenna-azimuth pair needed for local spin cancellation in Appendix B. We now insert the auxiliary measurement τ3ζ(Φm) _3^ζ( _m), just as xζx^ζ was inserted in Eq. (46a). The two radiation integrals are 1(ζ,b) _1(ζ;b) =∫dμ4bR4τ3ζ(Φ4), = _4^b\,R_4 _3^ζ( _4), (48) 2(ζ,b) _2(ζ;b) =∑a∫dμ5,ab⟨Raτ3ζ(Φ5)−Saτ3ζ(Φ4)⟩ϕ = _a _5,a^b R_a _3^ζ( _5)-S_a _3^ζ( _4) _φ +∫dμ4bK4τ3ζ(Φ4), + _4^b\,K_4 _3^ζ( _4), The real and counterterm measurements in the second line are evaluated on their respective states, within the same local pair. Expanding in the formal exponent ζ separates the poles in ζ generated by the singular radiation integrals from the constant radiation term k(b)=[ζ0]k(ζ,b)A_k(b)=[ζ^0]Z_k(ζ;b). The contacts are cM^(b) c_ M(b) =H1(b)+1(b), =H_1(b)+A_1(b), (49) c(2)(b) c_P^(2)(b) =H2(b)+2(b). =H_2(b)+A_2(b). The subscripts specify the projection prescription. Here M M is the NLO map of Appendix A, while P includes the outer map, NLO subtraction maps and real-channel assignment of Appendix B. The plus sign has the same origin as in Eq. (43). The radiation integral above t cancels the logarithms of Eq. (47) and leaves the constant kA_k. The following construction represents these constants through finite differences and compensation intervals, generalizing Eq. (45). The exact radiation kernels include all nonsingular terms. No resolved event is approximated by its small-jettiness limit, and neither a small numerical ζ nor an additional dimensionally regulated integral is required. The explicit NLO result is evaluated in the main text. The NNLO construction below includes the two-jet application and the proposed three-jet extension. D.3 Explicit NLO construction Set Q=1Q=1, with momenta expressed in units of Q, and use the inverse in Eqs. (16)–(17). Its variable λ rescales the Born spectators, while θ and ϕφ are the polar and azimuthal angles of the two-body decay. Recall that χh=1 _h=1 when h=(i,j)h=(i,j) is the pair minimizing dij=Ei+Ej−|i+j|d_ij=E_i+E_j-|p_i+p_j|, and χh=0 _h=0 otherwise. Define δ=1−λδ=1-λ and u=1−|cosθ|u=1-| θ|. The inverse gives dij=δd_ij=δ and implies χh=0 _h=0 whenever u<δu<δ. Thus every four-parton phase-space configuration with χh=1 _h=1 can be parametrized, for σ=±1σ=± 1, by δ=uv,cosθ=σ(1−u),0<u,v<1.δ=uv, θ=σ(1-u), 0<u,v<1. (50) Here τ3(Φ4)=dij=uv _3( _4)=d_ij=uv for χh=1 _h=1. To verify χh=0 _h=0 for u<δu<δ, use the momenta in Eq. (16b). Here e1=p¯10e_1= p_1^0, k=λe1k=λ e_1 and n is the direction of the Born leg replaced by the pair. If u<δu<δ, the lower-energy daughter q has Eq=(δ+ku)/2<δE_q=(δ+ku)/2<δ and ⋅=[(k+δ)u−δ]/2<0q·n=[(k+δ)u-δ]/2<0. The spectator momenta sum to −k-kn, so at least one spectator r satisfies ⋅r>0q·p_r>0. Its pair with q then obeys dqr=Eq+Er−|+r|≤Eq<dijd_qr=E_q+E_r-|q+p_r|≤ E_q<d_ij. Hence (i,j)(i,j) cannot minimize d, and χh=0 _h=0 there. Equation (50) therefore covers all configurations with χh=1 _h=1. Include the exact matrix element, radiation measure at fixed b and branch selection in g, defined by R4dμ4b R_4d _4^b =∑h,σdudvuvdϕgh,σ(u,v,ϕ,b), = _h,σ du\,dvuv\,dφ\,g_h,σ(u,v,φ;b), (51) gh,σ g_h,σ =χhe1(1−uv)3u2v32π3R4. = _h e_1(1-uv)^3u^2v32π^3R_4. The branch label h also includes the assignments of the Born momenta to parton labels. Flavour and identical-particle factors are inherited from R4R_4. Use the common average over ϕ+nπ/2φ+nπ/2, n=0,1,2,3n=0,1,2,3, in g and its limits. Suppressing the branch and angular arguments, define gs(v)=g(0,v)g_s(v)=g(0,v), gc(u)=g(u,0)g_c(u)=g(u,0) and g0=g(0,0)g_0=g(0,0). The finite radiation contribution is the two-variable generalization of the second line of Eq. (43). We subtract the soft and collinear limiting functions gsg_s and gcg_c and add back their common value g0g_0, which would otherwise be subtracted twice. Since the resolution variable is uvuv, the boundaries also contribute finite logarithmic terms. The result is 1=∑h,σ∫dϕ _1= _h,σ φ\, \ ∫01dudvuv[g−gs−gc+g0] _0^1 du\,dvuv[g-g_s-g_c+g_0] (52) +∫01dvvlnv[gs(v)−g0] + _0^1 dvv v\,[g_s(v)-g_0] +∫01duulnu[gc(u)−g0]. + _0^1 duu u\,[g_c(u)-g_0] \. The limiting functions follow from the soft and collinear limits of the tree amplitude. The double difference vanishes with a positive power in each unresolved coordinate; changes of the selected pair are confined to shrinking angular intervals near these limits. These properties make all three integrals finite. Equation (52) can therefore be evaluated directly once these limiting functions are known. Equation (52) uses the limiting functions gsg_s, gcg_c and g0g_0, just as Eq. (43) uses g0g_0. To avoid using these limiting functions explicitly, we apply the finite rescaling in Eq. (46b) to 1Z_1 in Eq. (48), using τ3(Φ4)=uv _3( _4)=uv and the density in Eq. (51). Taking the constant term will then give the two-variable counterpart of Eq. (45). The resulting integral uses only interior values of g and is more readily extended to higher orders, where overlapping limits make explicit soft and collinear limiting functions more involved. Choose 0<ρ<10<ρ<1, set L=lnρL= ρ, and define Δug(u,v)=g(u,v)−g(ρu,v) _ug(u,v)=g(u,v)-g(ρ u,v) and Δvg(u,v)=g(u,v)−g(u,ρv) _vg(u,v)=g(u,v)-g(u,ρ v). Applying both differences gives g(u,v)−g(ρu,v)−g(u,ρv)+g(ρu,ρv)g(u,v)-g(ρ u,v)-g(u,ρ v)+g(ρ u,ρ v). Changing variables in one shifted term gives, at Reζ>0Reζ>0, (1−ρ−ζ)∫01duuζg(u,v) (1-ρ^-ζ) _0^1 duuu^ζg(u,v) =∫01duuζΔug(u,v) = _0^1 duuu^ζ _ug(u,v) (53) −ρ−ζ∫ρ1duζg(u,v). -ρ^-ζ _ρ^1 duuu^ζg(u,v). This is the general-ρ version of Eq. (46b). Applying it in both u and v gives the denominator (1−ρ−ζ)2(1-ρ^-ζ)^2. Taking its constant term, as in Eq. (46c), now gives the polynomial weight P1(ℒ) P_1(L) =[ζ0]eℒζ(1−e−Lζ)2 =[ζ^0] e^Lζ(1-e^-Lζ)^2 (54) =ℒ22L2+ℒL+512. = L^22L^2+ LL+ 512. Here ℒL is the logarithm multiplying ζ in the numerator. It is ln(uv) (uv) for the double difference; each compensation contributes a factor ρ−ζρ^-ζ and hence shifts ℒL by −L-L. The second line follows from (1−e−Lζ)−2=1/(L2ζ2)+1/(Lζ)+5/12+O(ζ)(1-e^-Lζ)^-2=1/(L^2ζ^2)+1/(Lζ)+5/12+O(ζ). The complete result is 1=∑h,σ∫dϕ _1= _h,σ φ\, \ ∫01du∫01dvP1(lnuv)ΔuΔvg _0^1 duu _0^1 dvvP_1( uv)\, _u _vg (55) −∫ρ1du∫01dvvP1(lnuv−L)Δvg - _ρ^1 duu _0^1 dvvP_1( uv-L)\, _vg −∫01du∫ρ1dvvP1(lnuv−L)Δug - _0^1 duu _ρ^1 dvvP_1( uv-L)\, _ug +∫ρ1du∫ρ1dvvP1(lnuv−2L)g. + _ρ^1 duu _ρ^1 dvvP_1( uv-2L)\,g \. The first term cancels the singularities at both u=0u=0 and v=0v=0 through four interior evaluations. The next two restore the finite intervals left by each rescaling, and the last restores their overlap. All four terms are required. The exact change of variables in Eq. (53) ensures that this expression gives the same 1A_1 as Eq. (52), now using only interior values of g. It is the two-variable counterpart of Eq. (45). The remainders that vanish with positive powers of the unresolved coordinates make the differences integrable with the displayed logarithmic weights. The parameter ρ can be chosen freely in (0,1)(0,1). It changes the finite integral representation, not the contact, and is not a cutoff; the exact compensations restore the full phase space. To complete the contact, combine this finite radiation integral with the EFT coefficient of δ(τ3)δ( _3). In the same Born-differential normalization, the NLO singular three-jettiness spectrum is dΣ3,sing(1)dτ3(b,τ3) d _3,sing^(1)d _3(b, _3) =H1(b)δ(τ3)+a1,0(b)[1τ3]+ =H_1(b)δ( _3)+a_1,0(b) [ 1 _3 ]_+ (56) +a1,1(b)[lnτ3τ3]+. +a_1,1(b) [ _3 _3 ]_+. The plus distributions use the interval [0,1][0,1], as in Eq. (47). Thus H1(b)H_1(b) is the coefficient multiplying δ(τ3)δ( _3) in the complete EFT prediction, with the hard, jet and soft contributions combined in this convention. The directly constructed NLO contact is cM^(b)=H1(b)+1(b), c_ M(b)=H_1(b)+A_1(b), (57) where 1A_1 is the explicit finite radiation integral in Eq. (52), or the interior expression in Eq. (55). Both terms are densities at the same fixed Born point and use the same jettiness definition, scale and normalization. Algorithm 2 gives the NLO evaluation used in Sec. VII.2. It evaluates all four terms of Eq. (55) at common radiation coordinates, using indicators for the restored intervals. The calculation averages over four azimuthal orientations and uses independent samples at ρ=0.1ρ=0.1 and 0.010.01. Algorithm 2: Local NLO contact evaluation used in Sec. VII.2. The density g is defined in Eq. (51); the weight combines all four terms of Eq. (55) at the same sampled coordinates. The sampling density includes both the continuous radiation density and the probability of the chosen branch. ⬇ INPUT: resolved Born momenta b; the four-parton tree coefficient R4; the complete EFT delta(tau_3) coefficient H1(b). Use the common normalization and Born measure of Appendix D. Choose rho = 0.1 or 0.01, as in the two independent production campaigns. Set L = log(rho) and P1(log_tau) = log_tau^2/(2*L^2) + log_tau/L + 5/12. For each sample in each independently scrambled Sobol batch: Draw (u,v,phi,h,sigma) with 0<u,v<1 and 0<=phi<2*pi. Use the frozen sampling density q, including the branch probability. h includes all six pair choices and six Born assignments; sigma = +/-1 gives the two polar signs. For i,j in 0,1: g[i,j] = 0 For n in 0,1,2,3: us = rho^i*u; vs = rho^j*v; angle = phi+n*pi/2 Generate four momenta at the same b, with 1-lambda = us*vs and cos(theta) = sigma*(1-us). Recompute the selected-pair condition on this state. Evaluate the radiation density g_h,sigma(us,vs,angle;b), including the tree coefficient, Jacobian and selection. g[i,j] += g_h,sigma(us,vs,angle;b)/4 log_tau = log(u*v) w = P1(log_tau)*(g[0,0]-g[1,0]-g[0,1]+g[1,1]) If u >= rho: w -= P1(log_tau-L)*(g[0,0]-g[0,1]) If v >= rho: w -= P1(log_tau-L)*(g[0,0]-g[1,0]) If u >= rho and v >= rho: w += P1(log_tau-2*L)*g[0,0] w /= u*v*q Accumulate the complete weight in the batch mean. Estimate A1(b) and its error from the independent batch means. Return c_M(b) = H1(b) + A1(b). Every scaled evaluation uses sufficient numerical precision; no large weight or small radiation coordinate is discarded. rho rescales the coordinates and leaves the full integration domain in the estimator. For the distributions, also integrate over b with the physical Born measurement and include the covariance with the P2B radiation term. The direct limiting-function expression is evaluated on the same samples for the contact comparison, including the shared-sample covariance. D.4 Reusing NLO subtraction for NNLO contact extraction The contact construction can reuse any complete compatible extra-jet NLO calculation, including its real, virtual and integrated-subtraction contributions. In Fig. 8, the radiation term uses the EERAD3 antenna NLO input, while the two-jet contact uses our NLO P2B construction recursively. The antenna representation below provides an alternative NLO input for the three-jet extension. Both contact constructions use the same finite-part relation; the maps and compensations follow the chosen NLO representation. For the two-jet application, let 2(ζ)Z_2(ζ) denote the three-jet NLO radiation integral with measurement τ2ζ _2^ζ, integrated over the Born orientation and normalized as in Eq. (13). Inserting the three-jet NLO contact H1+1H_1+A_1 cancels the projected real contribution against its counterpart in 1A_1 before integration, giving 2(ζ) _2(ζ) =[ζ10]∫dΦ4R4τ3(Φ4)ζ1τ2(Φ4)ζ =[ _1^0] _4\,R_4 _3( _4) _1 _2( _4)^ζ (58) +∫dΦ3H1(Φ3)τ2(Φ3)ζ, + _3\,H_1( _3) _2( _3)^ζ, 2 _2 =[ζ0]2(ζ),c2=H2+2. =[ζ^0]Z_2(ζ), c_2=H_2+A_2. Here τ3 _3 is the globally minimized three-jettiness of Eq. (10), τ2=1−T _2=1-T, and H1(Φ3)H_1( _3) is the three-jet EFT δ(τ3)δ( _3) coefficient. The two-jet coefficient H2H_2 includes the hard, jet and soft contributions defined in Sec. VI. The inner finite coefficient in ζ1 _1 is taken first at fixed ζ. As at NLO, both exponents are formal devices; the numerical calculation uses finite differences and their compensation integrals in four dimensions. Terms sharing an unresolved limit are combined before integration. Learned control functions then reduce the variance through Eq. (14), leaving this contact formula unchanged. For the three-jet extension, Eq. (49) defines 2(b)=[ζ0]2(ζ,b)A_2(b)=[ζ^0]Z_2(ζ;b). The integration weight includes the four-parton term, the five-parton real–counterterm differences and all finite compensations at fixed Born momenta. A useful reorganization transfers a finite part of an antenna counterterm between the two multiplicities. At a fixed four-parton state x, choose a dimensionless profile ϑa _a such that its antenna integral is finite, and define ΔK4(x)=∑a∫dΦant,a⟨ϑaSa⟩ϕ. K_4(x)= _a _ant,a\, _aS_a _φ. (59) The profile, counterterm and measure use the same inverse map p=Gaant(x,r,s,ϕ)p=G_a^ant(x;r,s,φ). Adding this contribution to the five-parton integrand and subtracting its integral from K4K_4 gives the exact identity 2(ζ,b) _2(ζ;b) =∑a∫dμ5,ab⟨Raτ3ζ(Φ5) = _a _5,a^b R_a _3^ζ( _5) (60) −(1−ϑa)Saτ3ζ(Φ4)⟩ϕ -(1- _a)S_a _3^ζ( _4) _φ +∫dμ4b[K4−ΔK4]τ3ζ(Φ4). + _4^b\,[K_4- K_4] _3^ζ( _4). The transfer and its compensation use the same four-parton state and together cover the full phase space. The actual τ3(Φ5) _3( _5) and τ3(Φ4) _3( _4) are evaluated on their respective states, with both azimuthal orientations included. D.5 Finite differences and exact compensations The finite-rescaling identity used at NLO also applies to the complete expression in Eq. (60). In each radiation region, choose coordinates zi∈(0,1)z_i∈(0,1) that vanish in its soft or collinear limits, and denote the remaining variables and measures by y and dydy. Include the phase-space Jacobians, branch selections and channel weights in gℓg_ . After extracting the logarithmic measure ∏idzi/zi _idz_i/z_i, write f(,y,ζ)=∑ℓgℓ(,y)τ3,ℓ(,y)ζ,τ3,ℓ=∏iziaiℓUℓ.f( z,y;ζ)= _ g_ ( z,y) _3, ( z,y)^ζ, _3, = _iz_i^a_i U_ . (61) Here τ3,ℓ _3, is three-jettiness on the state assigned to contribution ℓ , including the transferred terms. The powers aiℓ≥0a_i ≥ 0 describe how jettiness vanishes; UℓU_ contains its remaining dependence. Its logarithms contribute to the finite coefficient and must be kept. Additional zeros of UℓU_ require further coordinates. Each coordinate region retains its physical map and domain, including the complementary regions of any partition. For each coordinate, put the distinct positive powers aiℓa_i in a list ℰiE_i. Repeat a power p+1p+1 times when its leading coefficient contains logarithms through degree p. For example, a quadratic logarithm requires three successive differences. Choose 0<ρi<10< _i<1 and define if(zi)=f(ρizi) T_if(z_i)=f( _iz_i). The required finite differences and their polynomial weight are i(ζ) _i(ζ) =∏a∈ℰi(1−ρi−aζi)=∑j,Aci,j,Aρi−Aζij, = _a _i(1- _i^-aζ T_i)= _j,Ac_i,j,A _i^-Aζ T_i^j, (62) Di(ζ) D_i(ζ) =∏a∈ℰi(1−ρi−aζ), = _a _i(1- _i^-aζ), P(ℒ) P(L) =[ζ0]eℒζ∏iDi(ζ). =[ζ^0] e^Lζ _iD_i(ζ). Here j counts rescalings, A is the sum of selected powers, and ci,j,Ac_i,j,A is their signed multiplicity. At NLO, ℰu=ℰv=1E_u=E_v=\1\ reproduces P1P_1 in Eq. (54). Coordinates with no positive jettiness power retain the original NLO subtraction and require control of its remaining unresolved limits. Each change of variables produces a finite interval, which is restored through the identity Di∫01dzizif D_i _0^1 dz_iz_if =∫01dziziif = _0^1 dz_iz_iQ_if (63) +∑j≥1,Aci,j,Aρi−Aζ∫ρij1dzizif. + _j≥ 1,Ac_i,j,A _i^-Aζ _ _i^j^1 dz_iz_if. Applying this identity to all rescaled coordinates includes every combination of full and restored intervals. Their overlaps are included, as in the four NLO terms of Eq. (55). The constant coefficient is obtained algebraically, without evaluating a small numerical ζ. One possible parametrization maps every interval to 0<vi<10<v_i<1. Let ki=0k_i=0 denote the full interval and ki>0k_i>0 the restored interval [ρiki,1][ _i^k_i,1], with ki=0,…,|ℰi|k_i=0,…,|E_i|. Here |ℰi||E_i| counts entries including repetitions, and ci=ci,ji,Aic_i=c_i,j_i,A_i. For ki=0k_i=0, sum all terms in iQ_i; for ki>0k_i>0, retain those with ji=kij_i=k_i. The evaluation coordinate and integration factor are (zi,νi)=(viρiji,ci/vi),ki=0,(ρiki(1−vi),−cikilnρi),ki>0.(z_i, _i)= cases(v_i _i^j_i,\ c_i/v_i),&k_i=0,\\ ( _i^k_i(1-v_i),\ -c_ik_i _i),&k_i>0. cases (64) For a given radiation region and interval combination, the resulting weight is W(,y) W_ k( v,y) =∑ji,Ai(∏iνi)∑ℓgℓ(,y) = _\j_i,A_i\ ( _i _i ) _ g_ ( z,y) (65) ×P(lnτ3,ℓ−∑iAilnρi). × P\! ( _3, - _iA_i _i ). At each rescaled point, regenerate the physical momenta, their mapped counterterms and all selection factors. The Born momenta remain fixed. All signed terms and azimuthal averages are formed before integration, with the transfer and its compensation treated consistently. The cancellation of leading logarithms follows from the finite-difference identity (1−ρi−aζi)p+1[ziaζp(lnzi)]=0,(1- _i^-aζ T_i)^p+1[z_i^aζp_p( z_i)]=0, (66) where p_p is a polynomial of degree p. The complete weight includes the finite interval terms and the logarithms of UℓU_ . D.6 Combined integration and numerical status An explicit calculation requires the coordinates, endpoint powers and logarithms in every radiation region. The four- and five-parton weights, including all restored intervals and transfer compensations, must then be combined on common variables while preserving the local cancellation. Each term retains its map, Jacobian and Born projection. Lower-dimensional terms can be included by adding variables whose integrals equal one. If the combined density 2W_2 on the six-dimensional unit cube is integrable, it gives 2(b)=∫(0,1)6d62(,b).A_2(b)= _(0,1)^6d^6 x\,W_2( x;b). (67) For a normalized sampling density q with support on the full integration domain, one complete sample has weight w2(,b)=2(,b)q(),∼q.w_2( x;b)= W_2( x;b)q( x), x q. (68) The numerical mean and error are computed from this complete weight. The proposed local extraction then takes the form c(2)(b)=H2(b)+2(b). c_P^(2)(b)=H_2(b)+A_2(b). (69) Here H2H_2 is the EFT coefficient of δ(τ3)δ( _3) in the same Born variables, scale and normalization. Equations (60)–(68) describe one route to its radiation contribution using four-dimensional integrals over the full phase space. Possible numerical implementations include finite rescalings, explicit limiting functions where available, and integration over angles and finite radiation ratios before sampling the remaining unresolved scales. All use the complete weight with its finite compensations. Figures 10 and 11 test the radiation subtraction. Completing this three-jet contact construction and evaluating it for differential predictions is left for future work. D.7 Higher multiplicities and orders The finite-difference identity is independent of the number of Born legs. Appendix C gives radiation maps at fixed Born momenta for any fixed massless final-state jet multiplicity. Extending contact extraction requires the complete combined weight and coordinates for radiation approaching each additional Born direction, together with the corresponding EFT singular distribution. At N3LO, P2B uses an extra-jet NNLO calculation and the third-order EFT singular distribution. Repeated finite differences can remove higher logarithmic powers. Their application would require maps preserving the NNLO input’s local cancellation and a treatment of the simultaneous unresolved limits of three emissions. These extensions and their numerical study remain future work.