Two-loop QCD corrections to $Wb\bar{b}$ production at hadron colliders

We present an analytic computation of the two-loop QCD corrections to $u\bar{d}\to W^+b\bar{b}$ for an on-shell $W$-boson using the leading colour and massless bottom quark approximations. We perform an integration-by-parts reduction of the unpolarised squared matrix element using finite field reconstruction techniques and identify an independent basis of special functions that allows an analytic subtraction of the infrared and ultraviolet poles. This basis is valid for all planar topologies for five-particle scattering with an off-shell leg.


INTRODUCTION
The production of a W -boson in association with a pair of b-quarks at hadron colliders is of fundamental importance as a background to Higgs production in association with a vector boson. The process is one of a prioritised list of 2 → 3 scattering problems for which higher order corrections are necessary to keep theory in line with data. These amplitudes are related to a large class of processes contributing to pp → W + 2j production and the work presented here represents a significant step towards achieving a complete classification of the missing two-loop amplitudes.
The process has been studied extensively at next-toleading order (NLO) [1][2][3][4][5] and was the first in a set of off-shell five-particle amplitudes to be studied using the unitarity method [6,7]. The present state of the art in phenomenological studies allows full mass effects, shower matching, electro-weak corrections and the inclusion additional QCD jets [8][9][10].
A numerical computation of the two-loop helicity amplitudes [11] demonstrated the importance of an efficient analytic form with a well understood basis of special functions. Major steps forward came via efficient numerical evaluation of the differential equations [12] and analytic evaluation in terms the Goncharov Polylogarithms (GPLs) [13,14]. These results opened the door for a fully analytic amplitude computation yet significant challenges remain. The complexity of the external kinematics represents a challenge for integral reduction techniques and the identification of a minimal basis of special functions is required to find analytic simplifications after subtracting universal infrared and ultraviolet divergences.
In this short letter we outline the extension of this method to processes with an additional mass scale.
The leading order process consists of two simple Feynman diagrams as shown in Fig. 1. We label our process as follows, where p 2 1 = p 2 2 = p 2 3 = p 2 4 = 0 and p 2 5 = m 2 W . The colour decomposition at leading colour is 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 interfere the L-loop partial amplitudes A (L) in Eq. (2) with the tree-level partial amplitude A (0) to obtain the unrenormalised L-loop unpolarised squared partial amplitude, After the interference with the tree-level amplitude the analytic expression can be written in terms of scalar invariants, and a parity-odd quantity, tr 5 = 4i µνρσ p µ 1 p ν 2 p ρ 3 p σ 4 . Our results are the so-called finite remainders F (L) , obtained after subtraction of infrared and ultraviolet divergences, where P (L) takes the well known form [50][51][52][53]. The explicit form for our process using the same conventions can be found in Ref. [11].

AMPLITUDE REDUCTION
Feynman diagrams for the ud → W + bb scattering are generated using Qgraf [54]. In the leading colour approximation, there are 2, 16 and 210 diagrams contributing to the tree level, 1-loop and 2-loop amplitudes, respectively.
where p are the external momenta which live in four dimensions, and k i are the loop momenta. We work in the conventional dimensional regularisation (CDR) scheme, where we have d = 4 − 2 dimensions.
The W -boson polarisation sum is performed in the unitary gauge, The terms containing traces with a single γ 5 are treated using Larin's prescription [55], while those with two γ 5 's are treated using the anti-commuting γ 5 prescription. Larin's prescription has been employed in a wide variety of multi-loop computations and a detailed discussion can be found, for example, in Ref. [56]. We have checked  that using Larin's scheme throughout gives the same results for F (L) . We can then split the squared partial amplitude into parity-even and parity-odd parts, M (L) even receives contribution from the terms with no or two γ 5 's while tr 5 M (L) odd is made up of terms with a single γ 5 . We note that the parity-odd part vanishes at tree level, M To perform the reduction of the 2-loop amplitude onto a basis of master integrals we first map each topology T to a set of 15 maximal cut or master topologies as shown in Fig. 3. The master topologies are then defined with a spanning set of 11 propagators and, after tracking shifts in the loop momentum, the change of variables for each topology T can be computed. The resulting squared partial amplitude is now written as a linear combination of scalar integrals I where k ∈ {even, odd}. Analytic forms of the unreduced squared matrix elements above are derived using a collection of Form [57,58] and Mathematica routines. The integrals appearing in Eq. (8) are not all independent.
Relations between integrals I can be found using IBP identities and the squared amplitude can be written in terms of an independent set of master integrals as follows The reduction to master integral basis is then performed within the FiniteFlow framework [23], separately for M (2) even and M odd . We use LiteRed [59] to generate the IBP relations in Mathematica, together with the Laporta algorithm [60] to solve them numerically over finite fields. We note that only master topologies T 1 − T 10 are included in the IBP system since the integrals belonging to master topologies T 11 − T 15 can be mapped onto master topologies T 6 −T 10 . The procedure for performing the reduction onto master integrals using IBP relations is of course extremely well known, the challenge in this example is one of enormous algebraic complexity. By encoding the problem within a numeric sampling modular arithmetic we are able to efficiently solve the Laporta system with tensor integral ranks of up to five, avoiding all large intermediate expressions. For planar topologies such as the ones appearing here the application of syzygy relations [61][62][63] to optimise the IBP reduction would likely lead to a substantial speed-up in computation time, although in our case it was not found to be necessary. We did not perform an analytic reconstruction after completing the set up of the reduction in FiniteFlow graphs. Instead we continued to map the amplitude onto a basis of special functions.

A BASIS OF SPECIAL FUNCTIONS FOR THE FINITE REMAINDER
There are 202 master integrals contributing to the amplitude, 196 of them are covered by the 3 independent pentabox master integral topologies, while 6 are of oneloop squared type that involve one-loop massive on-shell bubble integral. We choose the canonical bases of master integrals constructed in Ref. [12]. They satisfy differential equations (DEs) [64][65][66][67] in the canonical form [68], where − → MI is the set of canonical master integrals for any of the involved topologies, the a i are constant rational matrices, while {w i } 58 i=1 is a set of algebraic functions of the external kinematics called letters (see Ref. [12] for their definition). The alphabet, i.e. the set of all letters, is the same for all planar one-mass five-particle integrals up to two loops, whereas the constant matrices a i depend on the topology. In Ref. [12], the authors also discuss a strategy to evaluate the master integrals numerically, based on the solution of the DEs (10) in terms of generalised power series [69]. More recently, analytic expressions of the canonical master integrals in terms of GPLs [70][71][72] have become available [13,14]. Both approaches allow for the numerical evaluation of the master integrals in any kinematic region and with arbitrary precision. Both approaches, however, also share certain drawbacks. Whether we reconstruct the prefactors of the -components of the master integrals in Eq. (9) or we map the latter onto monomials of GPLs, we cannot subtract the infrared and ultraviolet poles analytically and reconstruct directly the finite remainder.
We overcome these issues by constructing a basis out of the -components of the canonical master integrals up to order 4 . The crucial tool we employ in this construction are Chen's iterated integrals [73]. We can define them iteratively through where s denotes cumulatively the kinematic invariants, s 0 is an arbitrary boundary point, and the iteration starts from [] s0 (s) = 1. The depth n of the iterated integral is called transcendental weight. We refer to the notes [74] for a thorough discussion. All GPLs can be rewritten in terms of iterated integrals. The latter however offer two useful advantages. The first is that -conjecturally -they implement automatically all the functional relations. Once a GPL expression is rewritten in terms of iterated integrals in a given alphabet {w i }, finding the functional relations becomes a linear algebra problem, as 'words' [w i1 , . . . , w in ] with different letters are linearly independent. The second is that it is completely straightforward to write out the solution of the canonical DEs (10) in terms of iterated integrals. Eq. (10) in fact implies the following differential relation between consecutive components of the expansion of the master integrals, where − → MI (k) is the O( k ) term of the master integrals. Comparing Eq. (12) to Eq. (11), we see that the iterated integral expressions of − → MI (k) are obtained by adding a letter to the right of those of the previous order, multiplying them by the constant matrices a i , and adding the boundary values. The master integrals are normalised to start from O( 0 ) and so the O( k ) components have transcendental weight k.
We used the GPL expressions of Refs. [13,14] to compute the values of the master integrals in an arbitrary point s 0 with 1100-digit precision. Using the PSLQ algorithm [75], we determined the integer relations among the boundary values, and rewrote them in terms of a basis of transcendental constants. Next, we used the differential equations provided by Ref. [12] to express the relevant master integrals in terms of iterated integrals. This allowed us to determine a minimal set of linearly independent integral components, order by order in up to 4 . We denote these functions by {f corresponds to an component of the master integrals, we can evaluate them numerically using the methods of Refs. [12][13][14], with the additional advantage that they are linearly independent.
In order to subtract the poles analytically, we need to be able to write in the same basis also the subtraction term. From the transcendental point of view, the latter is given by the product of certain logarithms and transcendental constants coming from the anomalous dimensions -π 2 and ζ 3 -times the one-loop amplitude. In order to accommodate this in the basis, we add the transcendental constants as elements, and work out the relations between the functions at each weight and products of lower-weight ones using the shuffle algebra of the iterated integrals. As a result, the functions in the basis {f (w) i } are indecomposable, i.e. they cannot be rewritten in terms of lower-weight elements of the basis.
Armed with this function basis, we can proceed with the reconstruction of the two-loop finite remainders. We map the master integrals appearing in Eq. (9) onto a monomial basis of the functions {f (w) i }, which we denote by {m(f )}, and perform a Laurent expansion in up to O( 0 ). We do the same for the subtraction term P (2) . The resulting finite remainder, is indeed free of poles. We set s 12 = 1 to simplify the reconstruction. The dependence can be recovered a posteriori through dimensional analysis. The coefficients are not all independent. We find the linear relations between them and a set of additional coefficients we supply as ansatz. We used tree-level expressions, coefficients from the one-loop amplitude and from the unreduced scalar integrals. Through these linear relations we rewrite the complicated coefficients in F (2) in terms of known coefficients from the ansatz and simpler ones, which finally have to be reconstructed. Moreover, we simplify the reconstruction of the remaining coefficients by partial fractioning them with respect to s 23 . First we determine the denominator factors by computing a univariate slice and matching it against an ansatz made of letters w i . Using the information about the denominator and the polynomial degree in the numerator, we construct an ansatz for the partial-fractioned expressions of the coefficients. Then we fit the ansatz with a numerical sampling. See Refs. [38,76,77] for recent work on multivariate partial fractioning. To emphasize the effectiveness of our strategy, we note that the coefficients of the parity-even (-odd) two-loop amplitude written in terms of GPL monomials have maximal degree 62 (63). The maximal degree drops to 54 (54) when we use the basis of special functions {f (w) i } in the finite remainder, and then to 31 (32) in the remaining 4 variables after partial fractioning. The reconstruction finally required 38663 (45263) sample points over 2 prime fields, gaining a factor of 7 in the reconstruction time with respect to the GPL-based approach [78]. The reconstructed analytic expressions are further simplified using the Multi-variateApart package [77].
The iterated integrals expression of the f (w) i functions allow us to study the analytic structure of the finite re-mainder in a very convenient way. Interestingly, we observe that certain letters do not appear. As it was already noted in Ref. [12], the last 9 letters do not show up in any two-loop amplitude up to order 0 . Out of the relevant 49 letters, 6 (w i with i ∈ {16, 17, 27, 28, 29, 30}) appear in the master integrals but cancel out in the two-loop amplitude truncated at O( 0 ). Finally, the letter w 49 = tr 5 is present in the two-loop amplitude, but cancels out in the finite remainder. This letter has already been observed to exhibit the same behaviour in all the known massless two-loop five-particle amplitudes [36][37][38][39][40][41][42][43][44][45]79], which has spawned interest in the context of cluster algebras [80].
As regards the numerical evaluation, we propose a strategy based on the generalised power series solution of the DEs [69] applied not to the master integrals, but directly to the basis of special functions. If we rescale each special function in the basis f , 1} satisfies a system of DEs in the canonical form (10). This follows from the differential property of the iterated integrals (11). Differently from the DEs for the master integrals, the DEs for the special functions contain only the minimal amount of information necessary to evaluate the finite remainder. For instance, instead of evaluating all the weight-4 functions that may appear in any one-mass two-loop five-particle amplitude, we can restrict ourselves to evaluating only the 19 linear combinations that actually appear in F i 's, and so on down to weight zero. The resulting DEs are much simpler that those for the master integrals. For instance, they are by-construction free of the letters which do not appear in the finite-remainder, and their dimension is smaller than the number of master integrals for all the relevant families. Finally, we evaluate the g (w) i 's by solving the corresponding DEs using the Mathematica package DiffExp [81]. We compute the boundary values in an arbitrary point in the physical scattering region through the correspondence between the g (w) i 's and the master integral components.
The complete analytic expression of the two-loop finite remainder in terms of rational coefficients and special functions is included in the ancillary files, together with the differential equation and the boundary values necessary to evaluate the latter numerically. We performed Ward identity checks at the level of master integrals for M (2) even and at the level of the finite remainder for M (2) odd : we modified the numerator functions by replacing the loop and tree-level amplitude polarisation vectors with p 5 and p 1 respectively, and observed that M (2) even and F (2) odd vanish. We also compared numerically the finite remain- k , and finite remainder, F k , normalised to the tree level squared partial amplitude in 4 dimensions, M ders derived in this work against results from an independent helicity amplitude computation in the t'Hooft-Veltman scheme using the framework of Ref. [11]. For the convenience of future cross-checks, we provide the numerical values of M  Table I.

DISCUSSION AND OUTLOOK
The results we have obtained represent a major step forward and open the door to phenomenological applications. The identification of a basis of special functions has resulted in a substantial speed up over previous studies as well as uncovering explicit cancellations and reduction in complexity. To demonstrate the suitability for phenomenological applications we present the evaluation on a univariate slice of the physical phase space. For this we use a parametrisation in terms of energy fractions and angles of the final state, where p 1 and p 2 are taken back-to-back along the z-axis with a total centre-of-mass energy of s. We have chosen p 3 to be produced at an elevation of π 2 from the z-axis and the on-shell phase space conditions impose cos θ = 1 + 2 In Fig. 4 we plot values of the one-and two-loop finite remainders against x 2 for a configuration with φ = 0.1, m W = 0.1, s = 1 and x 1 = 0.6. The special functions were evaluated with DiffExp [81] using rationalized values of the invariants. An average evaluation time of 260s over 1000 points was observed and the function is smooth and stable over the whole region. This demonstrates that even with a basic setup in Mathematica a reasonable evaluation time can be achieved and that realistic phenomenology can now be performed.
The results obtained here pave the way for a broader class of 2 → 3 scattering problems. The solution of the IBP system and the basis of special functions do not depend on the on-shell approximation of the W -boson and apply equally to the planar sectors of pp → W/Z +2j (including decays) and pp → H + 2j. Going beyond leading colour for pp → W/Z + 2j or any complete pp → H + 2j amplitudes at two-loops still requires missing information on the non-planar master integrals, nevertheless we believe they can be easily incorporated into the strategy we introduce here.
We thank Herschel Chawdhry, Thomas Gehrmann, Johannes Henn, Alexander Mitov, Tiziano Peraro and Rene Poncelet for numerous insightful discussions. We also thank Nikolaos Syrrakos for kindly providing the results of Ref. [14] prior to its publication. This project has received funding from the European Union's Horizon 2020 research and innovation programmes New level of theoretical precision for LHC Run 2 and beyond (grant agreement No 683211), High precision multi-jet dynamics at the LHC (grant agreement No 772009), and Novel structures in scattering amplitudes (grant agreement No 725110). HBH has been partially supported by STFC consolidated HEP theory grant ST/T000694/1. SZ gratefully acknowledges the computing resources provided by the Max Planck Institute for Physics.