A numerical evaluation of planar two-loop helicity amplitudes for a W-boson plus four partons

We present the first numerical results for the two-loop helicity amplitudes for the scattering of four partons and a W-boson in QCD. We use a finite field sampling method to reduce directly from Feynman diagrams to the coefficients of a set of master integrals after applying integration-by-parts identities. Since the basis of master integrals is not yet fully known analytically, we identify a set of master integrals with a simple divergence structure using local numerator insertions. This allows for accurate numerical evaluation of the amplitude using sector decomposition methods.


Introduction
The growing precision of high energy collider experiments puts increasing strain on our ability to make reliable theoretical predictions. Standard techniques for the computation of perturbative scattering amplitudes often fail when applied to the multi-loop and multi-leg processes currently produced in abundance at the LHC. Constantly evolving methods have led to differential predictions at next-to-next-to-leading (NNLO) for 2 → 2 scattering and N 3 LO for 2 → 1 processes 1 .
The need to match experimental precision has led to increasing efforts from the theoretical community to develop new techniques for 2 → 3 predictions at NNLO. The first hurdle has been to compute unknown two-loop amplitudes in which the analytic and algebraic complexity causes conventional approaches to integral reduction to fail. Major new advances that exploit numerical evaluations over finite fields [3][4][5] have recently produced the first analytic results for five-parton amplitudes in the leading colour approximation. Combined with the recently computed analytic master integrals [6,7] using the canonical basis approach to differential equations, a form suitable for combination with the unresolved contributions to the cross section has been obtained [8][9][10][11].
The production of a W -boson together with jets at hadron colliders are important signatures that can be used as precision probes of the Standard Model. QCD corrections to W +jets have been a traditional testing ground for new technology. pp → W + j was among the first 2 → 2 process computed at NLO [12]. The amplitudes for pp → W + 2j were computed using the recently developed on-shell unitarity method [13][14][15][16] and were implemented into MCFM to provide differential cross-section predictions [17]. NLO results for associated W -boson production with three or more jets are accessible through automation and the use of the generalised unitarity method [18][19][20][21][22][23].
NNLO corrections to pp → W + 2j will open up possibilities for further precision tests of the Standard Model. The two-loop amplitudes are obvious targets for the new technology developed for massless five-point amplitudes, yet the off-shell vector boson adds an extra scale and therefore a new layer of complexity. The first step towards a complete analytic computation is to set up a procedure that could evaluate the amplitudes numerically using rational kinematics. It is this benchmark evaluation of the amplitudes that is the subject of this paper.
The computation of higher order corrections to perturbative scattering amplitudes is a well studied problem. Amplitudes with two or more loops have relied on the technology of integration-by-parts (IBP) [24,25] reduction, which in recent times has involved following Laporta's algorithm [26], together with numerical or analytic methods for the evaluation of the resulting basis of master integrals. For these multi-scale basis integrals with massless internal propagators the differential equation technique [27][28][29] has been employed to find analytic expressions, most recently for the complete set of planar [7,8] and non-planar integrals [30][31][32][33][34]. For the case, in which the amplitudes considered here fall, only one of the three planar families has been evaluated [6]. Combining the master integrals into complete amplitudes requires the solution of increasingly complicated linear systems of IBP equations. Considerable effort has led to a variety of efficient solutions [3,[35][36][37][38] and public implementations [39][40][41][42][43]. Applications to five-particle problems have been possible though yielded large IBP reduction tables [44,45]. In this paper, we only perform the IBP reduction numerically over finite fields in order obtain the coefficients of the amplitude in terms of master integrals. As shown e.g. in refs. [5,9], when combined with functional reconstruction techniques, this approach also allows to directly reconstruct analytic results for amplitudes, sidestepping the need of computing and using large analytic IBP tables, which are often significantly more complicated.
Another important ingredient has been the development of efficient methods to construct on-shell integrands and integral coefficients. Integrand reduction techniques [46] combined with the use of a Feynman diagram approach or generalised unitarity have been very successful for the computation of one-loop amplitudes, in particular to construct scalar integral coefficients numerically. These techniques have been extended to two loops [47][48][49][50][51][52][53] and methods to employ unitarity cuts [13,14] to build amplitudes by directly incorporating IBP decomposition have been established [37,[54][55][56][57].
The first steps towards helicity amplitudes for five-point amplitudes were taken through numerical evaluations of two-loop five point amplitudes in QCD using modular arithmetic [58][59][60][61]. These algorithms have been generalised to allow for a full reconstruction of the coefficients of the pentagon functions classified in [7] leading to an analytic form of the single-minus helicity amplitudes [9] and the complete leading colour five-parton helicity amplitudes within the numerical unitarity framework [10,11]. The success of computations in the planar sector has shifted focus to the non-planar sector of massless two-loop fivepoint amplitudes with a series of new results in super-symmetric Yang-Mills [32,62] and gravity [63,64] as well as in the all-plus sector of QCD [65].
In this paper we consider the case of planar amplitudes with an off-shell external leg. We apply the recently developed technology for the computation of two-loop five-particle amplitudes using sampling of Feynman diagrams over finite fields. Using a modular approach, recently presented as part of the FiniteFlow algorithms [5], we are able to numerically evaluate the diagrams and perform an integrand reduction, subsequently reducing the resulting integrals using integration-by-parts identities. Since the complete set of analytic master integrals is not known, some of the integrals were evaluated numerically using sector decomposition [66][67][68]. Analytic results for the following classes of master integral are available: one of the three families of the off-shell five-point pentagon-box [6] and fourpoint functions with one [69] and two off-shell [70][71][72][73] legs. For master integral topologies for which a numerical evaluation through sector decomposition is challenging, we identified a basis of master integrals using local numerators [74,75] with simplified divergence structure and therefore easier numerical evaluation. We consider both the qQQq ν and qggq ν sub-processes in our computation where the decay of the W -boson is also incorporated.
We describe our integrand reduction setup that is subsequently interfaced to IBP reduction in Section 2. In Section 3 we discuss the structure of the leading colour W +4 parton amplitude at two loops including its singularity structure. The identification of a master integral basis with local numerator insertions is elaborated in Section 4. Finally, we present numerical benchmark results for both sub-processes in Section 5 and draw our conclusions in Section 6.

Calculational framework
The framework described in this section is a modification of the numerical algorithm for two-loop amplitudes presented in [9] to allow for the use of integrands built from Feynman diagrams as an alternative to generalised unitarity cuts in six dimensions. While the unitarity method can be very efficient, a fully numerical approach, with rational reconstruction, is also able to avoid the traditional problems associated with the Feynman diagram approach.
We start by generating a set of Feynman diagrams using Qgraf [76] and performing colour decomposition to separate the colour parts of the amplitude from the kinematic parts that depend only on external momenta {p}. We obtain where A where k i is the loop momenta, d = 4 − 2 is the space-time dimension and d s = g µ µ is the spin dimension. The loop amplitude in t'Hooft-Veltman (HV) scheme [77] can be obtained by setting d s = d, while the Four-Dimensional-Helicity (FDH) scheme [78] can be achieved by setting d s = 4. Each numerator function, N T (d s , {k}, {p}), that contains numerators of Feynman diagrams that share the same diagram topology, is processed by applying the t'Hooft algebra. This is carried out with the help of Form [79,80] and the Spinney library [81]. In general, the explicit functional dependence of the numerator function at this point is given by is a spinor string made up of slashed momenta ( / q i and / k i ). The d-dimensional loop momenta can be decomposed into a four-dimensional part and an extra-dimensional part Due to rotational invariance in the extra dimensions,k i can only appear in the numerator function as µ ij = −k i ·k j . We obtain helicity amplitudes by specifying the helicity/polarisation of each external particles and we further parametrise the dependence on the external kinematics by using momentum twistor variables, x i [82]. This allows us to express the spinor products of external momenta ( ij , [ij]) and Mandelstam invariants (s ij ) uniformly in terms of momentum twistor variables, where momentum conservation and spinor product relations like Schouten identities are already built in. At this point we are considering the two-loop n-point helicity amplitude where the explicit functional dependence on the helicity-dependent numerator function In processing the algebraic expressions in Eqs. (2.1), (2.2) and (2.5), we have used in-house Form [79,80] and Mathematica scripts. In this form the helicity-dependent numerator functions in Eq. (2.5) must be reexpressed in terms of integral families that can later be reduced using IBP equations. To achieve this we apply an integrand reduction algorithm to obtain where ∆ is the irreducible numerator for an independent topology T . In order to determine ∆, we first need to construct a basis of irreducible scalar products (ISPs). We opt to use a basis of ISPs in terms of auxiliary propagators that is suitable for IBP reduction.
To build an IBP compatible integrand basis, we define an integral family G a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9 a 10 a 11 = where p ij···k = p i + p j + · · · + p k . Up to cyclic permutations of the external legs, all integrals appearing in Eq. (2.7) can be written in the form of Eq. (2.8). We follow the conventions used in [9] where negative exponents, a i < 0, correspond to the ISPs of the irreducible numerator, ∆ h T . The irreducible numerators are the most general polynomials in the ISPs with exponents bounded by renormalisability conditions. As an example, the parametrisation for the two-mass double-box topology is where the figure represents the topology T . The bounds on the exponents are −6 ≤ a 2 + a 9 + a 10 + a 10 ≤ 0. (2.12) The helicity-dependent coefficients are functions of the spin dimension, d s , and the external kinematics, To determine the coefficients we express the numerators in terms of the propagators and ISPs. This is achieved by expanding the loop momenta in terms of external momentak The coefficients of the spanning vectors, p j , are functions of the inverse propagators and ISPs, a ij = a ij (D α , ISPs). a ij can be determined by solving a linear system of equations constructed by contracting Eq. (2.13) with the spanning vectors, p j . All variables in the numerators Eq. (2.6) can then be expressed in terms of these coefficients. For example 14) The last relation is obtained by squaring Eq. (2.4). The variables are straightforwardly evaluated on generalised unitarity cuts by setting all propagators to zero without relying on explicit loop momenta solutions to the cut constraints. We observe that this form is the starting point for the derivation of the Baikov representation [83], which is obtained by integrating out angular dependence in the space transverse to the external momenta. We note that the change of variables to rewrite the numerators in terms of propagators could be performed directly. However, the choice to apply the substitution using the integrand reduction approach breaks the problem into a series of linear systems with fewer parameters, rather than one large system. At this point we can solve for the coefficients of the integrand parametrisations by equating them to the diagram numerators. Using the two-mass double-box as example again, we have the cut equation which is valid only when both sides are evaluated on the hepta-cut for the two-mass double box, i.e. D α = 0 for α = 1, 3, . . . , 8.
After setting up the integrand reduction system, we write the helicity amplitude as a linear combination of integrals in the integral family of Eq. (2.8), where we sum over tuples a = (a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 , a 8 , a 9 , a 10 , a 11 ). At this stage IBP reduction can be applied to the integrals G a . This basis choice allows for a simple interface to the IBP reduction. The integrand is never reconstructed analytically but only sampled numerically. The final results are independent of the particular parametrisation and therefore the exact form of the integrand is not the main concern here. Nevertheless, alternative representations of the integrand, e.g. [84], may lead to improved efficiency. The integrals with non-zero coefficients after numerical sampling of the integrand are reduced to a set of master integrals via IBP identities. The IBP relations are generated in Mathematica using the Laporta approach [26] with the aid of LiteRed [39], and solved numerically over finite fields within the FiniteFlow framework [5]. We can finally write the helicity amplitudes in the master integral basis J k (2.18)

Planar two-loop W plus four parton scattering
The number of Feynman diagrams contributing to qQQq ν and qggq ν processes at leading colour are 210 and 603, respectively, and the leading colour partial amplitudes are extracted according to, where n = m N c α s /(4π), α s = g 2 s /(4π) and m = i(4π) e − γ E . g s and g W are the strong and weak coupling constants respectively. We note that the vector boson only couples to the quark line connecting q andq and does not couple to the equal flavour quark pair Q,Q.
We choose a rational parametrisation of the massless 2 → 4 kinematics using the momentum twistor parametrisation [82] where y i = i j=1 j k=1 1 x k and We stress that while we generate a complete parametrisation for the 2 → 4 scattering process, analytic expressions could be obtained with only six independent parameters since the decay of the W boson completely factorises. Since it is easy to generate a rational parametrisation for n-particle scattering of massless particles, it is simplest to start from a configuration including the decay of the W boson. The leading colour partial amplitude is passed through an integrand reduction stage which projects onto a basis of 453 topologies with irreducible numerators written into the basis of the 15 maximal cuts shown in Figure 1. The remaining integrals are then passed through a Laporta style IBP reduction to find a basis of 202 master integrals (including the 5 cyclic permutations). The distinct master integral topologies are shown in Figures 2 and 3. The most complicated integrals that need to be reduced are rank 5 pentagon-boxes, e.g. G 11111111−3−1−1 according to the notation defined in Eq. (2.8).
Once the amplitude is decomposed in terms of master integrals and the evaluations of master integrals are available (either analytically or numerically), we can perform a Laurent expansion in the dimensional regularisation parameter, . The -expanded partial amplitude contains a divergent part, manifested by the poles in , and a finite part. The infra-red (IR) divergent part of the partial amplitude, obtained after removing the ultra-violet (UV) divergences by introducing a set of counter-terms, is universally known [85][86][87][88]. Figure 1: Independent maximal cut topologies contributing to planar W + 4 parton scattering at two-loops. The full set of 15 maximal cuts can be obtained by including 2 permutations of A 1 , A 3 , B 1 , B 2 , C 1 and C 2 topologies.
The pole structure of the unrenormalised amplitude in the HV scheme at one and two loops is given by whereÂ (1) is the unrenormalised one-loop amplitude normalised to the tree-level amplitude. The I 2 ( ) operator is defined by while the I 1 ( ) operators for the W + 4 parton process at leading colour are where N ( ) = e γ E /Γ(1 − ) and Note that the H qQQq ν ( ) and H (2) qggq ν ( ) functions are given in the leading colour limit. The β function coefficients and anomalous dimensions without the contribution from closed fermion loops N f are 14) . With the pole structures available, we can check that the divergent part of the twoloop amplitude agrees with the predicted UV and IR poles. On the other hand, we can obtain the so-called two-loop finite remainder by subtracting the UV and IR poles from the two-loop amplitude where the pole structure in Eq. (3.6) is expanded to O( 0 ). The analytic form of the finite remainder is in general much simpler than the finite part, as demonstrated in [9][10][11]65].

Local master integrals
The master integrals appearing in the two-loop leading colour qQQq ν and qggq ν amplitudes, that are not known analytically, are evaluated numerically using sector decomposition [66][67][68]. In general, it is challenging to obtain results with good numerical accuracy for complicated master integral topologies (e.g. topologies with 6, 7 or 8 propagators in Figs. 2 and 3) within a reasonable amount of time, even for an evaluation in the Euclidean region. Having numerical accuracy under control is particularly essential when large cancellations occur between different terms in the amplitude. One way in which this can be achieved is to use a basis of master integral with local numerator insertions [74,75] to regulate divergences. This is the approach we explore in this work. Another approach well suited to numerical evaluation is to use a quasi-finite basis of integrals [89,90]. We did not attempt to compare the two approaches but note that at least two factors are involved: firstly, in the reorganisation of the amplitude through the change of basis and secondly, the improved convergence of the resulting master integrals when evaluated with sector decomposition.
As an example, we consider one of the master integral topologies with 8 propagators, A 1 pentagon-box topology in Figure 1. There are three master integrals for this topology labelled according to Eq. (2.8). In Eq. (4.1), we use a notation for the integral with a numerator insertion N (k i , p i , µ ij ) where T is the diagram topology that we will specify by drawing it, and D α is a set of massless propagator denominator for a given topology T . The Laurent expansion for those master integrals starts at O( −4 ) for G 11111111000 and O( −3 ) for G 11111111−100 and G 111111110−10 , hence, getting accurate numerical results for the finite part is a demanding task. Here we follow the notation used in [91] where the local numerators were applied to six-gluon all-plus helicity amplitudes in Yang-Mills theory. For the master integral topology under consideration, we choose the following basis of master integral that exhibits improved IR behaviour tr where tr ± (ijkl) = 1 2 tr((1 ± γ 5 ) / p i / p j / p k / p l ). The first two integrals in Eq. (4.4) evaluate to O( ) and do not contribute to the two-loop amplitude, while the last integral is finite. Note that, for this topology the master integral coefficients do not contain any poles in , therefore, the master integrals need to be expanded to O(1). To evaluate the last integral in Eq. (4.4) using sector decomposition method, we first need to write the numerator insertion in terms of scalar products and momentum twistor variables (k i · k j , k i · p j , x i ) using loop momentum decomposition given in Eqs. f (k i · k j , k i · p j , x i ) .

(4.5)
The integral on the RHS can be directly evaluated using pySecDec [68] by passing the whole numerator into sector decomposition algorithm. For topologies with four point kinematics an extra stage of transverse integration is necessary to convert the integrals into a form compatible with the sector decomposition approach.
In Tables 1 and 2, we present a list of master integral topologies with local numerator insertions, where we use the following shorthand notation for the numerator insertions ) We also include in Tables 1 and 2 the order at which the expansion of the master integral starts, O(δ 1 ), the highest order in needed from the master integral for the amplitude evaluation, O(δ 2 ), as well as the list of possible permutations/crossings to obtain the full list of master integral topologies with local numerator insertions. δ 1 > δ 2 indicates that an integral does not contribute to the finite part of the amplitude, while δ 1 = δ 2 means that the integral contributes only to the finite part of the amplitude. In this latter case, it means only the leading order term in the expansion of the integral is required.  Table 1: Master integral topologies made up of five external legs with local numerator insertions. The topology T and numerator N (k i , p i , µ ij ) correspond to the integral definition in Eq. (4.3). Numerator building blocks Ψ, Φ and Ω are defined in Eq. (5). δ 1 is the order at which the expansion of the master integral starts, while δ 2 the highest order in needed from the master integral for the amplitude evaluation.   (5). δ 1 is the order at which the expansion of the master integral starts, while δ 2 the highest order in needed from the master integral for the amplitude evaluation.

Numerical results
We select a random Euclidean phase-space point by choosing rational values in the momentum twistor parametrisation from Eq. The numerical results of the leading colour partial amplitudes are obtained by evaluating the master integrals in three different ways: 1. We make use of the master integral solutions that are known analytically and readily available for evaluation in the Euclidean kinematics [6,69]. Tables 1 and 2, are numerically evaluated using pySecDec [68].  Note that there is a class of four-point master integrals with two off shell legs, that are not covered in [6], but available in [70][71][72], where the solutions are derived in the physical region. Instead of taking those results and perform analytic continuations to the Euclidean region for the numerical evaluation, we choose finite local master integral bases for the 6-and 7-propagator two-loop topologies, as shown in Table 2, and directly evaluate the rest of the integrals in this class numerically 3 . To assess uncertainties from the numerical evaluations via the sector decomposition method, we perform numerical integrations with three different random number seeds. The final results are obtained by taking the average and the error is computed by averaging the difference among the three results. We present the results for unrenormalised qQQq ν and qggq ν helicity amplitudes, normalised to the tree level amplitudê

Master integrals defined with the local numerator insertions, shown in
with helicities λ i . We can further split the amplitude into components of d s − 2 These integrals belong to the third type of evaluation discussed above. We see no problem in performing the analytic continuation on the expressions for the double off-shell 2 → 2 integrals. However, since the approach with local numerators and the sector decomposition worked with sufficient accuracy it was unnecessary to do this for our example.    In Tables 3 and 4 we display numerical evaluations of the helicity amplitudes for qggq ν and qQQq ν channels, respectively, using the kinematic point given in Eq. (5.1) for each (d s − 2) component. In Tables 5 and 6, we compare the divergent part of our numerical results against the universal two-loop pole structure in Eq. (3.6) in the HV scheme. To obtain the two-loop pole structure of Eq. (3.6) up to the single pole, we need to compute the one-loop amplitude up to O( ) for both the qggq ν and qQQq ν processes. The oneloop amplitude is computed by processing the Feynman diagrams through our integrand reduction setup, followed by IBP decomposition into the one-loop master integral basis consisting of six bubbles, a three-mass triangle, two one-mass boxes, a two-mass easy box, two two-mass hard boxes and a one-mass pentagon, in a similar fashion to the two-loop case that is discussed in Section 2. The O( ) part of the two-mass easy box and one-mass pentagon integrals are evaluated numerically using Fiesta/pySecDec, while the rest are obtained from available analytic expressions [92,93]. Therefore, the numerical values quoted for the poles in Tables 5 and 6 are exact up to O( −2 ). We assess the uncertainty of the O( −1 ) part of the pole structure using the same method as in the two-loop numerical evaluations.
We additionally perform another check of our results by independently processing the two-loop Feynman diagrams through a numerical diagram-based integrand reduction into an integrand representation consisting of four-dimensional ISPs and the extra-dimensional part of the loop momenta µ ij . The integrals containing µ ij are first written as dimensionshifted integrals and all integrals appearing in the amplitude are evaluated directly using Fiesta/pySecDec. This approach is similar to the method we employed for the numerical evaluation of the planar two-loop five gluon amplitude [58]. The results obtained using this method are in perfect agreement with the results reported in this paper.

Conclusions
In this article we have presented numerical results for the planar two-loop helicity amplitudes for the scattering of a W -boson with four partons for the first time. This computation is the first step towards obtaining analytic expressions using the reconstruction of rational functions with finite field arithmetic.
A number of important steps remain to be completed in order for this to be a feasible target. Firstly, the complete list of master integrals should be evaluated analytically since the sector decomposition approach is still too CPU intensive for phenomenological applications. This seems a reasonable aim owing to the recent success of the planar and non-planar pentagon functions [32-34, 62-64, 94], though the efficient numerical implementation of the resulting analytic functions will require further study. Secondly, the coefficients of the master integrals still have a high degree of algebraic complexity. As shown in applications to five-parton scattering, direct reconstruction of the finite remainder after subtraction of UV and IR poles leads to a substantial reduction in complexity [9][10][11].
We hope that the work presented here presents valuable information that can be used to achieve these goals, as well as providing encouragement that they are realistic in the near future.