Local Unitarity: a representation of differential cross-sections that is locally free of infrared singularities at any order

We propose a novel representation of differential scattering cross-sections that locally realises the direct cancellation of infrared singularities exhibited by its so-called real-emission and virtual degrees of freedom. We take advantage of the Loop-Tree Duality representation of each individual forward-scattering diagram and we prove that the ensuing expression is locally free of infrared divergences, applies at any perturbative order and for any process without initial-state collinear singularities. Divergences for loop momenta with large magnitudes are regulated using local ultraviolet counterterms that reproduce the usual Lagrangian renormalisation procedure of quantum field theories. Our representation is especially suited for a numerical implementation and we demonstrate its practical potential by computing fully numerically and without any IR counterterm the next-to-leading order accurate differential cross-section for the process $e^+ e^- \rightarrow d \bar{d}$. We also show first results beyond next-to-leading order by computing interference terms part of the N4LO-accurate inclusive cross-section of a $1\rightarrow 2+X$ scalar scattering process.


Introduction
The ever-increasing need for generic and accurate Monte-Carlo simulations for collider experiments spurred the emergence of an entire subfield of the high energy physics community whose research activities are to a large extent motivated by fulfilling this demand. Over the last three decades, physicists pursued this goal by developing new computational techniques and studying the mathematical structure of perturbative computations in Quantum Field Theories (QFTs). The 1990s saw the development of the first tools 1 for the automatic evaluation of tree-level scattering amplitudes. These efforts were naturally followed up by independent groups working on the corresponding phase-space integration programs, also called event generators, for automatically computing the associated differential cross-sections 2 . The following decade then witnessed the push for the automation of the computation of Next-to-Leading-Order (NLO) corrections by many members of the aforecited groups. This resulted in the so-called "NLO revolution" that cemented the arguably artificial divide between the task of computing loop amplitudes 3 and that of regulating the infrared (IR) divergences of the phase-space integral by means of dedicated counterterms. Over the years, various strategies have been elaborated to this end, first at NLO [19][20][21] and then at Next-to-Next-to-Leading Order (NNLO) [22][23][24][25][26][27][28][29][30][31][32][33]. These efforts have been complemented by advances in the analytic computation of (multi-)loop amplitudes that mostly follow a pipeline of distinct processing steps. Amplitudes are first projected onto scalar form factors that are reduced using integrationby-parts identities  so as to be expressed in terms of a smaller set of master integrals. These irreducible loop integrals are then computed by means of differential equations [56][57][58] which under certain conditions can be solved numerically [59][60][61] or by leveraging a detailed understanding of the mathematical structure of iterated integrals [62][63][64][65][66][67][68][69][70] yielding special functions with efficient numerical representations.
The large body of work cited in this concise summary of the state-of-the-art for the computation of higher-order corrections to differential cross-sections is a testament to its many successes and importance for collider phenomenology. However, its more recent progression also signals that the traditional pipeline is arguably facing a complexity barrier that is unlikely to be overcome by incremental progress. Instead, the situation calls for a radical change in methodology for addressing the root cause of this complexity increase, that is IR divergences, more efficiently than the canonical approach. One such alternative is to part ways with the IR subtraction paradigm and aim at accommodating a more direct cancellation of real and virtual degrees of freedom. A possible avenue in that regard is that of the Reverse Unitarity [71][72][73][74][75][76] approach. This technique turns the phasespace integrals of real-emission contributions into loop integrals so as to reduce the complete set of higher order contributions down to a small set of master integrals and eventually combine their divergences at the integrated level. This method produced milestone predictions for the N 3 LO corrections to the Higgs hadro-production cross-section [77] as well as N 4 LO corrections to inclusive two-body decays [78]. However, because of its reliance on a completely analytical treatment of the loop integrals, reverse-unitarity cannot accommodate arbitrary observable functions or complicated multi-scale processes.
Another possible direction is to turn all loop integrals into phase-space integrals that can be performed numerically if the phase-space measure of resolved and unresolved degrees of freedom can be aligned so as to guarantee a local cancellation of all infrared divergences. This is the main goal of our paper and, perhaps contrary to the expectations of many, we show that it is possible to write the differential cross-section for an arbitrary process (without initial-state singularities) and at any perturbative order as an expression that is locally free of any IR singularities. We refer to this rewriting of the differential cross-section as its Local Unitarity (LU) representation. Its Local aspect stresses the applicability to arbitrary (IR-safe) observables, whereas Unitarity highlights the direct combination of real and virtual degrees of freedom, i.e. emission and no-emission evolutions, into finite transition probabilities.
A first attempt at this programme at NLO accuracy goes back to 1999 with the avant-gardist work of D. Soper [79][80][81][82][83][84] who demonstrated its practical viability by applying it to a particular differential observable of the lepto-production of up to three partonic jets. However, the success of the competing approaches discussed earlier in this introduction diverted attention away from this limited form of Local Unitarity. More importantly, its generalisation to arbitrary perturbative orders and processes as well as the computational techniques and resources necessary for its practical implementation were not available at the time. Our work aims at leveraging recent progresses in the generalisation of the multi-Loop Tree Duality (LTD) relation [85][86][87][88][89][90][91][92][93], in order to provide solid theoretical foundations to Local Unitarity and demonstrate its potential for numerical applications. In particular, the alternative Manifestly Causal LTD (cLTD) representation recently presented in ref. [93] plays a key role in our proof of local IR cancellations in Local Unitarity. In what follows, we will discuss how the radically different perspective adopted by Local Unitarity allows to not solve but rather avoid altogether traditional core issues of perturbative computations, of which the regularisation of infrared singularities is a prime example.
Since R. Feynman drew his first eponymous diagram in 1948, his graphical representations have been the cornerstone of perturbative computations. Even seventy years later, the theoretical physics community still leans on Feynman diagrams to guide its intuition, although modern techniques use them mostly as a starting point for building the amplitude integrand which is then subject to heavy symbolic manipulations. Instead, our Local Unitarity representation applies individually to each forward scattering diagram, called a supergraph, which encompasses a particular subset of interferences of Feynman diagrams. The LU representation enjoys a straightforward graphical interpretation as it is constructed entirely from the identification and characterisation of the singular surfaces of each supergraph. The type of post-processing and analytic transformations of the integrand of supergraphs in LU implies that it retains a direct correspondence to the diagrammatic origin of each contribution. For example, this helps identify contributions that are non-abelian, n fdependent, planar or of electroweak origin. Moreover, LU is formulated in momentum space which facilitates the study of any particular kinematic regime. we apply LTD to analytically integrate over all loop energies of the supergraph. We show that the resulting expression for each supergraph can be made locally finite, even though it is not gauge invariant nor does it satisfy any of the usual universal factorisation properties [94][95][96][97][98][99][100][101][102][103][104][105][106][107] relating processes of different multiplicities close to infrared limits. By relying solely on energy-momentum conservation, Local Unitarity exposes that the mechanism underlying the cancellation of infrared divergences [108][109][110] in QFT is of purely kinematic origin.
The LU representation realises local IR cancellation by construction so that there is no need for explicitly listing and regularising all singular limits and their overlaps, thereby rendering it de facto valid for arbitrary perturbative orders, both conceptually and in the practical context of numerical computations. We view this characteristic as unique to Local Unitarity and it is directly responsible for LU's universal applicability as the formalism does not depend on the particular theory, scattering process or observables considered. Another defining property of Local Unitarity is that it does not require dimensional regularisation except for its unavoidable introduction when computing the integrated countertpart of the local UV counterterms and imposing that they reproduce results in the commonly used MS renormalisation scheme. Dimensional regularisation [111][112][113] is a fundamental development from the 1970s which played an essential role in providing a convenient regulator preserving most symmetries. While it allows one to consistently characterise the divergent pieces entering perturbative computations, considering infinitesimal dimensions obfuscates some core properties of amplitudes such as chiral symmetry and the definition of asymptotic states. It also significantly complicates intermediate results and it is not amenable to numerical computations.
The unorthodox approach of Local Unitarity not only provides new theoretical insights on many aspects of the perturbative expansion of scattering cross-sections but also offers practical prospects for computing deeper perturbative corrections. For each relevant section, we first provide a detailed account of the general concepts introduced for the specific illustrative example of the NLO correction to the scattering process e + e − → dd. In sect. 2, we set our notation and introduce the core concept of organising the computation in terms of supergraphs. We then present our main result in eq. (3.38) of sect. 3 together with the proof of the cancellation of IR singularities that do not involve initial-state splittings. In sect. 4 we discuss the regularisation of the remaining UV singularities in the Local Unitarity representation. Sect. 5 provides details about the future generalisation and automated implementation of our method. In sect. 6, we give quantitative results about the performance of our numerical implementation of LU for the computation of a next-toleading accurate differential cross-section. We also explicitly verify local IR cancellations beyond next-to-leading order by computing several scalar supergraphs featuring up to five loops. Finally, we give our conclusion in sect. 7.

Foundations of the Local Unitarity representation
The fundamental object underlying any QFT observable calculated perturbatively is the supergraph [79], an entity that generalizes interference diagrams used in scattering amplitudes. More specifically, the supergraph can be seen as the representative of a class of interference diagrams which we will prove to be IR-finite when summed together.
In the following section we will show • how to rewrite any differential cross-section as a sum over cuts of supergraphs, and give a local representation of each individual supergraph that provides some first insight on how an equivalent local and IR-finite object can be defined. In particular, we will show how to accommodate a global routing common to all interference diagrams referring to the same supergraph class.
• that there is a direct correspondence between connected subgraphs of a supergraph and its threshold (pinched or non-pinched) singularities. This implies that the singularities of an amplitude can be fully characterised by couplets of connected subgraphs of a supergraph (loosely speaking, one subgraph identifies the Cutkosky cut and the other one the actual singular surface).
• that it is possible to fully characterise the notion of pinching in the LTD formalism by identifying singular surfaces for which no valid causal deformation of the loop kinematics is possible because it would either violate the continuity constraints or the singular surface has no welldefined oriented normal.
• that is it possible to accurately study the pattern of IR cancellations, described in more detail in sect. 3, by considering a simple functional analogue of interference diagrams. We will show that the sum of this analogue of each interference diagram corresponding to the same supergraph is locally finite.
In each section of our paper, including this one, we choose to first address its content by applying it to the explicit example of the NLO computation of the differential cross-section of the scattering process e + e − → dd. We find that this approach facilitates the introduction of our notation and the more abstract concepts presented in later parts of each section. This also provides the reader with some first intuition about the inner workings of the method and of its graphical interpretation.

Illustrative example: NLO correction to e + e − → dd
Our example process only involves two supergraphs from which arise all interferences of Feynman diagrams contributing to the amplitudes. We then analyse the threshold structure of these two supergraphs in their LTD representation, and write a formula for cross-sections using the LTD formalism. A locally regularised expression for the NLO correction to the e + e − → dd cross-section has been studied before in ref. [79] and in refs. [114][115][116].

The double-triangle and self-energy supergraphs
Our first task is to enumerate all supergraphs contributing to our process of interest. Each supergraph encompasses a number of interference diagrams, and the collection of all interference diagrams stemming from all supergraphs reproduces the complete set of interference diagrams whose sum -4 -  Figure 1. The two supergraphs Self-Energy (SE, 1a) and Double-Triangle (DT, 1b) contributing to the NLO QCD correction to the cross-section of the process e + e − → γ → dd.
yields the scattering probability. In the specific case of the NLO correction to the e + e − → γ → dd cross-section, we identify only two distinct supergraphs, which we refer to as the nested Self-Energy (SE) supergraph and the Double-Triangle (DT) supergraph , shown in fig. 1. We note that the selfenergy supergraph has two isomorphic occurences which we combine into a single representative weighted by a symmetry factor of two (see details in sect. 5.3).
The two supergraphs Γ se and Γ dt can be described formally by a couplet of a graph and a set of incoming edges. More precisely, Γ ∈ {Γ se , Γ dt } can be rewritten as Γ = (G Γ , a Γ ), where G Γ = (v Γ , e Γ ) is a graph and a Γ is a set of initial states. We encode the graph as a couplet of a set of vertices and oriented edges connecting them. We then write e Γ = e Γ int ∪ e Γ ext so as to distinguish between exterior edges that are connected to a degree-1 exterior vertex and the other interior edges. Since degree-1 vertices are in one-to-one correspondence with external edges, we will exclude them from v Γ . In summary, both self-energy and double-triangle supergraphs can be characterised as follows:  (7,8), (8,7), (8,9), (9,6), (9,10)} if Γ = Γ se {(5, 6), (6,7), (8,7), (8,6), (7,9), (9,8), (9,10) Whereas edges in e can in principle be identified from the two vertex labels they connect, we instead choose the more concise single label specified in fig. 1. Each oriented edge e i is assigned a four-vector q ei specifying the momentum it carries. We also define the characteristic vector of a set of edges x as follows: as well as the following sets of edges for a given set s of vertices: The set δ ± (s) consists of all edges of the supergraph with only one of the two vertices it connects being part of the set s. δ ± (s) can loosely be defined as the list of edges "entering" (resp. "exiting") the set of edges s. As mentioned earlier, each supergraph is effectively the representative of a class of interference -6 -diagrams, each associated to a particular Cutkosky cut of their reference supergraph. Each Cutkosky cut of the supergraph Γ is then defined using a subset of vertices s ⊂ v Γ that identifies two connected subgraphs δ • (s) and δ • (v Γ \ s), with the extra constraint that the initial states are contained in s, that is a Γ ∩ δ • (s) = a Γ . Perhaps more intuitively, the Cutkosky cut can equivalently be identified using the set of internal "cut" edges c Γ s = δ(s) \ a Γ whose removal divides the supergraph into two connected amplitude graphs.
Let E Γ s-ch be the set of all possible Cutkosky cuts of a given supergraph Γ. This set reads as follows for the two distinct supergraphs of our example process: Observe that we have intentionally left out the Cutkosky cuts {5} and v \ {10} as these two contributions are vanishing because, on top of having no phase-space support, they can be thought of as being subject to an observable function that is in this case identically zero on cuts containing γ. Using the definition of c Γ s , graphically represented as a line crossing all the edges contained in it, we list in fig. 2 the three (resp. four) Cutkosky cuts of the self-energy (resp. double-triangle) supergraph.
In most sections featuring this illustrative example process, we focus on the double-triangle supergraph only for simplicity in which case we suppress the upper index in Γ dt when not ambiguous.

LTD representation and thresholds of the double-triangle supergraph
Our goal is to demonstrate the local cancellation of the IR soft and collinear divergences of the double-triangle supergraph. To this end, we first identify its IR limits by re-expressing the doubletriangle supergraph integral using its Loop-Tree Duality (LTD) representation [86,93], where the energy components of the loop momenta are integrated out analytically using residue theorem: where the set B enumerates all possible momentum bases (or equivalently spanning trees) of the double-triangle supergraph: and σ bi are the cut-structure signs (see ref. [88]) assigned to each of the edges in b i . We show in fig. 3 all eight momenta bases b ∈ B of the double-triangle LTD representation together with their corresponding cut structure σ b4 . Each element in the basis corresponds to a particle becoming on-shell, i.e. it is cut. The cut structure σ b e sign that appears as superscript in the Dirac delta determines on which sheet the on-shell particle resides (the positive or negative energy solution), and is represented in fig. 3 as a plus or minus sign associated to each cut. We note that we opted here to use the original LTD representation of ref. [88], and not the Manifestly Causal (cLTD) variant of ref. [93]. While we will see that the latter plays an important Figure 3. The eight loop momentum bases (the inverse of spanning trees) resulting from applying the LTD expression of ref. [88] to the double-triangle topology. The cut structure σ b i is reported both in the drawings as signs placed next to the LTD cuts as well in the captions as a superscript of the edge labels.
role both in the proof of local IR cancellations and for the numerical stability of our implementation of LU, the former is better suited to highlight the connection of LTD with Cutkosky cuts. The LTD representation consists of a sum of tree-like graphs that is only singular on bounded, convex (ellipsoid-like) surfaces called E-surfaces in ref. [88]. Even though each individual summand (referred to as dual integrand ) building this representation is also singular on (hyperboloid-like) H-surfaces, their sum is regular on these surfaces in virtue of a mechanism known as dual cancellations [85,93,117].
We can now relate the E-surfaces of the LTD representation of the double-triangle supergraph with its Cutkosky cuts. To this end, we first recall what the elements of E Γ DT s-ch , i.e. the set of all Cutkosky cuts of the double-triangle topology (represented in figs. 2d-2g) are: We also introduce the following notation for the on-shell energy of edge e: where for this particular double-triangle supergraph the propagator masses m e are 0. An element s ∈ E s-ch can be associated to a threshold which, in Minkowski space, corresponds to a singularity of the integrand of the supergraph when energies of the corresponding cut edges take the following on-shell values: whereas in the LTD representation, the same singularity is characterised by the following E-surface: We can thus list the LTD representations of all thresholds of the double-triangle supergraph: (2.13) -8 -As we shall see, the particular signs selected for the on-shell energies of the edges being cut stems from the fact that these correspond to the only singular surfaces of the multi-loop LTD representation of the double-triangle supergraph (other than soft configurations). The LTD expression of the double-triangle also involves two additional E-surfaces that we refer to as internal as they do not involve Q 0 . Their implicit equation is: which corresponds to the two (non-Cutkosky) cuts identified from the subgraph with the set of vertices {7} and {8} respectively. The defining equations (2.14) can only be satisfied at soft points.
We define In order to show more precisely how this correspondence between Cutkosky cuts and thresholds naturally arises within the LTD formalism, we now consider the eight terms of fig. 3 whose sum M corresponds to the LTD expression of the double-triangle given in eq. (2.8).
We remind the reader that in virtue of dual-cancellations, only the E-surfaces of the form given in eq. (2.13) are singularities of M . Graphically, we identify in fig. 4 which LTD summands involve each of these four thresholds by separately highlighting the on-shell cuts due to the LTD treatment and the cut indicating a vanishing propagator. This clearly illustrates the correspondence between the Cutkosky cuts and the terms of the LTD representation of the double-triangle supergraph. We observe that the LTD terms containing the thresholds singularities corresponding to the Cutkosky cuts c s v 1 and c s v 2 already express the remaining triangle loop as a one-loop LTD expression, also including the correct energy sign for the cut e + 8,7 in the first term of fig. 4c. We stress that the key fact underlying the supergraph expression of eq. (2.7) is that the parametric equations of the four Cutkosky cuts, η s = 0 with s ∈ E s-ch , given in eq. (2.13) are the only (non-spurious) singular threshold surfaces of the 2-loop LTD representation of the double-triangle, due to the aformentioned dual cancellations. The only other divergences arise from internal propagators having a zero on-shell energy, typically leading to an integrable singularity.

Construction of the cross-section
In this section, we will define the cross section of the process in relation to the singular structure of the two supergraphs. Starting from the golden rule for cross-sections, we collect all possible supergraphs {Γ se , Γ dt } and their s-channel thresholds E Γ se s-ch , E Γ dt s-ch . We recall that s-channel thresholds of the supergraphs are in one-to-one correspondence with Cutkosky cuts and thus with interference diagrams contributing to the NLO correction to the cross-section of the process.
For every interference diagram, we express each of the two amplitudes (to the left and right of the Cutkosky cut) in their LTD representation. Furthermore, we choose a consistent loop momentum routing for all the interference diagrams corresponding to the same supergraph. As a consequence, all interference diagrams can be expressed using a common integration measure. This procedure is carefully carried out in sect. 2.2. For the moment, we state that after having performed these operations, the NLO correction to the inclusive cross-section of the process e + e − → dd can then be written as σ with where B s is the collection of all the loop momentum basis of the subgraph (s, δ • (s)) and N Γ a is the appropriately defined numerator, which may include non-trivial symmetry factors. The vector σ sbb fixes the energy flow of the on-shell particles consistently with fig. 4, by identifying δ(s)∪b∪b with the respective set of edges which are crossed by a cut or a circle in one of the diagrams of fig. 4, and the components σ sbb i ∈ {±1} with the signs associated to those circles or cuts. The identification between the terms summed in σ Γ a ,s and the cut diagrams of fig. 4 is now clear. More specifically, σ Γ dt ,s r 1 corresponds to the sum of diagrams in fig. 4a, which involves one term only, since in this case the Cutkosky cut does not leave loops on either side. Thus B s = B v\s = ∅. σ Γ dt ,s v 1 is the sum of three terms, corresponding to the three loop momentum basis of the remaining triangle loop, as depicted in fig. 4d. Note that in eq. (2.17), the propagator denominators in δ • (v \ s) \ b take the causal prescription −i because of the complex conjugation applied to the amplitude to the right of the Cutkosky cut.
We now generalise the concepts of supergraph and E-surface identification beyond this example process.

Supergraphs
is a directed and connected graph with a set e ext as external legs, such that |a| = |e ext |/2 and δ(v int ) = e ext . Roughly speaking, the edges in a correspond to the incoming particles of the process, and e ext contains two copies of the incoming particles of the process at hand.
It is convenient to describe features of the supergraph in terms of cuts, that is subsets of the set of all vertices of the graph, their boundary, that is the collections of edges connecting vertices in the cuts with vertices not in the cut and their interior, that is the subgraphs identified by the edges whose vertices are contained in the cut. Given any subset of vertices s ⊆ v, we define the following operators: As for most of the notation introduced in this section, δ always implicitly carries a dependency on the super graph Γ, which our notation will often omit for brevity. Furthermore, each subset of the edges can be fully characterised by a binary vector whose entries are 1 if the corresponding edges are in the subset and 0 otherwise. Given any subset of edges e ⊆ e, the characteristic vector of e , is defined as χ e ∈ {0, 1} |e| with Characteristic vectors allow to compactly write momentum-conservation conditions, which are interpreted as a conserving network flow. Each edge of the supergraph can be assigned with a weight q e ∈ R 4 , e ∈ e that corresponds to the momentum carried by the edge, and momentum conservation constraints can then be written as follows: Most of the constraints of eq. (2.20) are linearly dependent. In order to eliminate redundancies and obtain a minimal set of momentum-conservation constraints, it is possible to reduce the system and obtain the set of minimal constraints that is in one-to-one correspondence with the edges of a spanning tree. This fact alone is sufficient to show that the kinematic space of the virtual momenta (the linear space where momentum conservation constraints hold) is spanned by the momentum weights associated to the edges not contained in a given spanning tree. The basis corresponding to a given spanning tree can be mapped via a totally unimodular matrix to a basis corresponding to a different spanning tree.
In this framework, a Cutkosky cut admits an especially simple representation as a connected subset of the vertices, which thus allows to divide the supergraph in two connected subgraphs of it (i.e. two interfering amplitudes). The energy flow across the Cutkosky cut is enforced to be such that every edge in the boundary of the Cutkosky cut and not in e ext has an energy flow that is opposite to that of the edges in e ext themselves. More precisely, a Cutkosky cut on the FSR supergraph is a subset s ⊆ v with the following properties: -11 -• the graphs G = (s, e s ) and G = (v \ s, e v\s ), with e s = v∈s δ(v) are connected, The Cutkosky cut can be equivalently identified with the subset of edges c s = δ(s) \ δ(v int ), and graphically represented as a line crossing all edges in c s . Removal of the edges in c s from the graph G yields two connected components δ • (s) ∪ a and δ • (s c ) ∪ (e ext \ a) with s c = v \ s containing the external edges a and e ext \ a respectively. As already mentioned, these correspond to the two interfering amplitudes that form the supergraph when stitched back together. We observe that there is an apparent two-fold degeneracy since c s and c v\s identify the same Cutkosky cut. Let E s-ch be the set of all Cutkosky cuts modulo this two-fold symmetry (the name of this set will become clear later when we relate it to a subset of the threshold singularities of the supergraph in its LTD representation).
The couplet formed by an FSR supergraph and one of its Cutkosky cut c s is an interference diagram. In the perturbative formulation of relativistic quantum mechanics, the all-order crosssection is obtained by summing all supergraphs for which one sums over all possible interference diagrams arising from its Cutkosky cuts, each weighted by a density of states (i.e. observable): where O s : R 4|cs| → R is the observable function, and The formula for σ O matches the usual formula for semi-differential cross-sections, rewritten in terms of supergraphs. The phase-space integration over the final state particles is constrained by the onshell conditions associated to them. The amplitudes, however, in the Minkowski representation, feature unconstrained four-dimensional integrations for each of the loops of the graph. This apparent asymmetry in the treatment of virtual and real particles is partially lifted within the LTD representation, in which the amplitude is re-expressed as follows: where σ b ∈ {−1, 1} |b| is the cut structure for the basis b ∈ B (B is the set of all possible loop momenta basis of the subgraph δ • (s)) with respect to the reference basis b ∅ , and N µ s is a tensor polynomial in the loop variables whose (colour, spinor and Lorentz) indices are collected in the symbol µ.
When rewriting each amplitude using their LTD representations and factoring out the common loop integration measures, we find that the virtual and real degrees of freedom now appear undifferentiated. In particular, the procedure defines an extended Lorentz-invariant measure that encompasses both loop and phase-space integration, the only difference residing in the presence of an extra energy conservation condition that the external momenta of each Cutkosky cut must -12 -satisfy: This shows that quantum corrections are equivalent to performing phase-space integrals of tree processes in which virtual particles are substituted by on-shell particles whose momentum is allowed to range over all possible kinematic values (contrary to external particles, whose momenta are naturally constrained by the collision energy). Energy-momentum conservation conditions can be solved jointly for both graphs on the left and right of the Cutkosky cut by directly considering one loop momentum basis b ∅ of the complete supergraph. This allows one to write each of the extended phase-space integration measures, for varying s, b, b , as arising naturally from the loop integration measure of the supergraph, plus an extra energy conservation delta on the momenta crossing the Cutkosky cuts. It is then possible to solve all but one of the Dirac deltas and adopt a common basis for all spatial degrees of freedom: and given a reference momentum basis b Γ ∅ of the supergraph Γ. We have aligned the integration measures in b and b , which is a first important step towards proving local cancellation of IR singularities. We observe that, if a reference basis of the supergraph is fixed, the only dependence on s left in the measure element is due to the Dirac delta enforcing the conservation of on-shell energies across each Cutkosky cut. Aligning this last element of the measures requires a more advanced mathematical construction, presented in detail in sect. 3.2.4. In the rest of this section, we first discuss the singular structure of amplitudes and derive a heuristic cancellation argument, as both of these aspects do not strictly rely on the complete alignment of the measures [dΠ s ].

Identification of E-surfaces with cuts
The singular surfaces of interference diagrams are also singular surfaces of the supergraph. More specifically, E-surfaces of amplitudes correspond to intersections of E-surfaces of the supergraph, as one can think of E-surfaces of the amplitudes as E-surfaces of the supergraph evaluated on the E-surface corresponding to the energy conservation delta. Indeed, we observe that after solving the energy conservation Dirac delta in eq. (2.29), the integrand is evaluated on the zeros of the following E-surface: which is itself an E-surface of the supergraph in its LTD representation. Eq. (2.30) also suggests a useful identification of E-surfaces with connected subsets of the vertices. This graphical representation of thresholds, and the notion of connectivity, encodes the causal ordering of the scattering events (the vertices). The particles connecting two vertices on different sides of the cut generate a singularity by simultaneously lying on their respective mass shell. Cutkosky cuts, in particular, correspond to thresholds in which the incoming momenta enter the boundary of a subset of vertices and produce |c s | outgoing (that is, with opposite energy flow) on-shell particles. More specifically, let (2.31) Then an E-surface of the supergraph can be represented as a tuple (τ , α), with τ ∈ E and α ∈ {±1}. The parametric equation is which is a general equation establishing the generic form of a threshold singularity on the supergraph, for an on-shell energy flow assigned to the external edges e ext according to its bipartition in a and e ext \ a. The set E can be mapped surjectively to the set of singularities of the supergraph. This is one of the results of the Manifestly Causal LTD representation [91][92][93]118]. According to our earlier definition of Cutkosky cuts, the E-surfaces of the supergraph which correspond to interference diagrams are only those corresponding to cuts whose boundary contains all incoming edges. It is useful to decompose the set E as the disjoint union of the following four sets: • E s-ch = {τ ∈ E | δ(τ ) ∩ e ext = a or δ(τ ) ∩ e ext = e ext \ a}; an element τ ∈ E s-ch , is said to be s-channel-like.
is said to be t-channel-like.
; an element τ ∈ E isr , is said to be Initial-State-Radiation(ISR)-like.
so that E = E s-ch ∪ E t-ch ∪ E isr ∪ E int . We will consider all these sets to be modulo conjugation, that is a cut s and a cut v \ s are equivalent. The E-surfaces in E s-ch are the only ones that divide the supergraph into two graphs that have all incoming particles or none of them. The E-surfaces in -14 -E isr are typically not selected by common hadron collider observables and including them leads to interesting theoretical observations that we discuss further in sect. 5.1. The singularities of the interference diagrams with Cutkosky cut c τ , can be described as the intersection of two E-surfaces, one describing the on-shell energy conservation condition associated with the Cutkosky cut, and the other being a singularity of an amplitude. More specifically, the location of the singularities of the interference diagram is characterised by the set of points which, when embedded in R 3L , satisfies the following equations: This follows from a repeated application of the principle identifying connected cuts and E-surfaces to the amplitudes participating in the interference diagram. One interesting consequence of this claim is that no E-surface (s, α) can be a singularity of the interference diagram. This fact plays a crucial role in the determination of the cancellation pattern of IR singularities and precisely relates the singular structure of interference diagrams to those of the corresponding supergraph. For two s-channel E-surfaces corresponding to the cuts τ and s, one can also define which is an alternative representation for the E-surface η τ when evaluated at the zeros of η s . In the following, if the index α is suppressed, it is assumed to take the value α = 1. As we shall see in sect. 2.5, the property η (τ ,1) | s = −η (s,1) | τ plays a crucial role in the local cancellation of pinched singularities in LU since both c s and c τ are valid Cutkosky cuts of the supergraph.

Pinched E-surfaces and their properties
The notion of pinching has a direct physical interpretation as it is associated with the degeneracy of massless particles in collinear or soft configurations. A general definition of pinched E-surfaces characterises them as poles that cannot be deformed around via a complex contour. In order to study the conditions for which an E-surface becomes pinched, we provide each integral σ O s with its own deformation, constructed by analytically continuing the spatial degrees of freedom of the loop variables forming the momentum basis of the left and right subgraphs δ • (s) and δ • (v \ s). Given a basis b of the supergraph, we consider the contour where L j=1 α e j κ j = 0, ∀e ∈ δ(s), (2.36) which establishes that the momenta of the particles in the Cutkosky cut is kept real as they enter the observable function. It follows that understanding the pinching condition for interference diagrams is equivalent to understanding it in the context of amplitudes. More explicitly, we can analyze pinching within the object A µ s ( k − i κ). Since the deformation only affects the amplitudes, it must satisfy the four constraints laid out in ref. [89], where a general deformation for amplitudes in the LTD representation is constructed: the continuity constraint, the causal constraint, the expansion validity constraint, and the complex pole constraint.
The continuity constraint imposes that any valid deformation must go to zero faster than E i on the zeros of E i . We thus conclude that any surface in is a pinched surface due to a soft configuration (when m 2 i = 0). We now turn to the causal constraint. An E-surface is pinched if it is impossible to satisfy the causal constraint for points on it, that is with the understanding that there exists no contour deformation winding around the thresholds and that is compatible with the causal constraints. An intuitive understanding of this condition (e.g. see fig. 6 of ref. [89]) is obtained by analysing the complex zeros of the E-surface. The real space is entirely sandwiched between two complex surfaces denoting the zeros of η. The two surfaces become purely real simultaneously at the location of the pinches, as established by the complex pole constraint [119][120][121]. Any valid contour deformation is thus constrained to be identically zero at the location of the pinches.
We now turn to the identification of the E-surfaces within the (loop) amplitude A µ s on the left of a particular Cutkosky cut c s . In order to be able to characterise all possible E-surfaces of such amplitudes, we identify them with η cs (τ ,α) , and τ ⊂ s , which is part of an interference diagram corresponding to the supergraph (G, a), that is with α = ±1. We choose, without loss of generality, that all the edges in δ(τ ) have the same orientation. The imaginary part of such E-surface can be written conveniently by expanding to first order in the Taylor expansion. For any j ∈ δ(τ ) \ δ(s) identifying the dependent edge (for example, for the amplitude E-surface give in fig. 5, one can choose j = e 2 and write q 2 = q 1 − q 3 − q 4 ), we can define b τ = δ(τ ) \ δ(s) \ {j} and: The condition of vanishing imaginary part of eq. (2.38) is trivially satisfied for E-surfaces arising from real emissions only (e.g. η s | τ of fig. 5), as there are no loop variables that can be contour deformed in this case. Applying the pinched condition in eq. (2.38) to eq. (2.40), and using that the Q i 's are linearly independent by construction, we see that the requirement that the imaginary part must be identically zero for every value of the deformation field, implies that each of the vector that Q i multiplies must themselves be identically zero. That is: which can only be satisfied when q i and q j are anticollinear. Observe that for massless on-shell momenta, the two vectors are normalized to unity. On the other hand, in the massive case the two -16 -vectors have varying norms that range between zero and one. This prohibits pinching of E-surfaces with non-zero masses. We recall that soft emission of massless particles from massive ones still lead to pinching in virtue of eq. (2.37). Using eq. (2.42) and eq. (2.39) we provide a necessary and sufficient condition for η cs (τ ,α) to be pinched: • m j = 0, ∀j ∈ δ(τ ), that is, all particles participating in the threshold are massless, with all the x i ≥ 0, ∀i ∈ δ(τ ), that is, the external particles participating in the threshold are all simultaneously collinear to a given direction and all the virtual particles participating in the threshold are all simultaneously anti-collinear to that direction (when assuming that all edges in δ(τ ) are either flowing inward or outward of the set τ ).
These conditions allow to completely characterise pinched surfaces of amplitudes in their LTD representation, on top of providing the precise locations of the singularities. We recall that thresholds of interference diagrams can be seen as intersections of E-surfaces of the supergraph. Analogously, pinched points of interference diagrams, are contained in (but do not coincide with) such intersections. We show in fig. 5 and explicit example of the correspondence between the E-surface η τ | s of a supergraph and its counterpart η cs (τ ,1) in the amplitude obtained when imposing the Cutkosky cut c s . We stress that the set of vertices τ defining the amplitude E-surface η cs (τ ,1) relates to the two sets τ and s identifying its corresponding supergraph E-surface η τ | s with τ = s \ τ (when τ ⊂ s). , with τ = s \ τ = {2, 3}, which is pinched for − qi = xi p τ , qm = −xm p τ , q l = −x l p τ , ∀i ∈ {e1}, ∀m ∈ ∅, ∀l ∈ {e2, e3, e4}, with p τ = −q1. More concisely, for xi > 0, the pinch happens for q1// q2// q3// q4. The complexconjugate of the tree-level amplitude A µ v\τ on the right of the Cutkosky cut cτ has no amplitude E-surfaces, but instead features the familiar real-emission triple-collinear phase-space singularity which, within the LU framework, cancels A µ s on the pinch of its amplitude E-surface.
The t-channel type of amplitude E-surfaces never induce divergences in LU as their pinched configuration is either excluded by the observable [122] and/or regulated by the propagator width assigned to unstable particles. We are therefore interested in s-channel type of supergraph Esurfaces and assume each particle to be massless, then the set of points at which any E-surface of any interference diagram corresponding to a fixed supergraph Γ is pinched can be written as Adding masses to particles of the interference diagrams decreases the size of H and, consequently, that of P. This concludes the study of pinched E-surfaces in the LTD framework. It is worth mentioning that a proper classification of the pinched E-surfaces requires studying their intersections. However, this study is unnecessary to prove FSR cancellations and is outside the scope of this work.

Local cancellations for final-state radiation within a toy model
In this section, we render the cancellation shown fig. 5 more systematic by investigating it within a toy model. Such local cancellation requires a specific alignment of the integration measures in order to solve the Dirac deltas enforcing on-shell energy conservation. This treatment is carried out in detail in sect. 3.2.1 and we first discuss here the cancellation mechanism for a simplified model. In order to carry out the argument, we consider an analogue for the integrand of interference diagrams which is constructed in the following way: each interference diagram, identified by s ∈ E s-ch , shall be associated to a function ω s : R 3L → R which is the product of the LTD representations of the graphs δ • (s) and δ • (v \ s) times the product of the inverse energies of all the particles in the Cutkosky cut c s , where we take the the amplitude to have a scalar numerator, that is A µ s = A s , for simplicity. The singularities of ω s are the same as those of the two amplitudes integrands A s , A v\s , plus the integrable singularities due to the inverse energies of the particles in the Cutkosky cut. More specifically, the E-surfaces of ω s correspond to the zeros of E-surfaces satisfying eq. (2.33); thus ω s is a valid analogue of the integrand of an interference diagram. Furthermore, ω s exhibits an interesting local factorisation property in the neighbourhood of such singularities: following the convention set in figure 6, in the limit η τ | s → 0, each amplitude integrand can be factorised as a product of inverse energies for each element of c τ multiplied by the LTD representation of the two subgraphs δ • (s) and δ • (τ \ s). We observe that such a local factorisation property relates to analogous ones holding at the integrated level and also playing an important role in traditional computational methods.
Let A s be divergent at a location identified by the E-surface η τ . Then this local factorisation property simply reads: Figure 6. Two different interference diagrams, corresponding to Cutkosky cuts cs and cτ . The two Cutkosky cuts identify a singularity of the amplitudes, ητ |s = −ηs|τ . The direction of the edges identify the on-shell energy flow necessary for the particles in δ(τ \ s) to participate in a physical threshold.
We can now show the explicit cancellation pattern. Let us consider the sum of the functions ω s for all s ∈ E s-ch , which is the analogue of the sum of all interference diagrams of a single supergraph. Such sum, if τ ∈ E s-ch , also contains the term ω τ . We observe that Thus, the singularity at η τ | s of an interference diagram whose Cutkosky is identified by s cancels pairwise with the singularity at η s | τ of an interference diagram whose Cutkosky is identified by τ . This heuristic argument can easily be generalised to an arbitrary number of E-surfaces which vanish simultaneously, by iterating the factorisation argument and using the fact that each s-channel E-surface η s | τ can be written as the difference of two variables, each being the sum of energies in the Cutkosky cut c τ or c s . This mechanism will be studied in more detail in sect. 3.2.4. In the next section, we will show how the cancellations unfold when including observables and construct a proof. As already mentioned, this will require re-expressing the integrand of interference diagrams in a different way by solving the energy conservation delta explicitly (or equivalently, by considering a contour integration of the LTD representation of the supergraph). After the integrands are rewritten in this fashion, the cancellations will be realised algebraically, similarly as for dual cancellations.

Local cancellations of threshold (IR) singularities
In this section we present the major steps in defining a local representation of differential crosssections that is manifestly free of IR singularities. As the paper is focused on final-state radiation, we will show that such a representation is integrable on the whole R 3L excluding the initial state radiation (ISR) contributions (see sect. 5.1 for treatment of ISR).
The conceptual unfolding of the proof is summarised by the following prescriptions • Given a supergraph, construct a flow φ satisfying an ODE involving a reference vector field κ satisfying the causal constraints laid out in ref. [89]. κ is then called the causal vector field and φ is called the causal flow. Then change variables so as to make it possible to perform a contour deformation in the flow parameter t, that is along the flow, thus allowing the use of the ordinary one-dimensional residue theorem.
• Construct a local representation of differential cross-sections, that is a function σ d : Such a function is locally equivalent to summing over the discontinuities of s-channel Esurfaces of the supergraph along a flow line.
• Show that σ d allows for a cancellation mechanism analogous to the type described in sect. 2.5, and that is mathematically summarised by the partial fractioning identity • Use this cancellation pattern to derive analytic constraints on observables by requiring that σ d is finite on R 3L excluding all the regions at which an initial-state E-surface vanishes. Next, we show that they are satisfied by observables that cluster particles with energy or relative direction under a mathematically well-defined scale δ. In other words, these constraints match the usual requirement of IR-safety for collider observables.
• Derive the scaling of σ d near soft points, and show that it depends on the scaling of the deformation field around soft points, thus relating the request of integrability of σ d with a constraint on the scaling of the causal flow.
• Perform power-counting and argue that soft points are always integrable in physical theories.
Before detailing the proof in sect. 3.2, we construct the LU representation of the e + e − → dd example, which we use to unfold explicitly the steps presented above.

Illustrative example: NLO correction to e + e − → dd
We start by recalling the LTD representation of the double-triangle supergraph, and rewrite eq. (2.16) and eq. (2.17) directly in terms of the LTD representation of that supergraph.
Observe that in the rest frame of γ , the two threshold E-surfaces η s r 1 and η s r 1 have the same exact functional dependence in the loop variables. Accounting for this accidental degeneracy, we can write eq. (2.7) as .

(3.3)
where the Cutkosky cuts {s v 1 , s v 2 , s r 1 } of the double-triangle are defined in eq. (2.9) and the addi- We now apply Cutkosky cuts to the LTD representation of the double-triangle by substituting η −1 s with δ(η s )O s for each Cutkosky cut s ∈ E s-ch . This yields a representation that is trivially equivalent to that of eq. (2.16) and eq. (2.17), after all the Dirac deltas except for the one imposing the conservation of on-shell energy flowing across the Cutkosky cut s have been solved. We have where the sum runs over all possible Cutkosky cuts, and O s is, for now, a non-specified function whose functional form depends on the cut s. It is clear that eq. (3.4) can be obtained from the usual form of dσ dO by applying LTD to the energy integrals left after applying the phase-space cuts. At the same time, eq. (3.4) also expresses the well-known fact that Cutkosky cuts can be seen as the residues of the supergraph acquired by contour-deforming around its thresholds, which is also the core of the original derivation by Cutkosky [123].
We would like to solve the Dirac deltas on the right-hand side of eq. (3.4) simultaneously for all Cutkosky cuts of the double-triangle topology, in a way that allows to write where σ d now contains no Dirac delta. In particular, the contribution from all interference diagrams stemming from the double-triangle topology should be written as a single integral over R 3 × R 3 of a particular integrand. In following three sections we will discuss a general method to solve the remaining delta function encoding energy conservation. We denote with dσ Γ /dO the sum of all the interference diagrams arising from Cutkosky cuts of the supergraph Γ and we suppress the Γ superscript henceforth.

Soper's rescaling for solving conservation of on-shell energies
For 1 → N processes, such as (effectively) e + e − → dd, final-state singularities can be aligned at any perturbative order in QCD using Soper's δ solving strategy [79] which offers an easy way to rewrite the phase-space measure in a form where there is no Dirac delta anymore. It was presented for the first time for integrands with conformal symmetries. In the following we will generalise it to multi-scale integrands, for arbitrary masses, loop momentum routings and Lorentz frames. Consider the integral (3.4), and multiply it by the integral in t of a normalized function h(t): We can now change variables from ( k , l ) to φ(t; k , l ). We call φ the causal flow, because for any fixed k and l , φ denotes a curve which always flows outwards with respect to the Cutkosky -21 -cut E-surfaces. This new object will be described in more detail in sect. 3.2.1. In the case of the example process considered in this section, we can choose 5 φ(t; k , l ) = t( k , l ). An illustration of this causal flow is given in fig. 9. This change of variables then yields: and by solving Dirac's delta explicitly, we find where t s , is the unique value of t such that the the E-surface identified by the energy conservation delta vanishes. From that point onward, we will abbreviate our notation by defining k := ( k , l ). For a given k, t s is a function satisfying Observe that for every k, t s is the factor with which to rescale the loop momenta in order for t s k to lie on an E-surface. Alternatively, one can think of fixing a point in loop momentum space and dilate or contract the E-surface by a quantity 1/t s so that the point k lies on it. If the E-surface contains the origin, then for every k there is only one positive value of t s such that t s k lies on the E-surface. This is a consequence of the convexity of η as a function of t. Knowing that there is at most two solutions if t is allowed to take positive and negative values, and since a solution t s satisfies we conclude that the equation in t can be solved numerically by using Newton's method with seeds provided by the bounds of the inequality in (3.10). Since there are at most two solutions, Newton's method is guaranteed to converge. Thus, it follows that Soper's δ solving strategy is numerically straightforward to implement.

LU representation of double-triangle interferences
The next step is to relate the expression of the double-triangle supergraph given in eq. (2.7) to that of the traditional expression dσ dO of the NLO QCD accurate differential cross-section of the scattering process e + e − → dd. The correspondence between the contributing threshold singularities of the LTD expression of the double-triangle supergraph and its Cutkosky cuts ( fig. 4) shows that dσ dO can be obtained by computing the residues associated with each of the causal surfaces for which we must however make sure to assign the observable function with the appropriate dependence. The fact that each Cutkosky cut involves the observable function with a different functional dependence on the kinematics is the very reason why the residues from each of the singular surfaces of the supergraph must be computed separately. Indeed, when only interested in the fully inclusive cross-section of 1 → 2 processes, one can instead consider directly computing the imaginary part of the supergraphs and extract from it the inclusive cross-section via the optical theorem (see, e.g. ref. [78]). 5 In this section, contrary to sect. 3, we choose the argument of the causal flow to be log(t) and not t for simplicity .
-22 -Eq.(3.8) can then be written as where we recall that we define k := ( k , l ) and where we used the formal definition of the residue of a single pole located at t s , and The symbol E k,φ describes the ensemble of all threshold surfaces, i.e. Cutkosky cuts, which have a solution in t given the change of variables induced by the causal field φ(t; k) and for the particular sampling point k = ( k , l ) considered. Due to convexity, the number of solutions in t of the E-surface equation is limited to one or potentially zero. However, in our simple 1 → N case, this change of variable amounts to a trivial rescaling which always offers exactly one solution for each threshold, so that we can effectively consider a summation over the complete set where K 1 is the Bessel function of the second kind. This particular choice of a normalised function h(t) is motivated by the fact that it vanishes exponentially fast at zero and infinity while being maximal at t = 1. We fix the tunable parameter σ to 1. These features naturally drive the integrator to probe points in space that are close to at least one Cutkosky surface (see sect. 5.4.1 for more details regarding integration efficiency) and also avoids spurious integrable singularities at t = 0. This also implies that the norms of the loop momenta before and after the rescaling treatment are of similar order of magnitude, thereby maintaining the direct interpretation of the location of the IR and UV domain. We are now equipped to delve into the details of the cancellation of IR singularities and the numerical implementation of the double-triangle and self-energy supergraphs within the formalism of Local Unitarity.

E-surface cancellations for the double-triangle supergraph
In general, the direct evaluation of eq. (3.12) reads: where we substituted each Cutkosky cut threshold surface η τ , τ ∈ E s-ch in the denominator by its Taylor expansion around its on-shell solution t τ (so that the zeroth-order term of this expansion is necessarily zero). We use the short-hand notation ∂ t η s (t s k) := ∂ t η s (t)| t=t s k . Thanks to the simple functional form of the rescaling change of variable as well as the fact that all propagators of the double-triangle supergraph are massless and that we are considering an incoming momentum configuration at rest, i.e. Q µ = (Q 0 , 0), we can give a simple expression for the rescaling solution t s as well as the derivative function ∂ j t η τ (t k) of the Cutkosky surfaces: In our case, this expression of course reproduces the exact expression of η τ (t s k, Q 0 ) since the Taylor expansion terminates, but this is in general not the case for more complicated causal flows or when in presence of massive propagators.
Proving the local cancellation of IR singularities amounts to demonstrating that the expression of the differential cross section σ d in eq. (3.14) is integrable except for remaining UV divergences. From our earlier discussion, we can already argue that f is free of any singularity. The causal nature of the field φ inducing our change of variable also insures that |∂ t η s (t s )| = 0 (see sect. 3.2.1). In our case, we have: Finally the product of on-shell energies i∈ein E i and the product of inverse internal E-surfaces In conclusion, we can rewrite σ d by underlining that the denominator can be written as a polynomial in the energies and in the η s , s ∈ E s-ch : We can now rewrite σ d and show that the cancellations can be made explicit algebraically in the denominators (t i − t j ), thereby sidestepping their complicated kinematic dependence. One peculiar implication of this proof is then that it can be made for any parametric kinematic point k, that is without taking any limit. This is particularly convenient given that enumerating all possible IR divergent kinematic limits becomes cumbersome at high perturbative orders.
The summand of σ d corresponding to s, as given in eq. (3.21), is clearly singular at the locations η τ | s = η τ ( k, 0) − η s ( k, 0) = 0 and at soft points. That is, regions of kinematics space at which any energy vanishes. In order to make the notation less heavy, let us call x s = η s ( k, 0) and suppress the dependence of functions unless their dependence is itself dependent on s, which is an index that is summed over.
In the absence of numerators, we could immediately show that the sum in eq. (3.21) is identically zero thanks to the general partial fractioning identity given in eq. (3.1). If observables are nontrivial, instead, we have to expand the numerator in the variable t. Let us assume, in the following, for fixed > 0. This condition, which will be discussed in detail later, allows to state that in a neighbourhood of problematic points where g is a continuous function on R \ {0} in virtue of the continuity of the observable, of the numerator and of the normalising function. Its singularity at the origin is not problematic if we use that the integrand must initially be UV convergent, either on itself or with the aid of counterterms. With this in mind, we can write We will now rearrange the sum in a way that makes cancellations manifest.
Written in this form, it is clear that away from soft points, when A power-counting procedure can be set up to show that this integrand has at most integrable singularities at soft points. Such an analysis, however, requires studying the structure of g and is performed rigorously in sect. 3.2.7. The straightforward cancellation structure that manifests itself in eq. (3.24) already alludes to its generalisation. Also, except for considerations regarding the observable dependence, this cancellation pattern does not discriminate between types of singularities (pinched or non-pinched) and holds on intersection of singular surfaces.
Finally, we go back to the condition shown in eq. (3.22), which enforces the IR safety of the observable. We will assume that the observables only depend on the size of the cut and the masses of the particles belonging to the Cutkosky cut, other than their momentum, that is: (3.25) When rewriting (3.22) explicitly for a pinched singularity and the momentum routing shown in fig. 2d, with Q µ = (Q 0 , 0), we find: which is the familiar IR-safety condition that relates observable functions O n acting on kinematic configurations of different multiplicities n in the soft and/or collinear limits of its massless constituents. Thus the constraint on observables implies by the request of finiteness of σ d on pinched singularities can be satisfied in the usual way of constructing observables. We can also consider the potential singularity at x s v 1 − x s v 2 = 0, identifying the phase-space points satisfying | k | = | l |, (see eq. (2.13)). This corresponds to a configuration on the non-pinched E-surface of the one-loop triangle on the left (right) of the Cutkosky cut c s v 1 (c s v 2 ). The study of the regularity of the differential cross-section at these threshold singularities seems at first glance to follow completely analogously to that of the IR pinched singularities. Perhaps surprisingly, this implies that we also expect cancellation between the two threshold singularities of the Cutkosky virtual contributions illustrated in fig. 7. We can follow the exact same study of cancellation between the IR pinched surfaces in the sum of terms in the r.h.s of eq. (3.21), with the only qualitative difference being the resulting condition on the observable function: It is clear that this condition is of a completely different nature than the IR safety condition obtained in (3.26). Indeed, the limit lim | l | → | k | does not imply any degeneracy in the experimental signature since one can obviously resolve the directional information of the quark and anti-quark in the final state. Consequently, one should in general expect observable functions not to satisfy eq. (3.27), implying that for non-trivial observable functions, considering the contour deformation discussed in sect. 5.3.4 is necessary. A couple of key additional points are in order: • When considering fully inclusive cross-sections, observable conditions similar to eq. (3.27) are always satisfied, so that the computation can be performed without considering a deformation. Note however that omitting altogether certain Cutkosky cut contributions effectively amounts to setting the observable to exactly zero for them. This implies that even for the computation of the inclusive NLO cross-section of the scattering process e + e − → γ → Htt, a deformation would be warranted since all Cutkosky cuts involving only the two top quarks are not considered given that the observable demands a final-state Higgs.
• Because the two observable functions on each side of the condition eq. (3.27) do not share any loop momentum dependency, it is clear that in the case of the double-triangle supergraph, only a constant observable function can satisfy it (i.e. inclusive measurement).
• Observable functions often consists of only products of Heaviside functions (e.g. implementing phase-space cuts and/or binning into histograms), in which case this opens the possibility of investigating dynamically at run time and for each integration sampling point what are the Esurfaces whose pair of canceling Cutkosky cuts are not either both selected or both removed by the observable definition. Then, the contour deformation for this point only needs to consider those surviving E-surfaces.

Proof of local cancellations of threshold (IR) singularities
In the following sections we will construct the LU representation for a generic differential crosssection and require it to be free of non-integrable singularities. We show that the resulting constraint on observable functions is satisfied by IR-safe observables. This results in a systematic proof of local IR cancellations within the LU representation of differential cross-sections.

Causal flows
In section 2 we introduced cross-sections by referring to them as weighted sums of interference diagrams, each of which is associated to a well-defined Cutkosky cut, which in turn corresponds to a Dirac delta imposing the conservation of on-shell energies. For a fixed and positive center of mass energy, the equation imposing conservation of the on-shell energies is the equation of an E-surface. These deltas, in Cutkosky's original derivation [123], arise from contour deforming the energy around thresholds of the diagrams. However, such a derivation obscures the important subtlety that the energy variables of the particles in the Cutkosky cut are linearly dependent and thus cannot be used as independent integration variables. In the following, we will provide an alternative derivation of Cutkosky cuts, by reducing the integration along thresholds to a one-dimensional problem. This derivation of Cutkosky cuts will expose the local structure of the integrands and the location of their singularities. It is formulated such that the cancellations of all divergences related to E-surfaces are local, thus allowing one to construct an integrand which is locally free of divergences by aligning the integration measures supported by the different Cutkosky cuts. Furthermore, it shows how considering transition probabilities, rather than amplitudes, is required for the infrared structure of observables to be completely understood. Let be the LTD representation of a supergraph Γ, with The discontinuities of M Γ across its s-channel thresholds represent distinct summands that contribute to define the total probability of the process whose initial states are fixed to be the particles in a and whose final states are, for now, unspecified. The function f , which is the LTD representation of the supergraph multiplied by the product of all energies and all E-surfaces appearing in the representation itself, is finite for any value of the loop momenta. We introduce a dummy integration variable t by introducing unity as the integral of a normalized function (3. 30) and then introduce the change of variables k = φ(t, k ), with φ being the solution to the following first-order system of differential equations: where we introduce the vector field κ : R 3L → R 3L which we require to be Lipschitz-continuous. We will use the map φ t : k → φ(t, k) to change the parametrisation of the phase space integral. If the zeros of κ form a zero-measured subset of R 3L with respect to the Lebesgue measure in R 3L , then the change of variables is well-defined. Thus, we can exclude all the sinks, sources and ridges that the flow φ t may have from the integration. Therefore, we can write We are now interested in performing contour integration in the variable t. In order to do this, observe that for each k, the curve φ(t, k) crosses a number of E-surfaces. However this alone does not allow to determine how many times a specific E-surface is intersected by a curve in the flow, and if approaching an E-surface along a curve yields a simple pole in the integrand. If the pole is simple, then there is a well-defined principal value procedure associated with it, and the sign of the imaginary part acquired in the contour integration is fixed by the sign of the Feynman prescription. It is possible, however, to construct a flow whose properties make the answer to these two questions manifest. Specifically, consider solutions to flow ODEs where the vector field κ is chosen to be causal, that is If that is the case, then the flow will consequently have three properties: • ∀ k, there exists at most one value t s s.t. η s (φ(t s , k)) = 0, which follows from the fact that the curve φ(t, k) cannot flow outward and inward of the E-surface η s without violating the causal prescription of eq. (3.33).
• sign[∂ t η s ] = sign[∂ t η τ ], ∀φ(t, k) ∈ δη s ∩ δη τ , also trivially guaranteed by eq. (3.33), where we define δη s := k ∈ R 3L η s ( k) = 0 . With a slight abuse of notation, we write ∂ t η s (φ(t s , k)) = ∂ t η s (φ(t, k))| t=t s . The first property determines the number of intersections a curve has with a determined Esurface. The second determines that all poles appearing on the real t axis for fixed k are simple. The last one, although momentarily obscure, will be fundamental in order to realize local cancellations of pinch singularities. We stress that the first two conditions are not strictly necessary in order to build a valid contour integration in t. Indeed, regarding the first condition, there is nothing in the principal value procedure that forbids us to contour integrate around two distinct poles. Regarding the second condition, one could think of excluding from integration the regions of space at which ∂ t η(φ(t, k)) = 0, which lie on a zero-measured set, and then establish if the integration is finite.
Given the first property, for every k one can define the set which contains all the E-surfaces which are intersected once by the curve φ(t, k). Thus, for each s ∈ E k,φ , we can write the expansion of η s around its unique zero, t s : (3. 35) and the first order in the expansion is ensured to be non-zero by the second property. In conclusion, we observe that the existence of a causal vector field κ being causal on all the E-surfaces of the supergraph Γ is guaranteed by the work carried out in ref. [89]. Since the causal vector field, as constructed therein, is infinitely differentiable, Sard's theorem also ensures that its zeros lie on a zero-measured surface. For future convenience, given a set of points V ⊂ R 3L , we define a set containing all the points that can be mapped into V by the causal flow The inverse image of the causal flow is fundamental in determining the analytic properties of the supergraph and how they relate to those of residues in t of the supergraph as parametrised along the flow.

Visualisation of the causal flow
In figs. 8a, 8b and 8c the causal flow of a particular one-loop example, called Box_4E in ref. [89], is shown. Since the origin is not on the interior of all E-surfaces, the rescaling strategy defined in section 3.1.1 cannot be applied. For all 1 → N processes the simple rescaling flow is applicable, and we show the causal flow of Box_4E only for the purpose of illustrating the complications arising for a non-trivial overlap structure of E-surfaces, such as the one that can appear for challenging supergraph topologies of 2 → N processes. In that case, it is likely that the system of ODE defining the causal flow requires a numerical solution, but all properties of the LU representation discussed in this paper would still hold. The causal flow of the double-triangle supergraph is more delicate to represent than that of Box_4E, since there are two loop momenta in that case. We choose to project the six-dimensional input space ( k , l ) of the DT topology with Q µ = (2, 0, 0, 0) onto the parametric plane λ kêk + λ lêl withê k = ((0, 0, 1 2 ), (0, 0, 0)) andê l = ((0, 0, 0), (0, 1 2 , 1 2 )). This section involving a non-constant momentum component is convenient since it necessarily contains the image of a rescaling flow, thus allowing to render the flow within this same plane, like it was the case in fig. 9. An important remark is that the E-surfaces η s v 1 and η s v 2 are unbounded. This is a result of their independence of λ k (resp. λ l ) which only controls the one-loop integration volume of the triangle loop remaining on the left (resp. right) of the Cutkosky cuts η s v 1 and η s v 2 . This is what allows the LU representation of the DT supergraph to probe the UV regime with a rescaling parameter of O(1). In contrast, the volume described by the E-surface η s r 1 results from an equation involving the sum of three square roots and is therefore quite complicated but still bounded since it corresponds to a particular hyperplane of the three-body decay phase-space volume Q µ → (k µ , k µ − l µ , l µ − Q µ ), which we know to be necessarily contained within a sphere of radius Q 2 . Notice however, that even for a sampling point ( k , l ) with arbitrary large moduli | k | and | l |, the global rescaling flow will always yield a contribution for the real-emission η s r 1 and η s r 2 Cutkosky cuts. However these will be associated with very small values of the parameter t s r 1 and t s r 2 yielding a contribution exponentially suppressed for our choice of normalised function h(t) of eq. (3.6). For the particular projection chosen for fig. 9, the integrable singularities at k = 0 and l = 0 also correspond to vanishing values of the corresponding rescaling solution and will thus also be suppressed. These observations underline the appealing feature of LU that the change of variable k → φ(t, k) does not affect the interpretation of what kinematic region of the cross-section is probed since contributions from t i solutions far from one are exponentially suppressed. Moreover, in the particular case of a global rescaling causal flow, we even observe that the change of variable retains the collinear and soft properties of the -29 -  with V being the circled dots. The vector field depicts the orientation of the contour deformation vector constructed according to ref. [89] (using the four deformation sources indicated with blue dots), and its colour pertains to its relative magnitude.
momenta of the (pairs of) partons in the input configuration. This allows one to easily probe potentially singular regions, which we investigate in sect. 6.1.2 by choosing a different projection of the double-triangle sampling space.

The Local Unitarity representation of differential cross-sections
In his original 1960 paper [123], Cutkosky recognised that interference diagrams can be seen as the imaginary part acquired by contour deforming around specific thresholds of supergraphs. The equation of the threshold appears in the energy conserving delta as a result of the principal value identity that establishes the functional form of the imaginary part acquired by contour deforming around a pole on the real axis. Moreover, every supergraph Γ gives rise to a number of squared amplitudes obtained by summing over the imaginary parts acquired by contour integrating along all of its thresholds corresponding to cuts s ∈ E s-ch . Such a derivation does not detail the subtleties of setting up a contour deformation programme for the supergraph. For example, it does not specify what constraints the chosen contour needs to satisfy and what integration variables can easily be chosen to transform the real integration into a contour integration.
We will now show that is is possible to see Cutkosky cuts as the discontinuities of the LTD representation of the supergraph along the one-dimensional flow line of a causal flow. This allows us to explicitly construct a local representation of differential cross section from which the second golden rule arises naturally as a consequence of identifying picking up the residues of the LTD representation of the supergraph along the flow with explicitly solving the Dirac delta associated to Cutkosky cuts. Thus the work in this section should be considered to be equivalent to that of aligning the integration measures dΠ s and, specifically, resolving Dirac deltas expressing energy conservation across on-shell physical particles.
In the following, we will assume that we have a causal flow φ(t, k). For every fixed k, the unique curve φ(t, k) passing through it will intersect a subset of all E-surfaces. Recall that we defined E k,φ -31 -to be the set of such E-surfaces in eq. (3.34), which can be expanded as shown in eq. (3.35). For all s-channel E-surfaces not in E k,φ , that is for every τ ∈ E s-ch \ E k,φ , we have that η τ (φ(t, k)) = 0 for every t. Thus, every E-surface not contained in E k,φ never has a vanishing zeroth order in the expansion on the flow line. Thus, given a supergraph one can define a locally finite volume form as and O s defines the observable and is a function of k through the dependence on cut edge momenta, At the integrand level, σ Γ d reproduces what we obtained in sect. 2 through LTD and the golden rule for differential cross-sections in eq. (2.25). The index associated to the residue is fixed by the convention Ind s = +sign[∂ t η s (φ(t s , k))]. Fixing this choice is equivalent to choosing a Feynman prescription for the s-channel thresholds which are included in the definition of σ Γ d ; this arbitrary choice was already discussed in ref. [123]. The relation defines the Local Unitarity (LU) representation of the differential cross-section dσ/dO. In the following, we will sometimes refer to σ Γ d as σ d , without the explicit reference to the supergraph Γ. Yet another feature of eq. (3.38) is the presence of unregulated E-surfaces, that is E-surfaces for which the prescription can be set to zero. This is a consequence of σ d not being divergent at the location of these singularities, that is: This is related to the property of f that establishes a necessary condition for the multiple limits of the LTD expression when approaching the intersection of many E-surfaces to be zero. More specifically, given a set of E-surfaces E ⊆ E intersecting on a surface, one can ask how f ( k) behaves in a neighbourhood of any point on the intersection surface. We will show that approaching the intersection δη s ∩ δη τ , f does not vanish if and only if s ∩ τ ∈ {s, τ , ∅}, that is if one of the corresponding set of vertices defining the cuts is contained within the other, or their intersection is the empty set. This property was assumed in the usual construction of Cutkosky cuts, and was used extensively in sect. 2. In order to generalise and properly describe this property we define cross-free families of subsets. A cross-free family is a laminar family of subsets of the vertices such that no subset in it can be written as the union of two or more sets in the family itself. A laminar family is a family of sets such that for any two set, the intersection is either one of the two sets or the empty set. These two definitions imply that a cross-free family is a set of |v| − 1 subsets of the vertices. Let us define the set of all cross-free families of connected cuts of the graph: F = {F ⊆ E|F is a cross-free family}. (3.43) An example of a cross-free family is illustrated in fig. 10. We then write where c F is a polynomial in the energies and the spatial loop momenta. Observe that we do not claim c F to be non-zero; instead, the claim should be interpreted as excluding a number of contributions to M Γ . Determining the coefficients c F is a non-trivial task which can be carried out with the treatment of ref. [93], and is not strictly necessary for the proof. Finally, we observe that it is also arbitrary that only a subset of thresholds of the supergraph should be selected to contribute with their discontinuities to σ d . In the current way of defining differential cross-sections, all the residues associated to thresholds with δ(s) ∩ a = a are completely ignored: the sum in eq. (3.38) runs over elements which also are in E s-ch , and thus all the thresholds in E isr are left out of the definition.
The price to pay for this choice is the appearance of Initial-State Radiation (ISR) singularities in observables, which are commonly treated within the collinear mass factorisation paradigm. In sect. 5.1 we will consider how the construction of σ d changes and what could happen if observables which allow for all thresholds of the supergraph are included.
For now, as anticipated, we show that σ d is bounded on R 3L \ B isr, where An important subtlety in the definition of σ d , in eq. (3.38), is the possibility of taking the limit for (t − t s ) → 0 of a quantity that has raised propagators: this makes σ d ill-defined as the limit would yield infinity. Raised propagators and their treatment varies according to their origin, following a general principle: that the degenerate case is obtained as the limit of the non-degenerate case. More specifically: • Raised propagators that are due to external self-energy insertions can be dealt with by showing that after mass renormalisation in the on-shell scheme, the contribution of a self-energy diagram equates that of terms of order one or higher in a Taylor expansion around the on-shell condition. This realisation allows us to effectively lower the raised propagator by one power, thus solving the problem. An extensive analysis of this method is performed in sect. 4. Raised propagators that are due to internal energy insertions can be dealt with by adding a different fictitious shift in each of the raised propagators, calculating σ d with this configuration, and then taking the limit in which all the fictitious shits vanish.
• Two or more distinct E-surfaces intersecting in a non-empty region can give rise to a locally raised propagator (the first-order expansion of the E-surfaces in the flow variables t is proportional to (t − t ) for all the E-surfaces which vanish at the intersection, with t being the flow location of the intersection point). In this case we impose the behaviour of σ d to be the continuation of the non-intersecting case. More specifically, since the deformation field κ, related to the causal flow by eq. (3.31), is non-zero on any point lying on a non-pinched E-surface, we are ensured that the set of curves passing through an intersection point is zero-measured. If we define then I is zero-measured with respect to R 3L , as argued before, and thus proving that σ d is integrable on R 3L \ B isr, \ I itself proves that the integral of σ d on R 3L \ B isr, is finite.
• For a particular choice of the collision reference frame there can be E-surfaces corresponding to different Cutkosky cuts but which take the same functional form (e.g. see η s r 1 and η s r 2 in the example process discussed in sect. 3.1). These raised poles are spurious and can be eliminated by taking advantage of Lorentz invariance to change the frame.
Finally, we observe that σ d is continuous on R 3L \ B isr, only if soft configurations are also excluded, since at these locations the inverse energies (other than, possibly, elements of E int ) of massless particles diverge. We define where S x is defined as in eq. (2.37). The definition of W is manifestly a zero-measured set with respect to R 3L . We then consider σ d to be defined on R 3L \ B isr, \ (I ∪ W). We will always make a distinction between the sets I and W, as it makes the unfolding of the proof easier to understand.
We shall now show that σ d is integrable on the R 3L \ B isr, \ (I ∪ W) space, which implies that σ d has a finite integral on the whole R 3L \ B isr, space.

Cancellation of pinched surfaces
We now investigate the mechanism that leads to the cancellation of pinched E-surfaces (and Esurfaces in general, in the fully inclusive case). As anticipated previously, the mechanism relies on the generalisation of the straightforward algebraic relation in eq. (3.1): Cancellations are explicitly shown by expanding all the summands in σ d in the proximity parameters t s − t τ . This leads to large expressions, with multiple sums and apparently complicated coefficients; however, these expressions allow to factorise sums in the form of eq. (3.49), and thus to immediately prove cancellations. This proof thus has the remarkable advantage of applying to any type of singularity. Differentiating pinched and non-pinched thresholds is the only important distinction and is only instrumental for characterising what are the IR-safety constraints imposed on the observable definition. We stress that cancellation of pinched E-surfaces, by itself, does not guarantee that σ d is integrable on R 3L \ (I ∪ W) \ B isr, , and specifically at the locations of soft singularities. However, it shows that purely collinear singularities (that is, excluding the boundaries of a pinched singularity, the soft points) of the summands of σ d are subject to cancellations, and that at soft locations, the resulting power-counting is less severe. The proof is thus concluded only after a power-counting analysis of soft singularities (which is carried in subsect. 3.2.7).
Another important feature of σ d is that the E-surfaces are associated to Feynman prescriptions which could complicate the proof for complex kinematics, e.g., in the presence of a deformation. In the following we will drop the prescription and discuss the case of a non-zero deformation in sect. 5.3.4. The limits in the definition of σ d , for k ∈ R 3L \ (I ∪ W) \ B isr, , can be performed explicitly and yield where We remind the reader that the explicit unfolding of the dense eq. (3.50) for the double-triangle supergraph of our example process is given in eq. (3.14). We will now rewrite σ d in a more convenient form, in which the divided differences (3.49) are manifest: . (3.53) We observe that, under appropriate conditions on h, w(t s , k) is bounded for every k. Furthermore, observe that we have factored out the product of the absolute value of the first order derivative of the E-surfaces in E kφ , which can only be done when the following condition is satisfied: which is guaranteed by the choice of a causal flow. Specifically, the causal flow is constructed from a deformation which is explicitly constructed to satisfy the causal prescription on all E-surfaces and on their intersections. In this approach, we see how the realisation of IR cancellations is interlocked with causality, enforced through the definition of the Feynman prescription. Eq. (3.52) is exactly in the form needed to apply divided differences. We will now set the stage to derive the conditions under which the observable function O s preserves the IR cancellations established by the divided differences, and thus can be considered IR-safe. At first, let O s = 1, ∀s ∈ E s-ch , so that cancellations are trivially realised, since we can immediately apply the divided difference relation of eq. (3.49) to and argue that all the multiple singular limits other than soft limits yield a finite result. Thus, in the fully inclusive case, σ d is integrable on R 3L \ B isr, \ (I ∪ B (W)), that is on R 3L \ B isr, \ B (W).
The simple pattern of cancellations in eq. (3.49) immediately applies for an arbitrary perturbative order and in the case of fully inclusive observables. For non-trivial observables, the situation is more complicated. In particular, we observe that observables carry an explicit dependence on s, and not just indirectly through their dependence on t s : the functional form of the observables is dependent on the Cutkosky cut. Thus for arbitrary observables, the pattern of cancellations in eq. (3.49) does not necessarily hold. Instead, the request of IR-finiteness implies constraints on observables.

IR-safe observables and infrared scales
In this section we investigate the role that non-trivial observables play in the mechanism of cancelling divergences discussed in the previous section. Since observables discriminate Cutkosky cuts and are intrinsically local objects, proving that the sum is finite for non-trivial observables is not a purely algebraic matter like in the fully inclusive case, but can only be done in a neighbourhood of the problematic points. In principle, for what concerns the proof, an observable can depend on the supergraph Γ, the Cutkosky cut s it corresponds to, the whole set of loop variables of the supergraph k, and the kinematic configuration of the initial states. In practice however, we defined observables as having a dependence on the momentum, mass and quantum numbers of only the particles in the Cutkosky cut. Thus, we will first work out a general sufficient condition for the cancellation pattern to be realised, and then enquire if and how observables satisfying eq. (3.40) satisfy these general constraints. We will impose these constraints to be flow-independent, since the choice of a causal vector field should not affect the physical properties of observables.
-36 - We start by observing that σ d is bounded on This is required for the divided differences in eq. (3.49) to be applicable. This condition is manifestly dependent on the flow. We will use the following relation: if a function f satisfies a property on an open set V , then the function f • φ satisfies it on Φ −1 [V ]. Indeed the flow projects a set B (I) on a neighbourhood of an intersection δη τ 1 ∩ δη τ 2 , and we conclude that a sufficient condition for eq. (3.56) to hold is that where B (δη τ 1 ∩ δη τ 2 ) is a toroidal neighbourhood of the intersection between the two E-surfaces, and is a quantity dependent on the experimental setup. The interpretation of is that of a natural resolution below which the experimental setup is unable to distinguish degenerate kinematic configurations of final states, identified by τ 1 and τ 2 . In conclusion, the notion of degeneracy and its relation with experimental resolution is derived here as a result of the analytic properties of the local representation of differential cross-sections. We explicitly verified that constant observables of the form of eq. (3.40) can retain local cancellation of non-pinched E-surface in the entirety of the integration volume. Instead of performing the same analysis for pinched E-surfaces and showing that IR-safe observables can always be defined, we will limit ourselves to stating how the condition relates to the common definitions of d IR-safe collider observables.
It may now seem that for non-trivial observables it is not possible to make σ d integrable on the whole R 3L \ B isr, , as non-pinched threshold are not subject to cancellations. However, nonpinched E-surfaces can be integrated using an ad hoc contour deformation, which cannot be done for pinched E-surfaces (for which the cancellation pattern must still hold, because of IR-safety). The deformation of amplitudes was discussed in ref. [89] and we discuss how to extend it to cross sections in sect. 5.3.4. In the following, we will assume that such a contour deformation can be constructed and that proving σ d to be integrable on B (H) implies that σ d can be analytically continued and its non-pinched thresholds contour integrated so that it is regular on the entirety of We observe that in order to preserve the cancellations at the pinched points, it is sufficient to require that If we use the constraint in eq. (3.40), this equation takes a more illuminating form: with α i , β j , i subject to momentum conservation. Eq. (3.59) makes manifest that the ability of observables to satisfy IR-safety constraints is necessarily associated with momentum-conservation laws on the graph; it can only be satisfied if the observables themselves cluster together all the particles of their defining Cutkosky cut that are collinear or soft, identifying them with one object which has the sum of the momenta of the constituent particles. The dependence of the observable on each of the momenta of particles in the Cutkosky cut individually morphs into a unique dependence on the sum of the momenta of collinear or soft particles in the neighbourhood of a pinched point -37 -and the usual independence on individual momenta of the particles that are not collinear or soft: Ultimately, eq. (3.60) establishes whether an observable preserves the necessary conditions and can be used in defining a locally finite differential cross-section. The volume form σ d , with observables satisfying eq. (3.60) is thus locally finite on R 3L \ B isr, \ B (W). In order to conclude the proof, we study the soft scaling of σ d and determine under which conditions it has at most integrable singularities.

Soft scaling from the causal flow
We now study the integrability of σ d in a neighbourhood B (W). Since soft points are at the extrema of pinched singularities, the IR-safety of observables allows to write σ d as eq. (3.49) in a neighbourhood of the soft points. Furthermore, we see that in the limit of eq. (3.49) the remaining function has a scaling which, at worst, is the scaling of g or any of its derivatives in t, assuming h and O s are well behaved. Thus, in order to analyse the scaling of the integrand in the neighbourhood of soft singularities, we have to understand the scaling of g and its derivatives in t. We now identify extra constraints on the causal flow such that derivatives of g have the same scaling as g, assuming the observables and the h function are bounded on such neighbourhood.
We recall that g can be defined in terms of the manifestly causal (cLTD) representation, and particularly in terms of the simplified constrained form of eq. (3.44), such that (3.61) Let us start by observing that, by definition, g cannot be singular on any s-channel E-surface, as g is bounded on R 3L \ B isr, \ I \ B (W) (and excluding the UV region too). Near soft points the scaling of g is entirely determined by the inverse E-surfaces in E int and by the inverse energies (since g is not singular on any s-channel E-surface and since ISR singularities are excluded, thus leaving internal E-surfaces only).
We can now state a sufficient condition on the causal flow φ and the causal vector field κ for derivatives of g with respect to t to have the same (or less singular) scaling as g: if the deformation field satisfies the continuity constraint of eq. (3.37) of ref. [93]: This condition is obtained by iteratively applying Leibniz's rule on the n-th derivative of g, thus isolating products of multiple derivatives of inverse energies, of inverse internal E-surfaces and the numerator. Then, we require derivatives of inverse energies and of inverse internal E-surfaces to have the same or better soft scaling than inverse energies and inverse E-surfaces themselves. Observe

(3.64)
If the left-hand side and right-hand side of each of the two equation are required to have the same soft scaling, Q i ( κ) is required to vanish as fast or faster than the on-shell energy E i . In the construction of the deformation field of ref. [89] this constraint already appears and is called the continuity constraint. The role of this constraint in that work was to make sure that the deformation is continuous and does not cross a branch cut.
Since it is shown in ref. [89] that such causal vector field can be constructed, we conclude that thanks to eq. (3.62) the scaling of σ d is the same as the worst scaling between g and its derivatives. We are now able to show that σ d is integrable on a neighbourhood of soft points B (W). Because of the previous work, this is equivalent to showing that the scaling of g( k) at soft points only leads to integrable singularities. These findings can be summarised in the following claim: The scaling of g, on the other hand, follows from its explicit expression in eq. (3.61) and from the property in eq. (3.44). In order to establish that the LU representation is integrable on soft singularities, one must still show that the soft scaling of g, and therefore of σ d as well in virtue of eq. (3.65), is sufficiently tame for physical theories. This is the object of the next section.

Power-counting of soft singularities
In the following, we will analyse the scaling of σ d close to points at which connected clusters of particles simultaneously become soft (any configuration that is not the union of disconnected clusters can easily be shown to have a better scaling). Consider a connected subgraph of Γ, (s, e s = δ • (s)∪δ(s)) closed under momentum-conservation conditions (that is, if the momenta of edges in the subgraph uniquely determine any other edge's momentum, then that edge is also in the subgraph). Using eq. (3.61) and eq. (3.65), we can conclude that the worst possible scaling of σ d near a soft singularity can be at most σ d ∼ g ∼ δ −(|s|+|es|) , since each inverse energy corresponding to an edge in e s contributes one power and, in the worst-case scenario, there exists a term in the sum of eq. (3.44) which has |s| E-surfaces of E int vanishing (such a case would correspond to a term of eq. (3.44) whose respective cross-free family contains a cross-free family of subsets of (s, e s )).
Since the actual soft scaling of the numerator is theory-dependent, we expect integrability at soft points to be conditional on the specific theory one is considering. We will now simplify the discussion by proving integrability for the physical theory of Yang-Mills coupled with massless fermions theories. Since the Higgs boson is massive, including it does not change the reasoning, along with any type of massive particles. We can distinguish between fermion edges and vector boson edges by writing e s = e f ∪e b . Furthermore, we write s  for vertices at which three vector bosons meet, s (4) for vertices at which four vector bosons meet, so that s = s (4) . Finally, if all momenta in the subgraph are set to zero, the degree of divergence associated to this subgraph is d = 3L s − |e s | − |s| + |s (3.66) We substitute the following relations between the number of vertices, their degree and the edges of a graph, |e| = 3/2|s (3) f b | + |δ(s)|/2 and |e f | = |s (3) f b |, such that we obtain: which confirms the inexistence of non-physical soft singularities, and confirms that the powercounting at soft points, in the absence of pinched E-surfaces (which are subject to cancellation), yields integrable singularities. However, for a cubic scalar theory, the power-counting formula for soft clusters reads d = |δ(s)| − |s| (3.68) which can become negative. This reproduces the known fact [124] that the IR region of superrenormalisable theories such as cubic scalar theory exhibits soft divergences that do not cancel. As an example of a graph that is still divergent within LU, we give an example of a supergraph in cubic scalar theory that features a quadratic soft divergence (d = −2) on the grey lines: This concludes the proof, as we have now shown that σ d is integrable on R 3L \ B isr, for a physical theory.

Self-energies and IR cancellations in non-abelian gauge theories
External self-energy corrections, and specifically one particle irreducible insertions on edges that belong to a Cutkosky cut lead to the presence of spurious divergences that can be eliminated through careful on-shell renormalisation of the propagators. Furthermore, the nature of the IR cancellation pattern requires that the computation be done with non-truncated amplitudes.
The problem of self-energy corrections also makes very clear another issue in preserving IR cancellations: external physical boson propagators correspond to a different Feynman rule than internal boson propagators. The difference in tensor structure also leads to miscancellations. In order to address this issue, we allow for ghosts and unphysical bosons to be external particles.

Propagator renormalisation and IR cancellations
When a supergraph features raised propagators, which in physical theories appear as a result of a self-energy insertion, the naive substitution of propagators with Dirac deltas leads to manifestly ill-defined interference diagrams, as performing the substitution for one of the raised propagators exactly evaluates the remaining repeated propagators on their mass shell.
The first possible solution to this problem relies on recognising the origin of the Cutkosky rule as an analogue for performing contour integration. This suggests that one should use the residue formula for higher-order poles, and modify the Cutkosky rule accordingly so that each propagator, raised to the power n, is substituted by an appropriately normalised (n−1)-th derivative of Dirac's delta function. However, considering this higher-order residue formula would require taking derivatives of the observable function, which is not always possible. As we discuss in this section, it is instead possible to apply the derivative to the self-energy subgraph only.
In ref. [79], D. Soper circumvents the problem of taking derivatives for addressing raised propagators by engineering an alternative representation of self-energy correction. Through the use of a dispersion relation and algebraic manipulation of the numerator, he was able to show that the raised propagators can be re-absorbed into a propagator with power one. However, the extension of this procedure to more complicated topologies is challenging, especially when the self-energy diagram features UV divergences, since the alternate representation does not allow to construct UV counterterms in four dimensions, but requires to derive them directly from the three-dimensional representation.
In this section we derive a novel treatment of self-energy diagrams that operates in Minkowsi space and does not require taking derivatives of observables. As the raised propagator issue arises when one of the propagators adjacent to a self-energy corrections is cut (that is, it is set to be an external particle), it is clear that any solution should bear some relation to the LSZ formalism, Dyson resummation and propagator renormalization. We will show that raised propagators disappear after mass renormalization in the on-shell scheme.
Let us start by defining the renormalized 1PI graph as In the following, we consider Σ to be derived from a scalar theory. The extension to physical theories requires extra attention and is discussed in sect. 4.2. The extra issues due to considering physical theories are however related to the role of numerator algebra and gauge invariance in the IR cancellation mechanism and do not directly affect the general treatment of raised propagators hereby presented.
-41 -Consider now the Taylor expansion of the renormalized one-particle irreducible function The zero-th order of the expansion contains the raised propagator. By imposing the first renormalization condition of the on-shell scheme, Σ R (m 2 ) = 0, we obtain that Thus, each self-energy correction of the graph only brings one power of the external particle propagator. In practice, this allows us to perform Cutkosky cuts by using the rule associated to simple poles. We observe that all terms O (p 2 − m 2 ) 0 are not singular at p 2 = m 2 , and therefore their associated residue is zero. The modified Cutkosky rule, when applied to a propagator that has a self-energy insertion then reads: It is important to observe that performing renormalization selectively for the external self-energy corrections does not break the pattern of IR cancellations, since is free of IR singularities. This statement implies that the IR singularities of two-point functions are fully contained in its derivative when evaluated on the on-shell hyperboloid, and that its singular structure is equal to that of the sum over all possible Cutkosky cuts of the two-point function. One can also think of eq. (4.5) as establishing a counterterm for the self-energy correction. However, in doing so, the relation with on-shell renormalization is lost.
Finally, we can discuss how field renormalization participates in the clear separation of UV and IR divergences. The poles obtained through the dimensional regularization of the self-energy insertions are usually renormalized, indifferently of their origin (UV or IR), and reabsorbed in the coupling constant. Considering truncated amplitudes and renormalizing away IR singularities of external self-energies obscures the presence of a deeper cancellation mechanism for IR divergences. In the on-shell scheme, required for eq. (4.3) to hold, the two-point function is substituted with a local representation of its pure IR pole that eventually cancels with its real-emission counterpart, while the finite part is set to zero, as prescribed by the renormalization condition. We start with the on-shell scheme constraint on the finite part of the energy derivative: Eq. (4.6) implies that dΣ R dp0 must be a local representation of the pure IR pole (in dimensional regularisation). This leads to the following formula: dΣ R dp 0 = dΣ dp 0 − CT uv dΣ dp 0 d=4 − finite dΣ dp 0 − finite CT uv dΣ dp 0 d .

(4.7)
dΣ R /dp 0 is finite in the UV region. It is divergent in the IR region, so that its expression in dimensional regularization features IR poles. Finally, it has zero finite part. Because of eq. (4.5), its IR poles locally cancel with the real counterparts. Figure 11. On the left, the derivative of the two point function with cuts on both of the propagators adjacent to it, obtained from the modified Cutkosky rule (4.4). It has the same IR structure as the sum of contributions which can be written as the diagram on the right.

Preserving IR cancellations in non-abelian gauge theories
In the derivations of sect. 2 and sect. 3 we have assumed that the procedure of performing Cutkosky cuts on the supergraph does not affect the numerator, i.e., that the numerator of an interference diagram is the same as that of the supergraph with some energies set to their on-shell value. This property is fundamental for showing local IR cancellations, and was used in the heuristic argument in sect. 2.5 (where the theory was taken to be scalar) and in eq. (3.38). For gauge theories, the numerator may seem affected by performing Cutkosky cuts, since there is a distinction between physical photons, which are allowed to be asymptotic states, and unphysical photons, which are not. In the covariant gauge, the functional form of polarisation sums for spinone bosons differ from that of propagators for spin-one bosons. This leads to apparent divergences at the local level associated with the non-physical degrees of freedom. Ward identities can be used to regularise such divergences, but this makes the proof of cancellation and the construction of interference diagrams more laborious.
Instead, we use the numerator of the usual Feynman propagators instead of polarisation sums for cut gluons, and allow ghosts to appear as final-state particles. This procedure relies on a known alternative to assigning a different Feynman rule to external bosons and internal bosons, and uses the fact that asymptotic states are defined up to states of the form Q|ψ , where Q is the BRST operator. The vector Q|ψ can then be written as a sum of vectors representing particles including non-physical bosons and ghosts. Rewriting transition amplitudes in this fashion, the polarisation sum rule for a boson can be substituted by the rule for a cut Feynman propagator, at the cost of including all the discontinuities of all supergraphs including those cutting through ghosts. We explicitly tested for all supergraphs and each LO Cutkosky cut of complicated processes such as γ → ttggg that the LU integrand is locally identical when computed with final-state ghosts instead of physical external gluons. We note that a symmetrisation over the final-state gluon momenta is necessary in order to carry out this exercise.
The sum over discontinuities of a fixed supergraph, although not gauge invariant, can be made locally IR finite through the procedure of sect. 4.1. This sum includes thresholds corresponding to ghosts becoming on-shell. We now outline a procedure that determines the contributions necessary to reproduce gauge-invariant quantities: (1) generate all supergraphs with the chosen theory's particle content, including ghosts and Goldstone bosons if needed, (2) proceed to list all Cutkosky cuts of such supergraph, including those cutting through ghosts and Goldstone bosons and (3) express this sum locally, as done in eq. (3.38), by aligning the loop momentum routing consistently for interference diagrams originating from the same supergraph. We stress that in this framework the action of a Cutkosky cut is equivalent to that of substituting the denominator of the propagator by two times its on-shell energy, and evaluating the numerator at the on-shell value. No polarisation sum rule is ever needed, and thus the interference diagram has the same numerator structure as the supergraph.
Finally, we observe that calculations in the Feynman gauge, i.e, with the gauge parameter set to one, are naturally the most convenient in our formalism, since the gluon propagators will not contain a term proportional to p µ p ν /p 4 , which leads to raised propagators whose treatment is difficult to generalise and can lead to ad hoc manipulation of the local integrand. In the Feynman gauge, this means that we treat gluon propagators, both cut and uncut, as a scalar theory propagator multiplied with the tensor structure g µν .
The procedure outlined above is easily automated and shows that local IR cancellations become straightforward once physical and non-physical degrees of freedom, whether they belong to virtual or real contributions, are treated in a fully symmetrical manner.
- 44 -In this section, we discuss future directions for extending the scope of Local Unitarity and we detail the path towards its fully automated implementation for the numerical computation of higher fixedorder corrections to generic differential cross sections.

First steps towards Local Unitarity applicable to initial-state singularities
The general route for extending the mechanism of FSR cancellations to ISR involves generalising eq. (3.38) so that it includes Cutkosky cuts corresponding to initial-state thresholds. The guiding principle is that any extension of σ d should still be written as a weighted sum over the discontinuities of the supergraph that is locally IR-finite.
A first attempt at the solution is to consider extending the sum of eq. (3.38) to the initialstate thresholds. In doing so, we can keep our representation of supergraphs involving a fixed set of initial states, but then consider all possible thresholds, including those in E isr . Although this may provide a valid cancellation mechanism, it yields disconnected amplitudes and imposes a fixed initial-state multiplicity which is at odds with the notion of backward evolution of initial-state partons. We instead tentatively consider a paradigm in which degenerate initial state configurations are accounted for by kinematic configurations featuring multiple incoming asymptotic states. This corresponds to thresholds identified by N in (resp. N out ) incoming (resp. outgoing) momenta. So far, our construction of Local Unitarity considered a constant multiplicity N in = 2 or N in = 1 and we discuss here the first steps towards its generalisation for arbitrary (and varying) N in . The implicit equation of such generalised thresholds reads: for the multiplicity N in and N out relevant to the perturbative order considered. In order to give rise to such threshold, the N in external legs of our current supergraph construction must be closed onto themselves so as to form a vacuum diagram. A first technical difficulty in addressing vacuum supergraphs is that their LTD representation only features internal E-surfaces, that is E int = E, so that all thresholds of the form of eq. (5.1) are subject to dual cancellations. The origin of dual cancellation stems from the fact that all propagators of the original representation of the vacuum supergraph in Minkowski space are assign the identical causal prescription +i with a positive sign. For this reason, we hypothesise that the generalisation of LU for initial-state singularities requires a modification of the sign of the causal prescription of some of the original loop propagators. We therefore propose a candidate for an alternative definition of vacuum supergraphs that is locally finite for all singularities and exhibits degeneracy in the multiplicity of both initial and final states. This alternative expression M G [v, s] for the vacuum supergraph G = (v, e) (with set of vertices v and oriented edges e) depends on two new quantities: a set of signs v and a set of vertices s. It reads: where s ⊂ v and δ • (s) is a connected subgraph. We stress that it is important to properly account for the symmetry factor of G (part of the numerator N in eq. (5.2)) in order to avoid doublecounting. Furthermore we defined: -45 -which corresponds to the subset of edges within s and at its boundary when satisfying a particular choice for the sign of the energy flow across it. The sign of the causal prescriptions for the propagators in M G [v, s] are engineered such that after applying LTD, it contains (among other ones) a threshold corresponding to the particles in e + v being on-shell with positive energy and the particles in e − v being on-shell with negative energy. These two sets are then naturally identified with the incoming and outgoing particles respectively, and participate in a threshold of the form of eq. (5.1). We note that the choice of i -prescription for the edges not in δ(s) is set so as to reproduce complex conjugation (that is an opposite prescriptions for E-surfaces on either side of the Cutkosky cut).
We now construct a set containing all possible connected subgraphs and all physical thresholds associated to such subgraph: so that each element (s, v) ∈ H can be associated to a threshold defined by the zeros of the equation We show in tab. 1 the list of all twenty thresholds in H contributing to the process dd → Z for one vacuum supergraph G relevant to that process. The only missing component in order to construct an analogue of σ d that also contains degenerate initial-state configurations is the corresponding causal flow. It should be constructed as the solution of an ODE whose defining vector field is causal on all surfaces in H: When the energy flow vector v does not have all of its components set positive or negative, the set of requirements above is different than the case studied in ref. [89]. We will however not discuss here whether eq. (5.7) can be satisfied by a vector field κ; this investigation will be left as future work. Instead, we assume here that such a field exists and that it is thus possible to map any point of R 3L to points on distinct thresholds. We then construct the extension of eq. (3.38) by summing the discontinuity acquired by contour deforming around the threshold η (s,v) of M G [v, s] for all possible elements of (s, v) ∈ H. Along the causal flow one can expand the thresholds around the solutions η (s,v) (φ(t s,v , k)) = 0 and carry out a similar study as the one performed in sect. 3 and thus arrive at the same conclusion regarding the pairwise cancellation of s-channel E-surfaces. We observe that the residue of M G [v, s] at η (s,v) factorises into two ordinary amplitude graphs δ • (s) and δ • (v \ s) times a product of inverse energies for each element of δ(s). Thus, using insights from the factorisation property noted in sect. 2.5, we see that the candidate expression for vacuum supergraphs of eq. (5.2) reproduces the known structure of cross sections for hadronic observables. Finally, we can define a local representation of the differential cross section stemming from the vacuum supergraph G that is locally finite at locations of both initial-state and final-state pinches (but excluding t-channel singularities [122]): (5.8) Table 1. List of all twenty thresholds (s, v) in H contributing to dd → Z stemming from one of the two distinct supergraphs G relevant to this process. The vertices in black are part of the set s while the grey ones are part of v \ s. The edges ei coloured in black are in e + v,s , that is either in δ • (s) (non complex-conjugated internal propagators) or in δ(s) with ve i = +1 (final-states). Similarly, the edges ej coloured in grey are in e − v,v\s , that is either in δ • (v \ s) (complex-conjugated internal propagators) or in δ(v \ s) with ve j = −1 (initial-states). We stress that we only retain contributions for which the final state includes a Z-boson, as per the process definition. As for FSR supergraphs, the potential double-counting of interference diagrams is compensated for by the supergraph symmetry factor. The embedding is chosen that the central vertex of the vacuum supergraph is always placed to the right of the Cutkosky cut, so that when it is black (i.e. part of s) the resulting interference diagram needs to be flipped around the Cutkosky cut for complex conjugated propagators to sit on its right.
where H k,φ is the analogue of E k,φ for H defined in eq. (3.34). There is a crucial difference between the equation above and its final-state counterpart given in eq. (3.38) in which we include all the thresholds that appear in the LTD representation of a supergraph specified with a fixed i prescription. Instead, the expression of eq. (5.8) considers the sum over a set of different prescriptions assigned to the same supergraph due to the nature of the thresholds singularities of vacuum diagrams. Moreover, for each choice (s, v), only the residue of the corresponding threshold is selected (i.e. Cutkosky cut c (s,v) ).
The final expression of eq. (5.8) also introduces the new symbol O s,v which is the observable associated to the vacuum cut s for the chosen energy flow v. This observable can be set to zero for many of the (s,v) cuts that does not match the process of interest (e.g. all 1 → N out configurations in the case of hadronic scatterings). Even in the context of our original construction of Local Unitarity, we can already find observable functions depending on the details of the Cutkosky cut.
For example when computing the cross section for the process γ → Htt, the observable function would select out all Cutkosky cuts involving only the top and anti-top quarks. In the case of hadronic scatterings with initial-state singularities, the role of the observable function O s,v is even more important. Indeed, it should also encode the details of the structure of the constituents of the colliding incoming bound states, that is parton distribution functions (PDFs). In particular, it is responsible to confine the contributions from (possibly degenerate) incoming momenta to lie within the close vicinity of an (anti-collinear) back-to-back configuration. In order to achieve this, the "initial-state observable" needs to contribute to the LU integrand with an actual weight (as opposed to just the common Heaviside) so as to be able to input the initial states density. However, we note that in the framework proposed in this section, the treatment of parton distribution functions would likely be different than that of the traditional collinear factorisation paradigm [94][95][96]. Indeed, since σ d is free of all s-channel infrared singularities, there is no direct analogue to the traditional PDF counterterms nor to the common factorisation scale µ F driving the renormalisation group equations induced by the singularities of the PDF evolution kernels. Instead, we find an alternative perspective wherein initial states can vary in multiplicity, in the exact same way as final states do, and whose degeneracy scale is set dynamically from the requirement of an IR-safe definition of the initial-state observable O s,v .
This concludes the discussion of our preliminary investigation of the definition of a differential cross section that is locally finite also on s-channel initial-state thresholds We leave further developments of this construction to a future publication, including the regularisation of the remaining t-channel singularities and a quantitative investigation of its relation to traditional collinear factorisation and PDFs.

UV counterterms and renormalisation
The regularisation of the UV behaviour of the LU representation, together with the introduction of renormalisation contributions, lead to the following schematic modification of eq. (2.21) for the differential cross section: Γ,s,u are their counterparts integrated analytically using dimensional regularisation. Finally, the operator R applies to the selected set of higher-order correction supergraphs G and generates the set of all corresponding renormalisation supergraphs. For each of them, we can identify a renormalisation counterterm δZ Γ R factorising the LU representation of a lower-order supergraph. In this section, we discuss the construction of each ingredient of eq. (5.9) qualitatively.
By construction, interference diagrams in physical theories must have a UV singular structure that corresponds to the combined UV singularities of the two subgraph amplitudes to the left and right of the Cutkosky cut. We choose to construct the local UV counterterms I (UV) in Minkowski space, and then convert them to the (c)LTD representation. This allows for the use of traditional analytic techniques and dimensional regularisation when computing their integrated counterpart I (UV) .
Beyond NLO, the local UV regularisation must be performed carefully so that the remaining UV divergences in one-loop UV counterterms locally cancel against the two-loop UV counterterms. In other words, the overlapping UV divergences in Γ,s,u must cancel locally amongst the terms in this sum. We achieve this by implementing S (x) Γ,s and I (UV) Γ,s,u according to the BPHZ forest formula treatment [125][126][127]. This procedure involves identifying all possible UV singular subgraphs and constructing an appropriate approximant of the integrand on UV limits. Proper subtraction of the original integrand and, specifically, the treatment of overlapping subdivergences is guaranteed by the forest formula. In particular, the study of overlapping subdivergences requires defining spinneys, which are collections of disjoint UV divergent subgraphs. We stress that this entire procedure must be repeated iteratively when constructing the terms I (UV) Γ,s,u and I Γ R ,s so as to also locally subtract their remaining UV divergences (beyond one-loop).
The BPHZ forest formula is independent on the chosen renormalization scheme and only prescribes a Taylor expansion of the denominators of UV subgraphs in their external momenta. Whereas one traditionally also considers expanding the numerator around the UV point, we find it more convenient to expand propagator denominators only. This operation is local, since it relies on a Taylor expansion of the four-dimensional integrand of the amplitude. It follows that the integrand can also be subtracted locally.
We identify the UV subgraphs with scalar graphs and assign them their superficial degree of divergence. After Taylor expanding in its external momenta, the scalar graph can be represented as a power series in scalar vacuum diagrams. We assign to each vacuum diagram propagator a UV mass m UV . This UV counterterm is added to the integration and subtracts correctly the UV singularity represented by the subgraph. The analytically integrated counterpart I Γ,s,u , which always corresponds the a massive vacuum diagram integral that can easily be computed using integration-by-parts identities and a finite set of known master integrals. In order to obtain a scalar expression, we tensor reduce the UV subgraph so that it completely factorises. Finally, we multiply the finite part of the integrated UV counterterm into the remaining diagram.
Ultimately, the integrated counterterm will be a Laurent series in = 4−d. These UV poles are cancelled by the renormalization supergraphs Γ R listed in R(G) and constructed by contracting the disjoint UV subgraphs in the spinneys to vertices. This forms separate renormalisation supergraphs with lower loop counts, and we assign to each vertex obtained in this fashion their respective renormalization constant. We stress that, beyond NLO, it is important that the poles in the dimensional regulator = (d − 4)/2 of the renormalization components δZ Γ R multiply the O( ) parts of I Γ R ,s computed in d-dimensions. This dynamic generation R(G) from the list of selected supergraphs considered G ensures the consistency of the renormalization procedure. We can numerically test if the UV poles of the renormalisation graphs cancel with the ones from the integrated counterterm by numerically integrating the numerator dimensional regulator pole coefficients instead of its finite part. We stress that this pole cancellation check can in principle also be made local by consistently rerouting the momenta of all supergraphs in {Γ ∈ G | Γ R ∈ R({Γ}) contributing to a particular renormalisation supergraph Γ R . In the example of the UV regularisation of the double-triangle supergraph that will be investigated in sect. 6.1.1, we find a one-to-one correspondence between the supergraph and its renormalisation contribution, leading to a direct cancellation of the UV poles (see eq. (6.11)). Because this is in general not the case, we opt to account separately for higher order supergraphs and their renormalisation terms.
Often, an approach similar to the one by D. Soper [79] is chosen in order to construct a modified local countertermĨ Γ,s,u which combines both local and integrated UV counterterms under the same integral measure. This can be done by adjusting the higher-order terms of the original local UV counterterm of eq. (6.5) obtained from the strict UV limit. In this way,Ĩ (UV) Γ,s,u directly integrates into a finite part that is zero and automatically reproduces MS results, or alternatively allows one to add whichever finite part from the renormalisation constant of the preferred scheme without having to worry about the impact of the UV regularisation in Local Unitarity. Such an adjustment of the numerator is however typically not systematic and therefore not practical for automation. Instead, we multiply the integrated counterterm I (UV) Γ,s,u with a rational polynomial in the loop degrees of freedom that are analytically integrated over (e.g. a simple product of massive one-loop tadpoles) so that it can be systematically combined with I (UV) Γ,s,u . Since we separately account for the contribution of the integrated UV counterterms and renormalisation constants, we must pay particular attention to the particular dimensional regularisation scheme considered when computing both terms. We identify four different parts in the expression of the integrated UV counterterms: the numerator N UV stemming from the edges part of the UV subdiagram, the external factorised ones N E and similarly for the denominators, D UV and D E . As is well-known in the literature [111,[128][129][130][131], different choices regarding the dimensionality adopted for each of these pieces are possible and ultimately equivalent as long as they are performed consistently across terms whose poles cancel. We find it most convenient to build an automated procedure by treating the combination of N UV , N E and D UV in d-dimension and D E in four dimensions. This choice also implies that no particular attention needs to be paid to the rational parts [132][133][134][135][136] of the integral as they will then automatically be accounted for in this manner.

Automation and numerical efficiency
The relevance of the LU representation for collider phenomenology crucially depends on the performance of its implementation. Our implementation leverages the framework of MadGraph5_aMC@NLO [14] (MG5aMC henceforth) for the abstract representation of the user inputs such as the observable definition, process generation syntax and physics model considered. In particular, the Feynman rules are provided through a model file following the Universal Feynrules Output (UFO) conventions [137]. This choice is motivated by our longer-term goal of offering the user an automated environment similar to MG5aMC for steering its simulations. However, the very different nature of LU requires an independent program. We therefore decomposed the implementation of the LU representation into individual fundamental tasks that are efficiently carried out by the development of dedicated codes that we discuss in this section.

Process generation
Given the specification of a scattering process and UFO model using MG5aMC syntax, we steer QGRAF [138] to generate all contributing supergraphs and output them in a Python-readable format. Next, all isomorphic supergraphs are combined and for each resulting supergraph we enumerate all Cutkosky cuts whose final state particles are compatible with the process definition. We stress that it is important to also weight each supergraph with its symmetry factor provided -50 -by QGRAF as it corrects for the over-counting due to identical particles (i.e. it plays a similar role as the customary normalisation of amplitudes by the symmetry factor stemming from identical particles in their final states). We then identify external one-particle irreducible diagrams in the supergraphs (that is, all the 1PIs that are adjacent to a Cutkosky cut), and apply the propagator treatment of sect. 4.1. We also detect all UV divergent subgraphs and construct local and integrated UV counterterms for them, together with renormalisation contributions, by following the procedure laid out in sect. 5.2. Finally, we collect all information on the supergraph that must be computed into input files ready to be processed and optimised by FORM into C compiled libraries.

Graph evaluation
The treatment of the numerator is often the bottleneck for both generation and run-time efficiency when considering higher-loop Feynman diagrams. An efficient implementation in FORM [139] is therefore crucial. Our FORM code first substitutes the Feynman rules for each supergraph, contracts all indices, and applies polarisation sum rules. After simplifying all Dirac traces and colour algebra, the output of this procedure is a sum of scalar products of loop momenta and external momenta.
Next, a Python code constructs the general Loop-Tree Duality expression derived in ref. [88] and the Manifestly Causal LTD (cLTD) expression derived in ref. [93] in a FORM-readable format. This expression is added to the expression of the numerator. At this stage, the entire integrand per Cutskosky cut is represented in a FORM expression. Next, FORM performs all derivatives necessary for the propagator treatment of sect. 4.1 as well as for the UV counterterms. This integrand can be viewed as a polynomial in the scalar products and linear propagators. We use FORM's optimisation [140,141] algorithms to generate highly efficient C code that can evaluate the entire integrand per Cutkosky cut by finding an optimal Horner scheme for the resulting polynomial and grouping common subexpressions. This optimisation step typically reduces the number of arithmetic operations by several orders of magnitude. Additionally, the C compiler further improves the implementation by vectorising (parts of) the resulting C code. This finally results in evaluation times typically below 100 µs, even for complicated processes. For example, one particular supergraph with 11 Cutkosky cuts contributing to the NNLO correction to the process H → ttgg can be evaluated for one sample point in only 25 µs (without multi-channeling). This means that the evaluation time of the integrand is not a limiting factor.

Numerical stability
The LU integrand consists of terms that can feature large numerical cancellations, leading to instabilities when considering finite-size arithmetic. Such unstable points can be large in magnitude (see, e.g. fig. 2 of ref. [93] for numerical stabilities in the UV) and spoil the central value estimator of the Monte-Carlo integration. It is therefore necessary to monitor the numerical stability of the integrand evaluation for each sample point and properly address cases that are insufficiently precise. In order to estimate the numerical accuracy of the integrand for one particular sample point, we repeat its evaluation with all spatial directions subject to arbitrary rotations and compare the results obtained in this way, as they are analytically identical in virtue of the symmetries of the integrand. The original sample point is then deemed unstable if its rotated evaluations differ by more than a certain threshold. Since the LTD integrands may exhibit instabilities in the UV that are removed by cLTD (see ref. [93]), we perform the same procedure using cLTD instead of LTD. If the point is still considered unstable, we repeat this procedure using quadruple-precision arithmetic which is -51 -about 100 times slower but rescues the majority of unstable points. Exceptional sample points that are still deemed unstable at this stage are typically very close to infrared singularities or cancelling E-surfaces and we override the evaluation to zero. By varying the numerical stability requirements, we can verify that our procedure does not alter the result of the Monte-Carlo integration.

Contour deformation
When considering non-trivial observables or scattering process definitions that exclude certain Cutkosky cuts, the LU integrand may still be singular at the location of non-pinched thresholds. For these cases a complex contour deformation that follows the construction of ref. [89] is required. In the context of the LU representation, we deform the loop variables of the amplitude on the left and right side of the Cutkosky cut independently, to ensure that the observable function is evaluated with real-valued kinematics (see sect. 2.4). On pinched surfaces the deformation must be zero in order not to disrupt the cancellation of real and virtual Cutkosky cuts. We further dampen the deformation magnitude on pinched E-surfaces by multiplying the deformation field with the following dampening factor d IR : for each pinched E-surface with implicit equation η IR , in addition to the dampening enforced by the complex pole constraint (see sect. 3.3.2 of ref. [89]). The construction of ref. [89] requires a set of valid deformation sources to be constructed. This problem is NP-hard and requires finding points in the internal region of overlapping E-surfaces. The difference between determining the overlap structure for an amplitude with fixed external momenta and determining it for the amplitudes participating in a supergraph is that in the latter case, the external momenta of the loop subgraphs change with every Monte Carlo sample point. Consequently, the overlap structure of the E-surfaces must be recomputed at run-time for each sample point and Cutkosky cut. Thus, it is important that the overlap structure can be determined efficiently. Testing for overlap of E-surfaces in a set E can be achieved by finding a point in the interior of all E-surfaces in E. Such problems are known as a second-order cone program in the field of convex optimisation. We use the ECOS solver [142] to determine such points. ECOS takes about 0.5 ms for a test, and in the worst case 2 |E| tests are required. In order to reduce the number of necessary overlap tests, we only test groups of E-surfaces if each pair in the group pairwise overlaps, which can be determined using fast heuristics. After these optimisations, even for complicated cases such as the two-loop amplitude 2L6P.e K1 from ref. [89], which features 21 E-surfaces and 34 unique sources, the construction time is only about 18 milliseconds. We stress that the loop topologies most relevant to LU often involve many on-shell massless external momenta, for which the number of nonpinched E-surfaces is generally small since many E-surfaces are pinched in this case. Additionally, we only have to contour deform around E-surfaces whose pair-wise cancellation is spoiled by either the observable or process definition. We thus find that for cases of phenomenological relevance the determination of the overlap structure, and therefore the complete deformation, typically requires about one millisecond per sample point.

Phase space sampling and variance reduction
Even once a local representation of a differential cross section free of both IR and UV singularities is established, it remains to be shown that its numerical integration is feasible in practice. In this subsection, we discuss the main improvements that can help reduce the variance of the integrand.

Multi-channelling
Our computation of the differential cross section in Local Unitarity is organised in three nested levels, as made clear by eq. (2.21) which we rewrite here with a notation appropriate to this section: which features: • A discrete sum over each supergraph Γ ∈ G of the scattering process of interest.
• A discrete sum over each Cutkosky cut s Γ ∈ E Γ s-ch appearing in supergraph Γ ∈ G.
• An integral over the spatial loop degrees of freedom of the integrand I sΓ ( k), which contains all elements of the LU representation, in particular the observable function O and the solution t sΓ ( k) of the causal flow enforcing the conservation of on-shell energies involved in the Cutkosky cut c sΓ .
In order to improve the Monte-Carlo accuracy, one must reduce the variance of the estimator for σ O . This is achieved by using a discrete importance sampling technique for the sum of contributions and an educated choice for the parameterisation of the continuous degrees of freedom, whose Jacobian should approximate the (inverse of) integrand as closely as possible. We discuss both aspects in turn.
Optimising the sampling of the sum over supergraphs is straightforward, since each contribution can in principle be integrated completely separately. This aspect is unique to the LU representation, and so are the benefits discussed below stemming from a discrete importance sampling over supergraphs. The optimal discrete probability distribution of sample points over supergraphs is obtained by following the inverse of their relative contribution to the total cross section. One typically starts with a uniform discrete grid that is progressively updated according to the new estimates of the relative cross sections for each supergraph. While simple to implement, this sampling optimisation can be very powerful if the supergraph contributions are unevenly distributed. In fig. 12a, we show their distribution for the LO cross section from the scattering processes e + e − → γ → ttggg and e + e − → γ → ttg hḡh g where g h denotes a QCD ghost. We find that the relative contributions from supergraphs spans more than four orders of magnitude. While this observation depends on the scattering process, observable and gauge considered, we expect this uneven distribution to be quite generic and thus adaptive sampling can significantly help to counterbalance the increase in computation time resulting from the growth in the number of distinct supergraphs for more complicated processes. Fig. 12a shows no apparent large (gauge) cancellations between supergraphs, showing that the supergraphs can be integrated independently. While there may still be benefits from integrating all supergraphs in a correlated manner, identifying a common parameterisation [143] that leverages potential (gauge) cancellations is difficult and unlikely to outweigh the gain from discrete adaptive sampling over supergraphs. One exception being supergraphs involving QCD ghosts. These supergraphs always share the same topology as that of their gluonic counterpart, and can thus be integrated together with a common loop momentum basis. However, in all LO cross sections we investigated thus far, we have never found large cancellations within such groups of supergraphs. The large spread in the distribution of supergraph contributions is perhaps less surprising when placed in the context of more general computations of mixed QCD and EW corrections in the presence of resonant intermediate unstable particles. In that case, one expects supergraphs encoding pure QCD corrections to dominate those corresponding to pure EW corrections. Similarly, supergraphs involving Breit-Wigner enhancements are relevant for kinematic regions where the unstable particle is close to its mass shell but typically contribute less than non-resonant topologies to the inclusive cross section. The separation of these various regimes within traditional methods is complicated, since only subsets of terms that are gauge invariant and that can be rendered finite within their computational paradigm can be used. Instead, LU offers the perspective of systematically generating all contributions and efficiently and dynamically home in onto the relevant ones for a particular observable of interest. This procedure is appealing both from a numerical and phenomenological perspective as it helps the identification of the underlying dynamics at play for any collider signature.
Once a particular supergraph Γ is selected, we must establish a strategy for sampling its Cutkosky cut contributions s Γ and continuous degrees of freedom. The LU representation relies on pairwise cancellations between the terms I sΓ so that they must be integrated together with a common loop momentum basis and parameterisation. At the same time, the LU representation still features integrable singularities on soft kinematic configurations and when internal propagators become on-shell, which we must render bounded with the Jacobian of an appropriate parameterisation so as to make the variance of the LU integrand finite. Multi-channeling is a common technique (see, e.g. refs. [10,80,89,144,145] for examples in a similar context) offering a solution to this problem by combining multiple parameterisations each designed to address the integrable singularities and enhancements of a single Cutkosky cut at a time. In its most generic form, the multi-channeling -54 -approach applied to eq. (5.11) for a L-loop supergraph yields: where x is the vector of integration variables defined in the unit hypercube, Ω Γ is the set of all integration channels chosen for the supergraph Γ, Φ a x := k a ( x) is the parameterisation chosen for channel a, and J a := det ∂ x ∂ Φ a is the Jacobian of said parameterisation. The scalar functions ξ a can be arbitrarily chosen 6 . We also define the inverse mapping Ψ a k := Φ a −1 ( k). The rewriting in eq. (5.12) shows that if one chooses the parameterisations Φ i and the functions ξ i such that: all integrable singularities are removed and the variance of the overall resulting integrand is minimised. In ref. [89], the authors chose to have one channel per loop momentum basis of the graph Γ, so that Ω Γ ≡ B Γ and the integration channels are then labelled with a basis b i which map to the scalar index i. Then, in order to accommodate eq. (5.14), ξ i was built from the product of all (complex) on-shell energies of the edges in b i whereas the spherical parameterisation adopted for the momenta in b i yielded a Jacobian J i that dampens all integrable singularities in ξ i . In this work and for all the numerical results presented in sect. 6, we adopted a different strategy. We also consider spherical parameterisations for each of the integration channel associated to one loop momentum basis of the supergraph Γ, but we choose ξ i ≡ J i Ψ i k which dampens all integrable singularities as per eq. (5.14) but also exactly fulfils eq. (5.13). We note that we combine integration channels from loop momenta bases that are degenerate due to our particular choice of Lorentz frame, for which Q = 0.
Our particular realisation of the generic multi-channeling strategy of eq. (5.11) then reads: The effect of this multi-channeling on the integrand I MC Σ is visualised in sect. 6.1.2. One can consider a discrete adaptive importance sampling similar to the one already discussed for the sum over supergraphs. However, we see in fig. 12b that the distribution of the cross section stemming from each of the 16 multi-channel integrands I a Γ of the 3-loop scalar supergraph a.3) of fig. 3 is very even. Even though the gain in variance reduction is more moderate in this case, sampling the integration channels instead of taking their sum yields a faster evaluation of each sample point and therefore a Monte-Carlo error reduction by a factor of up to |B Γ | for a fixed run-time. Furthermore, performing a separate adaptive sampling for the kinematic dependency of each integration channel allows the integrator to better adapt to their individual shape as opposed to when training directly on their sum.
The distribution of the cross section from each channel and the resulting variance of I a Γ heavily depends on the particular choice of parameterisations Φ a . It is clear that the spherical parameterisations considered in this work are far from the optimal choice and prior work from D. Soper [80] showcases the benefits of adopting more elaborate parameterisations. More specifically, our choice does not account for the details of the particular topology of the subgraphs on each side of a particular Cutkosky cut. It also ignores the sharp shape of the term h t sΓ ( k) , part of I sΓ , which strongly enhances kinematic configurations k close to the Cutkosky cut surface η sΓ . We now discuss a future enhancement that improves on these aspects. We start by factorising the dimensions sampled into three subsets of variables: where L left (resp. L right ) is the loop count of the subgraphs Γ left (resp. Γ right ) to the left (resp. right) of the Cutkosky cut c sΓ , whose number of edges is equal to |c sΓ | = L − L left − L right + 1. The complete parameterisation Φ a can then be constructed in the following three successive steps: • First, the external momenta q sΓ := { k e |e ∈ c sΓ } with k ∈ δη sΓ is generated so as to directly lie on the Cutkosky surface η sΓ . To this end, we take advantage of the traditional phase-space factorisation already used by most event generators. This iterative construction of the phasespace parameterisation is organised such that the invariant mass of repeated propagators that appear 7 in both Γ left and Γ right is directly aligned to one variable x ext,k which is sampled with a probability density flattening the corresponding resonant enhancement (see e.g. ref. [144]). For "interfering" propagators, there is instead no good ansatz for the density with which they should be sampled (and those topologies are often suppressed), so that we can consider sampling them uniformly instead (see e.g. ref. [146]).
• Second, the configuration q sΓ must be "upscaled" so as to cover the entire volume R 3(|cs Γ |−1) . This aspect is crucial to the LU construction as the causal flow induces a mapping onto the Cutkosky surface of each contribution and is key to obtaining local IR cancellations between them. The most natural choice for achieving this upscaling is to use the causal flow of the supergraph Γ itself, with boundary conditions given by q sΓ and a value of t sΓ derived from the input variable x t and sampled with a density of dt sΓ /dx t = 1/h(1/t) = h(t). This ensures that the chosen parametrisation renders the integrand I sΓ ( Φ sΓ x ) independent of the choice of normalising function h(t). The choice of the σ-tunable function h(t) = 2σ π cos π 2σ t σ 1+t 2σ is particularly convenient so as to be able to analytically solve for the cumulative distribution function H(t) = t 0 dt h(t ) = 2σ π(1+σ) cos π 2σ 2 F 1 1, σ+1 2σ , 1 2 3 + 1 σ , −t 2σ , which can then easily be inverted numerically.
The other Cutkosky cut integrands I s Γ ( Φ sΓ x ), with s Γ = s Γ , still retain their dependence on h(t) and will be conveniently suppressed away from the intersection δη sΓ ∩ δη s Γ where pairwise E-surface cancellations occur. The spread of the function h(t) must then be adjusted so as to find the optimal balance between making it narrow enough so as to efficiently focus each channel onto its defining Cutkosky surface and wide enough so as to maintain a good cancellation of threshold singularities in the neighbourhood of their intersection.
• Last, the loop spatial degrees of freedom must be assigned using the variables x left and x right .
Note that the external kinematics of these loop subgraphs are now entirely determined from the quantities q sΓ and t sΓ built in the previous two steps. The situation is thus completely analogous to the sampling problem of normal loop integrals in their LTD representation already studied in ref. [89]. A minimal solution is then to consider a separate channel for each loop momentum basis of the reduced subgraphs Γ left and Γ right which can be sampled with spherical parameterisations. This is sufficient to dampen all integrable singularities, though it may be beneficial to construct additional channels using elliptical parameterisations whose Jacobians can flatten possible enhancements in the vicinity of non-pinched threshold singularities (E-surfaces) which may not always cancel when considering differential observables.
The procedure described above yields complicated parameterisations Φ i x that are nonetheless invertible and whose overall multi-channeling factor x . This improvement is also a necessity when considering processes featuring narrow Breit-Wigner resonances. At LO, one could recover the same performance as that of traditional integration techniques since the h(t) normalising function no longer impacts the variance of the integrands. Finally, notice that the number of integration channels are no longer |B Γ |, but instead |E Γ s-ch | × |B Γ left | × |B Γ right | (when ignoring the possibility of additional E-surface parameterisations for the sampling of Γ left and Γ right ), such that the sampling is automatically simplified for observables selecting only a subset of the possible Cutkosky cuts.

Advanced adaptive sampling
The techniques discussed in sect. 5.4.1 aim at leveraging our prior knowledge about the structure of the LU representation and the location of its enhancement so as to construct an integrand with the smallest possible variance. However, once we start sampling the integrand, we gain additional knowledge about its more detailed structure for the particular process and observable at hand. It is then possible to take advantage of this newly gained information to further reduce the integrand variance at run-time using machine-learning techniques implementing some form of iterative adaptive sampling which we discuss here.
The most straightforward way of improving on the parameterisation of the spatial loop degrees of freedom is to further refine it with a trainable integration grid assuming a factorised ansatz for the functional form of the integrand [147,148], which can be complemented with the use of coarse hypercells [149]. Our use case of LU requires considering an individual continuous grid for each possible choice of supergraph and integration channel, which are moreover subject to discrete importance sampling. We developed a new integrator for this purpose, dubbed Havana, that implements this hybrid adaptive sampling through a system of continuous grids nested within multiple layers of discrete ones. This will serve as a basis to accommodate further refinements and help to design custom parallelisation strategies in order to facilitate deployment on various computing architectures.
Discrete importance sampling over integration channels can also be used to improve upon the quality of the integrand approximation obtained from the multi-channeling factor of eq. (5.14). The a priori unknown relative magnitudes of the particular enhancements captured by the channel factors ξ i may significantly vary. One way to adjust for this is to consider weighted channels α i ξ i , with weights α i that can be trained [144,150] at run-time so as to minimise the resulting variance of the multi-channeled integrand I Γ .
More recently, the development of advanced machine learning frameworks such as TensorFlow and pyTorch facilitated the exploration of sampling methods based on various neural networks architectures [151][152][153][154]. In particular, the use of Normalizing Flows models, originally introduced in ref. [155] in the context of ray tracing and later refined in refs. [156,157], offers a promising novel approach [158][159][160][161] for improving adaptive sampling for Monte-Carlo integration beyond the factorised paradigm. The investigation of the potential of these recent developments in the context of Local Unitarity can ideally be carried out in conjunction with the vectorisation of our implementation so as to render it suitable for accelerated hardware, such as graphics and tensor processing units (see e.g. ref. [162]).

Numerical results
In these sections we showcase the performance of our implementation of the Local Unitarity representation in two different setups: • The computation of the NLO correction to the differential cross section of the scattering process e + e − → γ → dd.
• The computation of subsets of IR finite interference diagrams (i.e. supergraphs) contributing to the N 4 LO accurate cross-section of a 1 → 2 + X scalar scattering process.
All results reported in this section are obtained from our own implementation in a program called αLoop for the computation of differential cross-sections within the Local Unitarity framework.

NLO correction to e + e − → dd
In this section we present results for the NLO correction to e + e − → dd. In sect. 2.1 we have studied the two supergraphs of this process, namely the double-triangle supergraph and self-energy supergraph, and have shown local IR cancellations for the double-triangle supergraph in sect. 3.1. Contrary to those sections, we label here the loop momenta of these supergraphs ( k, l) and not ( k , l ). In sect. 6.1.1 we study the numerator and the renormalisation of the double-triangle, and visualise the integrand in sect. 6.1.2. In sect. 6.1.3 the self-energy treatment is applied to the self-energy supergraph. Finally, we present results for the differential cross-section in sect. 6.1.4.

The double-triangle supergraph
We write out the Feynman rules for the double-triangle supergraph and obtain the following numerator (common to all its Cutkosky cut contributions) in d = 4: where we included the averaging factor over incoming lepton helicity configurations as well as the flux factor. The additional factors of i and −i arising from the complex-conjugation of propagator -58 -numerators and vertices on the right of the Cutkosky cut always come in pairs so that we can ignore them for now 8 . The expression of the numerator of supergraphs such the one above will need to be computed repeatedly during the Monte-Carlo integration and are optimised using the procedure described in sect. 5.3.2.
Given that the LU representation realises the local cancellations of all infrared divergences of the integrand, we only need to consider counterterms for the regularisation of its UV behaviour. As is laid out in sect. 5.2, this is achieved by implementing a local UV subtraction procedure. The local UV counterterms are built from the original (unexpanded) numerator in (4-2 )-dimensions and from denominators expanded around the UV point. This expansion is made systematic by scaling all external momenta involved by a parameter λ and expanding around λ = 0. This procedure is well suited for automation, but to demonstrate the interplay between renormalisation and integrated UV counterterms, we factorise the triangle loop correction and apply the UV regularisation procedure to it only. We therefore rewrite the integrand stemming from the numerator of eq. (6.1) for the Cutkosky cut c s v 1 as follows (the case of c s v 2 is fully analogous): where K (s v 1 ,left) and K (s v 1 ,right) are the factorised pieces that are irrelevant for this discussion, while the one-loop corrected γdd-vertex Γ µ ij reads: Given that the superficial degree of UV divergence of the one-loop vertex correction is zero, the local UV counterterm denoted by I µ can be constructed by setting d = 4 in eq. (6.4) and expanding the denominators around their UV limit to first order. We must also assign a mass M 2 UV to the resulting denominators in order to prevent introducing new IR divergences in the UV counterterm: The denominator structure (k 2 − M 2 UV ) is now guaranteed to never yield a pole in its (c)LTD representation since the absences of any shift in the denominator implies that the existence conditions of the UV counterterm E-surfaces can never be fulfilled. This is different from similar UV local counterterms in four-dimensional Minkowski space where the squared mass term in (k 2 − M 2 UV ) needs to be set negative or imaginary (see, e.g. [166]) in order to avoid complications.
The integrated counterpart of eq. (6.5), denoted by ∆I s v 1 is straightforward and is derived in appendix. A.1 for completeness. We obtain: When considering a renormalisation scheme introducing complex-valued renormalised masses and couplings, such as the complex-mass scheme [163][164][165], the situation may be more involved. We however leave this aspect to a future investigation.
-59 - Traditionally, the wavefunction counterterms are provided in the on-shell scheme (see, e.g. eq. (6.8) of ref. [131]) and the renormalisation of the strong coupling in the mixed scheme where massless quarks are renormalised in MS and massive ones at zero momentum (see, e.g. ref. [167] or eq. B.5 of [13]). The renormalisation constant δα QED of the QED charge from QCD correction is of course zero. In order to convert the massless quark renormalisation constant from the on-shell scheme to the MS one, it is important to separately keep track of the dimensional regularisation pole of UV and IR origin in order obtain: We then find that the combination of the integrated local UV counterterm ∆I s v 1 in eq. (6.6) and of the renormalisation contribution δI (γdd) (left) is finite, as expected: (6.11)

Double-triangle supergraph integrand visualisation
In order to visualise the important features of the double-triangle supergraph integrands, we must choose a different projection than that of sect. 3.2.2 so as to allow the approach of the soft and collinear configurations of the gluon. We choose here again Q µ = p µ 1 + p µ 2 = (1, 0, 0, 1) + (1, 0, 0, −1) with the momentum routing of fig. 3 and set ( k, l) = ((0, k y , 1 √ 2 ), (0, 1 √ 2 , l z )) . (6.12) Moreover, we consider the integrand for the semi-inclusive cross-section defined by the following phase-space cut: 0.4 < p t,j1 < 0.8 (6.13) applied to the p t -leading jet reconstructed according to the anti-kt algorithm [168,169] with cone radius parameter ∆R set to 0.4. It is clear that this observable is not phenomenologically relevant for lepton collider experiments but it serves as case-study of a non-trivial observables. We start by showing in figs. 13 various combinations of the integrands of the interference diagrams produced by αLoop when turning off the complex contour deformation. Except when showing results for multi-channelling, the Jacobians from the loop momentum parameterisations is not included.
The first fig. 13a shows the contribution stemming solely from the Cutkosky cut c s v 1 . A first observation is the linear bands that correspond to the selection cuts of eq. (6.13). The observable -60 -function for this cut reads O( l, l − Q) so that the jet p t is dependent on l z only for our projection and the selection bands appear vertical. It may appear counter-intuitive at first that the transverse momentum depends on the z-component of the l loop momentum. This is however expected because the space presented in the figures corresponds to the dependency of the Local Unitarity integrand before the change of variables of eq. (3.7) (i.e. the sampling space probed by the Monte-Carlo integration). This change of variables results in a l z -dependent rescaling of the fixed component l y = 1 √ 2 which in turn induces an indirect dependence of p t,j1 on l z . The situation is opposite for the contribution of the cut c s v 2 , shown in fig. 13b, since in that case the observable reads O( k, k − Q) so that the jet p t depends only k y , thus yielding an horizontal band. The boundaries of the phase-space cuts are more complicated for the real-emission contribution, as expected from its complicated three-body observable dependence O( l, l − k, k − Q). It is interesting to note however that the central region is only populated by the real-emission kinematics, making it effectively LO accurate. This feature bears resemblance with that of other more common differential quantities in lepton collisions, such as for example the C jet-shape parameter [170,171] which only receives contributions at C = 0 from e + e − → jj and for C < 3 4 from e + e − → jjj. Another important feature from the set of figs. 13 is the behaviour of the integrands around the non-pinched E-surface located at k y = l z , since Soper's rescaling will simultaneously bring both norms | k| and | l| to Q 0 /2 = 1, thus placing the configuration on both the E-surface of the triangle loop and the Cutkosky cut c s v 1 or c s v 2 (see fig. 7). In the absence of a deformation, the singularity along the E-surface is left unregulated and this is reflected by the bright diagonal lines in figs. 13a and 13b. The local cancellation of these non-pinched E-surfaces is then made apparent in fig. 13c where the areas accessible to both Cutkosky cuts c s v 1 and c s v 2 is devoid of any particular features along the diagonal (point labelled "B"). However, the E-surface cancellation breaks down on points labelled "A" and "D" where the observable function (i.e. phase-space cut in our present case) retains only one of the two cancelling terms. This shows that fully 9 inclusive cross-sections can be computed without a contour deformation.
We now consider the soft and collinear gluon kinematic configuration. Due to Soper's rescaling, the collinear subspace does not only lie at the soft end-point k y = l z but instead lies along an arc denoted by "Coll" in the upper-left portion of figs. 13. Each of the two virtual Cuktkosky cuts covers only "half" of the collinear space since our particular choice of signs for the fixed spatial components of k and l only allows one of the two pinched E-surfaces of the loop triangle to be reached (either the qg or theqg collinear one). Together however, the virtual contributions of fig. 13c reproduce the complete collinear singularity to be locally cancelled by the sum of real-emission contributions shown in fig. 13d. The overall sum in fig. 13e presents no particular feature close to collinear regions any longer.
The soft singularity in fig. 13e appears dampened but is still unbounded, as expected from the power counting arguments of sect. 3.2.7. This integrable singularity is particularly harmful for a numerical implementation since for the choice of momentum routing of fig. 2, the soft singularity lies at k − l = 0 which spans a full three-dimensional volume of the phase-space parameterised with ( k, l). In our particular example, this can be remedied simply by adopting a different loop momentum basis containing the gluon propagator as a defining edge. This solution is not generic, since for a triple gluon vertex it is not possible to consider all three connected gluons to have an 9 "Fully" in this context implies that all possible Cutkosky cuts of each contributing supergraph is considered in the entirety of the phase-space available to it. By definition, this therefore excludes any 1 → N, N > 2 process in the Standard Model -61 -independent momentum. The general solution then involves combining multiple parameterisations using a multi-channeling technique discussed in sect. 5.4.1 and analogous to the one we already presented in sect. 5.2 of ref. [89]. Fig. 13f showcases the suppression of the soft-singularity from multi-channeling. Notice that in all figures not involving multi-channeling, the integrand weight considered did not include the parameterisation Jacobian, since the objective was to focus on the behaviour of the integrand only. For our default choice of momentum routing, the Jacobian is anyway relatively shallow within our chosen sampling space projection. For multi-channeling, the combination of the three non-degenerate loop momentum bases of the supergraph results in a parameterisation Jacobian which is precisely the quantity responsible for achieving the dampening of the integrable singularities and it is therefore crucial that it be included in the weights displayed by the multi-channeling density plot. Furthermore, the multi-channelling figures are generated using the exact αLoop integrand function exposed to the integrator, which implies that it includes the numerical stability estimate and rescue system described in sect. 5.3.3 which exclaims the sparse spots in the density plot where the integrand is sent to zero as it failed the numerical stability requirement threshold. We find the expected suppression of the soft integrable singularity in fig. 13f, but also see a complicated mangling resulting from the combined contributions of our input k and l re-interpreted for all three relevant loop momentum bases of the double-triangle supergraph. We stress that the visual complexity of the integrand can be deceptive since in practice highly discontinuous bounded functions with a lot of structures often converge significantly faster than smoother integrands that have integrable singularities, especially when considering adaptive Monte-Carlo sampling. We note however that, as shown in fig. 16c, the multi-channeling strategy renders the integrand bounded only in the absence of a deformation. When a deformation is enabled, there remains an integrable singularity along the non-pinched E-surface cancelation diagonal line, which is however much less severe. The local UV counterterm discussed in sect. 6.1.1 is active as part of the integrands shown, but has no visible impact for our chosen value of M 2 UV = 2Q 2 and kinematic range (k y , l z ) considered.
In figs. 14 we visualise the double-triangle integrands with the dynamic deformation of sect. 5.3.4 enabled. The integrand visualised corresponds to the integrand I Γ ( k) introduced in eq. (5.11). When multi-channeling is considered, the resulting integrand I  fig. 2d. We remind the reader that the contour deformation is not the one constructed for the entire double-triangle supergraph (see vector field of fig. 9), but it is constructed independently for each loop integral remaining in each Cutkosky cut contribution, so as to have real-valued momenta in the observable function and allow to separately accommodate the complex-conjugated causal prescription applying to loops appearing on the right-hand side of the Cutkosky cut. This implies that in the case of the cuts s v 1 and s v 2 the remaining loops are triangles and their contour deformation must be constructed dynamically for each three-point external kinematics induced by the particular sampling point considered. The general solution for this problem discussed in sect. 5.3.4 is greatly simplified in the case of the doubletriangle supergraph as one can show that the maximal overlap structure (see ref. [89]) contains only a single E-surface (the only one present) which moreover always admits the origin in its interior. Thus a radial deformation field with the origin as a source is a valid deformation. The integrands of figs. 14 share most of the features already seen in figs. 13 with the main difference being the expected regularisation of the E-surfaces, for example at the points labelled "A" and "D". When following the collinear arc, we find blue regions (i.e. values approaching zero) in the density plots of Im[I s v 1 ] (fig. 14d) and . 14e) resulting from our dampening of the deformation -62 -on infrared singularities so as to retain the cancellation with their real-emission counterpart. This results in a non-trivial structure of the complex phase of the complete integrand ( fig. 14g) close to the integrable soft singularity. The multi-channeling procedure successfully dampens that soft singularity and yields an integrand whose real part ( fig. 14i) is now bounded and suitable for numerical integration. Note that the many discontinuities of the integrand are now not only related to the phase-space cuts but also stem from the Jacobian of our contour deformation, which is discontinuous due to the presence of a min function in the functional form of the dynamic normalisation of our deformation (see eq. 3.33 of ref. [89]). Despite the introduction of a contour deformation, the differential cross-section is guaranteed to be real-valued as it eventually corresponds to the norm of an overall complex-valued amplitude. In the particular case of LU, this reality condition is realised through the existence of a complexconjugated partner of each loop present in each Cutkosky cut contribution of each supergraph. The cancellation of the imaginary component of the cross-section only holds at the integrated level since in general it even involves terms belonging to different supergraphs (contrary to the accidentally symmetric double-triangle case where the two complex-conjugated triangle-loop partners happen to belong to the same supergraph). We use this fact to our advantage as a strong cross-check of the correctness of the LU implementation in αLoop as well as an independent assessment of the reliability of the numerical accuracy reported by the Monte-Carlo integrator.
The series of figs. 14 capture the core features of the LU integrand that we want to highlight, but it does not reveal more detailed aspects such as the complex phase of the argument since we always considered absolute value of its real and imaginary component. We therefore complement our visualisation with figs. 15, where the colouring reveals the phase of the integrand and its shading marks isolines reveals the logarithmic progression of its magnitude across the plane. We find the expected phase-shift of π/2 of the virtual integrands when crossing the soft singularity along the collinear line as the deformation infrared dampening kicks in. When comparing the phase of I s v 2 in fig. 15a and I s v 1 in fig. 13b on their respective side of the collinear singularity where their deformation is active, we find a phase-shift of π stemming from the complex-conjugation of the causal prescription of the triangle loop placed on the left-hand side of the Cutkosky cut, which flips the sign of the deformation vector κ. The combination of the virtual contribution shown in fig. 15c has a constant phase along the collinear singularity since the real part is divergent and negative. Once the positive divergent real-emission contribution of fig. 15d is added, we find in fig. 15e that the cancellation of this large negative real part leaves off a remnant with a complicated phase structure. We observe in fig. 15f that the multi-channeling procedure dampens the integrable singularity without disrupting much of the phase-structure of the overall integrand I Σ around that region, as one expects since the multi-channeling Jacobian responsible for the dampening is purely real.
Finally, the interesting detailed structure around the soft configuration is hard to resolve within the selected range. For this reason, we also present alternative visualisations homing in on the region of interest in figs 16. This reveals many discontinuity lines that result from the interplay of the various competing constraints on the deformation normalisation realised through the different functional forms of the scaling factors λ i of eq. 3.33 of ref. [89] which eventually send the deformation to zero on the soft point. The clipped regions of the 3D visualisation stem from discontinuities of I Σ while the few "spots" on the multi-channeling plots come from evaluations deemed unstable by our strict precision requirements and consequently sent to zero. Comparing the normalisation of the z-axis scale between figs. 16b and 16c highlights the impact of the multi-channeling on the integrable soft singularity. We note however that when a deformation is enabled (and only then), the multi-channeling treatment still leaves the integrand unbounded when approaching the soft singularity from the diagonal E-surface cancelling direction. This remaining integrable singularity does not severely increase the variance of the integrand for the double-triangle supergraph.  ].
-65 -  Figure 14. Similar visualisations as in figs. 13 but for the absolute value of the real and imaginary parts of the Local Unitarity integrands of the double-triangle supergraph when enabling the contour deformation discussed in sect. 5.3.4. The contribution from the real-emission 2I s r 1 is not shown since it is not subject to a contour deformation and it is therefore identical to that of fig. 13d. For the same reason, we have that fig. 14f. -68 -

The self-energy supergraph
We now consider the self-energy supergraph. The main novelty for the self-energy supergraph compared to the double-triangle supergraph is the presence of a repeated propagator, which requires special treatment since both of them are put on-shell by the Cutkosky cut. We describe our general treatment of the self-energy renormalisation in sect. 4.1 which addresses this problem by considering a combination of the on-shell scheme renormalisation conditions together with a Taylor expansion of the self-energy around its on-shell point. We observe that the self-energy insertion does not require a complex contour deformation, since the pairwise cancellation of all E-surfaces of the self-energy corrections present on either side of the Cutkosky cut in c s v 1 (fig. 2a) and c s v 2 ( fig. 2b) can never be spoiled by the observable function since it has the same kinematic dependence for both cuts. That is: 14) The local UV counterterm is easily obtained from eq. (6.18), since its superficial degree of divergence of zero implies that only the leading term of the UV expansion needs to be retained while keeping the original numerator unexpanded: It is important to build the UV counterterm only after the on-shell Taylor expansion has been constructed, since it does not commute with the UV expansion.
In appendix A.2 we provide the details of the computation of the integrated counterpart of the local UV counterterm of eq. (6.19). In order to write the result with the LO contribution factorised, we sum the two integrated UV counterterms contributions arising from both self-energy loops to the left (c s1 ) and right (c s2 ) of the Cutkosky cuts: We then add the renormalisation contribution δI We stress that we must consider the wavefunction renormalisation in the MS scheme and not in the on-shell scheme. Our procedure still realises renormalisation in the on-shell scheme, in virtue of our particular choice for the self-energy on-shell expansion given in eq. (4.7). When combining the integrated UV counterterms and renormalisation of eq. (6.11) for the triangle loop appearing in each of the two virtual Cutkosky cuts of the double-triangle supergraph as well as eq. (6.21) for both isomorphic self-energy supergraphs, we find the overall contribution from the UV regularisation to be − 1 2 C F αs π I (LO) . The cancellation of the logarithmic contributions log M UV results directly from the gauge cancellation between the vertex and quark wavefunction corrections. It also implies that the corresponding NLO QCD correction to e + e − → dd has no residual dependency in the renormalisation scale as expected, since its Born contribution does not depend on α s (µ 2 r ).

Results for the differential cross-section
For all numerical results and visualisations presented in this section, we choose Q µ = p µ 1 + p µ 2 = (1, 0, 0, 1) + (1, 0, 0, −1), Q d = − 1 3 , Q e = −1, g EW = 0.307954 and g s = 1.21772. The choice of the renormalisation scale is irrelevant for the process e − e + → γ → dd studied in this section, but matters for the contribution of the individual self-energy and double-triangle supergraph which we report here as well. We therefore specify that we set µ r equal to our choice of M UV = Q 0 2 = 1, such that the integrated UV counterterms in eq. (6.11) and eq. (6.21) feature no logarithmic terms.
We present our results in table 2 for the inclusive and semi-inclusive cross-section within the acceptance cuts of eq. (6.13). When a deformation is enabled, we keep the defaults for the deformation of ref. [89] with an overall maximal deformation magnitude scaling λ max = 1, together with its dampening on infrared singularities described in sect. 5.3.4. Table 2. Numerical results for the NLOQCD accurate (semi-)inclusive cross-section for the process e + e − → γ → dd. The second column entitled "C.D." stands for "Contour Deformation".
The timing indicated by t /p min corresponds to the single-core time spent in the C routine for the double-precision evaluation of each numerator involved per sampling point. In this sense, this is the minimal time that one evaluation could possibly take. The timing indicated by avg corresponds to the effective actual timing per sampling point, including the overhead from the integrator, observable functions, the deformation construction, and most importantly additional processing in αLoop such as the stability tests that necessitate additional evaluations in double and sometimes quadruple precision arithmetic (see sect 5.3.3). At LO, the timing is dominated by the overhead in αLoop, which explains the significant increase when enabling the phase-space cut selection code. Both timings reported are for one sampling point including the multi-channeling treatment, which in the case of the process e + e − → γ → ddg implies three independent evaluation of the integrands for each of its three channels integrated together. The imaginary part of the integrated cross-section is not reported but was found to always be compatible with zero given the Monte-Carlo uncertainty. The comparison of the result with and without a contour deformation enabled (only when the latter is possible) shows a moderate impact of a factor of about three in the variance of the integrand. The increase in the runtime per sample when enabling the deformation is mostly due to the increase in fraction of unstable point that require a quadruple precision rescue. Note that since the heuristics described in sect. 5.3.4 for the computation of the deformation centres always apply for this simple process, there is no need for numerically solving a second-order cone program at run-time and therefore not a lot of time is spent in constructing the deformation. As can be seen in tab. 2, the infamous NLO K-factor of α s /π (for entry 2×SE+DT) is very well reproduced by our implementation of numerical Local Unitarity.
We complement the (semi-)inclusive results of tab. 2 with the binned distribution of the transverse momentum of the leading reconstructed jet in fig. 17. This observable is of little interest in the case of lepton collisions but was chosen to demonstrate the applicability of Local Unitarity to arbitrary observables. To the best of our knowledge, this result of humble appearance stands as the first NLOQCD accurate differential binned distribution ever computed without relying on IR  Figure 17. Comparison of the differential prediction for NLOQCD correction to the observable pt(j1) for the process e + e − → γ → dd as obtained from MadGraph5_aMC@NLO and αLoop . The Monte-Carlo sampling for this run with αLoop is 5 · 10 9 points. The dashed pattern indicates when a contribution is negative.
counterterm and/or dimensional regulator.
Overall, we find excellent agreement between αLoop and MadGraph5_aMC@NLO in both tab. 2 and fig. 17. In particular, discrepancies are always within a couple of Monte Carlo standard deviations, with a relative accuracy of the NLO correction always below the percent mark, except for the region 0.2 < p t (j 1 ) < 0.4 where the distribution changes sign and therefore neighbours zero. These results and especially the reliability of the Monte Carlo accuracy reported, are indicative of the predictive power of our numerical implementation of Local Unitarity.

Individual scalar supergraphs
An important verification of the local properties of the LU representation is by testing whether the IR cancellation pattern is realised in practice in a Monte Carlo integration.The most direct test consists of integrating the overall contributions obtained by summing together all the interference diagrams corresponding to a fixed, higher-loop scalar supergraph Γ in the LU representation, σ Γ = dσ Γ /dO with constant observable. Scalar graphs do not have a numerator that could improve the IR behaviour and are not helped by gauge invariance or any property specific to physical theories in order to realise IR cancellations.
In table 3, we enumerate the massless scalar supergraphs whose LU representation can be numerically integrated using a Monte Carlo procedure. We divide these scalar supergraphs in three series: the diagrams labelled by a correspond to all of the three-loop supergraphs whose σ Γ is nonzero and has no self-energy insertions; the diagrams labelled by b correspond to all of the four-loop supergraphs whose σ Γ is finite, non-zero and has no self-energy insertions; and we selected three five-loop supergraphs, labelled by c.
The integration results, obtained via our implementation named αLoop, are exhibited in table 4, along with the analytic result obtained through forcer [48]. We also report the total number of sample points used for the integration, N p , the minimum evaluation time per sample and per core (that is, the time for one single evaluation of the integrand in double precision on one Intel(R) Xeon(R) Gold 6136 CPU @ 3.00GH core) and the average evaluation time per sample, which includes the sum over all integration channels, the numerical stability checks and the dynamic stability rescue mechanisms. We also report the number of integration channels N ch and the ratio ∆ [σ] of the discrepancy w.r.t the exact result with the Monte Carlo error, and the relative size ∆ [%] of this discrepancy.
We find that all numerical integrations yield results that are within three sigma of the analytic result with Monte Carlo errors at the percent level or better. Overall, the results undoubtedly confirm that the integrands are stable and well-behaved. We observe that the b.16 supergraph requires double the sampling points to achieve a relative error below 2%.
The differences in the evaluation time across the topologies considered, both minimal and average, is related to a variety of factors, including the number of interference diagrams arising from the supergraph and the number of integration channels. Another determining factor is the overall numerical stability of the integrand: indeed, a worse stability implies that a larger fraction of sampling points trigger the numerical stability rescue mechanism, which significantly slows does the evaluation of the integrand which is then performed using quadruple precision arithmetics.  Table 3. List of scalar supergraphs whose corresponding cross-section contribution is computed in 4. They are divided in three series: a) 3-loop supergraphs, b) 4-loop supergraphs, c) 5-loop supergraphs. Grey edges and black edges do not intersect each other, but overlap due to the planar embedding.  Table 4. Values of σΓ (that is, the integrated sum of all interference diagrams arising from Γ, or inclusive cross-section per supergraph) for the supergraphs listed in 3 computed with αLoop for Q µ = (1, 0) and compared with benchmark results from forcer [48]. The table includes the total number of samples used for each integration, Np, the number of integration channels Nch (see sect. 5.4.1), the minimum and average per sample evaluation time of the integrand t/p, and the relative discrepancy with respect to the standard deviation ∆ [σ] and numerical result ∆ [%].
-74 - In this paper, we have presented a general and systematic way to construct a representation of differential cross-sections that is locally free of IR divergences. Given an observable and scattering process, the procedure yields an integrand free of (non-integrable) singularities and whose integral reproduces physical differential cross-sections. We call the new expression of the differential crosssection its Local Unitarity (LU) representation. Our framework ultimately achieves two objectives: first, it shows that the appearance of infrared singularities is a feature of traditional calculation methods and not a fundamental attribute of physical theories (a statement already implicitly contained in the KLN theorem) and second, it provides a concrete method well-suited for an automated and fully numerical computation of physical observables. Achieving both of these objectives inherently requires an in-depth study of the singular structure of Feynman integrals, of phase space integrals and their relationship. We use the Loop-Tree Duality (LTD framework) to express loop integration on the same footing as phase space integration. Once amplitudes are re-expressed using their LTD representation, a clear picture emerges: all the interference diagrams that correspond to a particular Cutkosky cut of the same forward scattering topology (i.e. supergraph) can be collected under the same integral sign, and their sum is finite. This principle effectively establishes how to collect and group interference diagrams into classes that are infrared finite. We even identify that cancellations between IR singularities happens pairwise. More specifically, if an interference diagram has an IR singularity, there exists a unique distinct other interference diagram which exhibits the exact same IR singularity so that they sum up to an integrable quantity on the location of said singularity.
One of the core issues in constructing an integrand that realises IR cancellation locally while simultaneously yielding differential cross-sections once integrated, is that of aligning the phase space measures. Part of this task is performed through the LTD formalism, along with the choice of a consistent loop momentum across all interference diagrams corresponding to the same supergraph. The other fundamental component allows one to solve simultaneously all of the Dirac delta functions imposing the energy conservation of on-shell external particles: the causal flow. The simplest form of this flow amounts to a trivial rescaling of all loop variables, but its most generic solution naturally follows from the general solution to a contour deformation encoding the Feynman prescription associated to each propagator. The causal flow is an essential ingredient to the construction of the function whose integral yields the physical cross-section, and is used for the all-order proof of local cancellations of final-state IR divergences in the LU representation. A corollary of our investigation reveals how gauge invariance, and more specifically the IR behaviour of diagrams featuring ghosts, together with the LSZ formalism accommodate IR cancellations in external self-energy corrections.
Aligning the integration measures also provides a solution to the problem of automating the calculation of differential cross-sections, as it allows for the systematic construction of an integrand especially well-suited for Monte Carlo integration. In order to corroborate and demonstrate the workings of the LU representation in practical cases, we provided the results of the numerical integration of the LU representation for a variety of example scalar supergraphs. Furthermore, we explicitly unfolded our theoretical construction for the NLO correction to the differential crosssection of the physical process e + e − → dd.
This paper opens up new horizons both in the understanding of the analytical properties of observable quantities calculated perturbatively in QFTs and in the ambitious goal of achieving the complete automation of the calculation of physical cross-sections for collider experiments beyond -75 -NLO accuracy. It goes without saying that, ultimately, these two goals are inextricable, and can potentially lead to a radically different understanding of QFTs as we know them.
We expect the foundational work laid out in this paper to serve as a stepping stone for its many implications, both regarding the theoretical understanding of initial and final state infrared singularities and its concrete applications to the computation of higher-order perturbative corrections for collider phenomenology.
A Integrated UV counterterms for the NLO correction to e + e − → dd

A.1 Double-triangle supergraph
In this section, we explicitly derive the result given in (6.6). We first drop the term linear in k that integrates to zero and arrive at: The tadpole tensor integral in (A.1) can be evaluated using the usual decomposition in the only manifestly Lorentz invariant structure g ρ1ρ2 : followed by an elementary reduction step: together with the general expression for a scalar tadpole integral (see for instance eq. (3) of ref. [172]): where we introduced the shorthand notation log M UV for the omnipresent logarithm involving the UV mass regulator M 2 UV since the tadpole nature of all integrated counterterms is such that it will always appear together with µ 2 r and never any other scale. The presence of the negative sign in −M 2 UV calls for a regulator in order to set the sign of the imaginary part of the logarithm. The local UV counterterm of (6.5) will be cast into its (c)LTD representation while assuming a causal Feynman prescription +iδ in the denominator of the UV propagator. This choice dictates what is the correct branch to select for the evaluation of log M UV . For this process at NLO we note that the particular sign of this imaginary part is irrelevant as it -77 -would always cancel between the two integrated counterterms arising from the loop integrals sitting to the left and right of the Cutkosky cut, since in the latter case the complex conjugation flips the sign of the causal prescription.

A.2 Self-energy supergraph
We provide here a detailed computation of the integrated counterpart of the local UV counterterms given in (6.19) for the self-energy supergraph appearing in the NLO QCD correction of the process e + e − → dd. To this end, it is useful to write it in a manifesly Lorentz invariant fashion by introducing the reference vector η µ = (1, 0, 0, 0, . . . ), where the dots denote the (d − 4)-dimensional orthogonal subspace. We find where we have immediately removed the linear contribution from the numerator since it integrates to zero and recycled the expression for the tensor integral found in A.2. In order to obtain the complete integrated UV counterterm, we must re-insert ∆I (se) into (6.18) and (6.17): where we used the identity γ µ γ ν γ µ = (2 − d)γ ν . In order to demonstrate the factorisation of the LO supergraph, we must consider the sum of the contributions from both Cutkosky cuts c s v 1 and c s v 2 . The cut c s v 1 yields a derivation fully analogous to the one carried out in this section, with the exception that the on-shell propagator being regulated is sitting to the right of self-energy, so that one naturally arrives at the same expression as (A.10) but with the Lorentz structure / η/ k in place of / k / η. We then observe that: where we used the fact that both Cutkosky cuts of the SE supergraph force the on-shell condition k 0 = k . Using the factt hat that K