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

Simon Badger<EMAIL_ADDRESS>Dipartimento di Fisica and Arnold-
Regge Center, Università di Torino, and INFN, Sezione di Torino, Via P. Giuria
1, I-10125 Torino, Italy  Heribertus Bayu Hartanto
<EMAIL_ADDRESS>Cavendish Laboratory, University of Cambridge,
Cambridge CB3 0HE, United Kingdom  Simone Zoia<EMAIL_ADDRESS>Dipartimento di Fisica and Arnold-Regge Center, Università di Torino, and
INFN, Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy

###### Abstract

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.

††preprint: CAVENDISH-HEP-21/01

##  I 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\to 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\to 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-to-leading 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.

Efficient amplitude and integration-by-parts reduction (IBP) [15, 16] using
finite field arithmetic [17, 18, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27] has
gained significant interest in recent years. Through multiple evaluations of a
numerical algorithm [28, 29, 30, 31], fully analytic forms for planar massless
five-particle amplitudes have been extracted using a rational parametrisation
of the kinematics [32]. Following a complete understanding of a pentagon
function basis [33, 34], a large number of two-loop amplitudes are now
available in compact analytic form [35, 36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47]. We have also seen the first phenomenological predictions at NNLO
in QCD for the production of three photons in hadron colliders after
combination with real-virtual and double real radiation [48, 49].

In this short letter we outline the extension of this method to processes with
an additional mass scale.

##  II Leading colour $u\bar{d}\to W^{+}b\bar{b}$ amplitudes

Figure 1: Leading order Feynman diagrams contributing to $u\bar{d}\to
W^{+}b\bar{b}$ . Figure 2: Sample Feynman diagrams in the leading colour two-
loop $u\bar{d}\to W^{+}b\bar{b}$ amplitude.

The leading order process consists of two simple Feynman diagrams as shown in
Fig. 1. We label our process as follows,

$\bar{d}(p_{1})+u(p_{2})\rightarrow b(p_{3})+\bar{b}(p_{4})+W^{+}(p_{5}),$ (1)

where $p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=p_{4}^{2}=0$ and $p_{5}^{2}=m_{W}^{2}$.
The colour decomposition at leading colour is

$\displaystyle\mathcal{A}^{(L)}(1_{\bar{d}},2_{u},3_{b},4_{\bar{b}},5_{W})=$
$\displaystyle\qquad
n^{L}g_{s}^{2}g_{W}\delta_{i_{1}}^{\;\;\bar{i}_{4}}\delta_{i_{3}}^{\;\;\bar{i}_{2}}A^{(L)}(1_{\bar{d}},2_{u},3_{b},4_{\bar{b}},5_{W}),$
(2)

where $n=m_{\epsilon}N_{c}\alpha_{s}/(4\pi),\ \alpha_{s}=g_{s}^{2}/(4\pi)$ and
$m_{\epsilon}=i(4\pi)^{\epsilon}e^{-\epsilon\gamma_{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,

$\displaystyle M^{(L)}=\sum_{\mathrm{spin}}A^{(0)*}A^{(L)}.$ (3)

After the interference with the tree-level amplitude the analytic expression
can be written in terms of scalar invariants,

$\displaystyle s_{12}=(p_{1}+p_{2})^{2}\,,\quad
s_{23}=(p_{2}-p_{3})^{2}\,,\quad s_{34}=(p_{3}+p_{4})^{2}\,,$ $\displaystyle
s_{45}=(p_{4}+p_{5})^{2}\,,\quad s_{15}=(p_{1}-p_{5})^{2}\,,\quad
s_{5}=p_{5}^{2}\,,$ (4)

and a parity-odd quantity,
$\mathrm{tr}_{5}=4i\epsilon_{\mu\nu\rho\sigma}p_{1}^{\mu}p_{2}^{\nu}p_{3}^{\rho}p_{4}^{\sigma}$.
Our results are the so-called finite remainders $F^{(L)}$, obtained after
subtraction of infrared and ultraviolet divergences,
$F^{(L)}=M^{(L)}-P^{(L)}$, 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].

##  III Amplitude reduction

Feynman diagrams for the $u\bar{d}\rightarrow W^{+}b\bar{b}$ 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. Example 2-loop diagrams are shown in Figure 2. Upon
interference of the $L$-loop partial amplitude $A^{(L)}$ with the tree level
partial amplitude $A^{(0)}$ according to Eq. (3), the squared partial
amplitude can be written as

$M^{(L)}(\\{p\\})=\int\prod_{i=1}^{L}\frac{d^{d}k_{i}}{i\pi^{d/2}e^{-\epsilon\gamma_{E}}}\sum_{T}\frac{N_{T}(d,\\{k\\},\\{p\\})}{\prod_{\alpha\in
T}D_{\alpha}(\\{k\\},\\{p\\})},$ (5)

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\epsilon$ dimensions.

The $W$-boson polarisation sum is performed in the unitary gauge,

$i\sum_{\lambda}\varepsilon_{W}^{\mu*}(p_{5},\lambda)\varepsilon_{W}^{\nu}(p_{5},\lambda)=-g^{\mu\nu}+\frac{p_{5}^{\mu}p_{5}^{\nu}}{m_{W}^{2}}.$
(6)

The terms containing traces with a single $\gamma_{5}$ are treated using
Larin’s prescription [55], while those with two $\gamma_{5}$’s are treated
using the anti-commuting $\gamma_{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)}(\\{p\\})=M_{\mathrm{even}}^{(L)}(\\{p\\})+\mathrm{tr}_{5}M_{\mathrm{odd}}^{(L)}(\\{p\\}).$
(7)

$M_{\mathrm{even}}^{(L)}$ receives contribution from the terms with no or two
$\gamma_{5}$’s while $\mathrm{tr}_{5}M_{\mathrm{odd}}^{(L)}$ is made up of
terms with a single $\gamma_{5}$. We note that the parity-odd part vanishes at
tree level, $M^{(0)}_{\mathrm{odd}}=0$.

Figure 3: Topologies with maximum number of propagators.

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 $\mathcal{I}$

$M_{k}^{(2)}(\\{p\\})=\sum_{i}c_{k,i}(\epsilon,\\{p\\})\;\mathcal{I}_{k,i}(\epsilon,\\{p\\}),$
(8)

where $k\in\\{\mathrm{even},\mathrm{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 $\mathcal{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

$M_{k}^{(2)}(\\{p\\})=\sum_{i}d_{k,i}(\epsilon,\\{p\\})\;\mathrm{MI}_{i}(\epsilon,\\{p\\}).$
(9)

The reduction to master integral basis is then performed within the FiniteFlow
framework [23], separately for $M_{\mathrm{even}}^{(2)}$ and
$M_{\mathrm{odd}}^{(2)}$. 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.

##  IV 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 one-loop 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],

$\displaystyle d\overrightarrow{\text{MI}}=\epsilon\sum_{i=1}^{58}a_{i}d\log
w_{i}\,\overrightarrow{\text{MI}}\,,$ (10)

where $\overrightarrow{\text{MI}}$ is the set of canonical master integrals
for any of the involved topologies, the $a_{i}$ are constant rational
matrices, while $\\{w_{i}\\}_{i=1}^{58}$ 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 $\epsilon$-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
$\epsilon$-components of the canonical master integrals up to order
$\epsilon^{4}$. The crucial tool we employ in this construction are Chen’s
iterated integrals [73]. We can define them iteratively through

$\displaystyle d[w_{i_{1}},\ldots,w_{i_{n}}]_{s_{0}}(s)=d\log
w_{i_{n}}[w_{i_{1}},\ldots,w_{i_{n-1}}]_{s_{0}}(s)\,,$ (11)
$\displaystyle[w_{i_{1}},\ldots,w_{i_{n}}]_{s_{0}}(s_{0})=0\,,$

where $s$ denotes cumulatively the kinematic invariants, $s_{0}$ is an
arbitrary boundary point, and the iteration starts from $[]_{s_{0}}(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_{i_{1}},\ldots,w_{i_{n}}]$
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 $\epsilon$
expansion of the master integrals,

$\displaystyle d\overrightarrow{\text{MI}}^{(k)}=\sum_{i=1}^{58}a_{i}d\log
w_{i}\,\overrightarrow{\text{MI}}^{(k-1)}\,,\quad\forall k\geq 1\,,$ (12)

where $\overrightarrow{\text{MI}}^{(k)}$ is the $\mathcal{O}(\epsilon^{k})$
term of the master integrals. Comparing Eq. (12) to Eq. (11), we see that the
iterated integral expressions of $\overrightarrow{\text{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
$\mathcal{O}(\epsilon^{0})$ and so the $\mathcal{O}(\epsilon^{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 $\epsilon$ up to $\epsilon^{4}$. We denote these
functions by $\\{f^{(w)}_{i}\\}$, where $w=1,\ldots,4$ labels the weight.
Since each $f^{(w)}_{i}$ corresponds to an $\epsilon$ 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 – $\pi^{2}$ and
$\zeta_{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 $\epsilon$ up to
$\mathcal{O}(\epsilon^{0})$. We do the same for the subtraction term
$P^{(2)}$. The resulting finite remainder,

$F_{k}^{(2)}(\\{p\\})=\sum_{i}e_{k,i}(\\{p\\})\;m_{k,i}(f)+\mathcal{O}(\epsilon)\,,$
(13)

is indeed free of $\epsilon$ poles. We set $s_{12}=1$ to simplify the
reconstruction. The dependence can be recovered a posteriori through
dimensional analysis. The coefficients $e_{k,i}(\\{p\\})$ in $F_{k}^{(2)}$ 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 111Since estimates
of evaluation time rely on system specific parameters, we have taken a
conservative value for the speed improvement.. The reconstructed analytic
expressions are further simplified using the MultivariateApart package [77].

The iterated integrals expression of the $f^{(w)}_{i}$ functions allow us to
study the analytic structure of the finite remainder 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 $\epsilon^{0}$. Out of the relevant 49 letters, 6
($w_{i}$ with $i\in\\{16,17,27,28,29,30\\}$) appear in the master integrals
but cancel out in the two-loop amplitude truncated at
$\mathcal{O}(\epsilon^{0})$. Finally, the letter $w_{49}=\mathrm{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, 41, 42, 43, 44, 38,
45, 79, 39, 40], 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^{(w)}_{i}$ by a power of $\epsilon$
corresponding to its weight, the ensuing list of functions
$\vec{v}=\\{\epsilon^{w}f^{(w)}_{i},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_{k}^{(2)}$. We therefore
define a new basis of special functions, $\\{g^{(w)}_{i}\\}$, which at weight
four includes the aforementioned 19 combinations of $f^{(4)}_{i}$’s, at weight
three contains only the $f^{(3)}_{i}$’s appearing in $F_{k}^{(2)}$ and in the
derivatives of the $g^{(4)}_{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_{\mathrm{even}}^{(2)}$ and at
the level of the finite remainder for $M_{\mathrm{odd}}^{(2)}$: 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_{\mathrm{even}}^{(2)}$ and $F_{\mathrm{odd}}^{(2)}$ vanish. We also
compared numerically the finite remainders 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_{k}^{(2)}$ and
$F_{k}^{(2)}$ at one phase space point in Table 1.

|  Re$\\{M_{\mathrm{even}}^{(2)}/M^{(0)}_{\epsilon=0}\\}$ |  Re$\\{\mathrm{tr}_{5}M_{\mathrm{odd}}^{(2)}/M^{(0)}_{\epsilon=0}\\}$
---|---|---
$\epsilon^{-4}$ |  2 |  0
$\epsilon^{-3}$ |  -2.19718713546 |  0
$\epsilon^{-2}$ |  -12.7892676147 |  -0.211614995129
$\epsilon^{-1}$ |  -7.77698255746 |  11.5990058250
$\epsilon^{0}$ |  116.073111075 |  28.7730449523
|  Re$\\{F_{\mathrm{even}}^{(2)}/M^{(0)}_{\epsilon=0}\\}$ |  Re$\\{\mathrm{tr}_{5}F_{\mathrm{odd}}^{(2)}/M^{(0)}_{\epsilon=0}\\}$
$\epsilon^{0}$ |  144.141227186 |  -7.34974777490
Table 1:  Numerical results for the leading colour two-loop squared partial
amplitude, $M_{k}^{(2)}$, and finite remainder, $F_{k}^{(2)}$, normalised to
the tree level squared partial amplitude in 4 dimensions,
$M^{(0)}_{\epsilon=0}$, at the physical point
$\\{s_{12}=5,s_{23}=-1/3,s_{34}=11/13,s_{45}=17/19,s_{15}=-23/29,s_{5}=1/7\\}$.

##  V 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,

$\displaystyle p_{3}$
$\displaystyle=\tfrac{x_{1}\sqrt{s}}{2}\left(1,1,0,0\right)\,,$ (14)
$\displaystyle p_{4}$
$\displaystyle=\tfrac{x_{2}\sqrt{s}}{2}\left(1,\cos\theta,-\sin\phi\sin\theta,-\cos\phi\sin\theta\right)\,,$
$\displaystyle p_{5}$
$\displaystyle=\sqrt{s}\left(1,0,0,0\right)-p_{3}-p_{4}\,,$

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 $\tfrac{\pi}{2}$ from the $z$-axis and the on-shell phase
space conditions impose
$\cos\theta=1+\tfrac{2}{x_{1}x_{2}}\left(1-x_{1}-x_{2}-\tfrac{m_{W}^{2}}{s}\right)$.
In Fig. 4 we plot values of the one- and two-loop finite remainders against
$x_{2}$ for a configuration with $\phi=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 $260$s 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.

Figure 4: The finite remainder
$F^{(L)}=F_{\mathrm{even}}^{(L)}+\mathrm{tr}_{5}F_{\mathrm{odd}}^{(L)}$ for
$u\bar{d}\to W^{+}b\bar{b}$ at one and two loops as a function of the variable
$x_{2}$ defined in Eq. (14).

The results obtained here pave the way for a broader class of $2\to 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\to W/Z+2j$ (including decays) and
$pp\to H+2j$. Going beyond leading colour for $pp\to W/Z+2j$ or any complete
$pp\to 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.

###### Acknowledgements.

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.

## References

  * Ellis and Veseli [1999] R. K. Ellis and S. Veseli, Strong radiative corrections to W b anti-b production in p anti-p collisions, Phys. Rev. D 60, 011501 (1999), arXiv:hep-ph/9810489 .
  * Febres Cordero  _et al._ [2009] F. Febres Cordero, L. Reina, and D. Wackeroth, W- and Z-boson production with a massive bottom-quark pair at the Large Hadron Collider, Phys. Rev. D 80, 034015 (2009), arXiv:0906.1923 [hep-ph] .
  * Badger  _et al._ [2011] S. Badger, J. M. Campbell, and R. K. Ellis, QCD Corrections to the Hadronic Production of a Heavy Quark Pair and a W-Boson Including Decay Correlations, JHEP 03, 027, arXiv:1011.6647 [hep-ph] .
  * Frederix  _et al._ [2011] R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, R. Pittau, and P. Torrielli, W and $Z/\gamma*$ boson production in association with a bottom-antibottom pair, JHEP 09, 061, arXiv:1106.6019 [hep-ph] .
  * Oleari and Reina [2011] C. Oleari and L. Reina, W +- b $\bar{b}$ production in POWHEG, JHEP 08, 061, [Erratum: JHEP 11, 040 (2011)], arXiv:1105.4488 [hep-ph] .
  * Bern  _et al._ [1997] Z. Bern, L. J. Dixon, D. A. Kosower, and S. Weinzierl, One loop amplitudes for e+ e- —$>$ anti-q q anti-Q Q, Nucl. Phys. B 489, 3 (1997), arXiv:hep-ph/9610370 .
  * Bern  _et al._ [1998] Z. Bern, L. J. Dixon, and D. A. Kosower, One loop amplitudes for e+ e- to four partons, Nucl. Phys. B 513, 3 (1998), arXiv:hep-ph/9708239 .
  * Luisoni  _et al._ [2015] G. Luisoni, C. Oleari, and F. Tramontano, $Wb\overline{b}j$ production at NLO with POWHEG+MiNLO, JHEP 04, 161, arXiv:1502.01213 [hep-ph] .
  * Kallweit  _et al._ [2015] S. Kallweit, J. M. Lindert, P. Maierhöfer, S. Pozzorini, and M. Schönherr, NLO electroweak automation and precise predictions for W+multijet production at the LHC, JHEP 04, 012, arXiv:1412.5157 [hep-ph] .
  * Anger  _et al._ [2018] F. R. Anger, F. Febres Cordero, H. Ita, and V. Sotnikov, NLO QCD predictions for $Wb\bar{b}$ production in association with up to three light jets at the LHC, Phys. Rev. D 97, 036018 (2018), arXiv:1712.05721 [hep-ph] .
  * Hartanto  _et al._ [2019] H. B. Hartanto, S. Badger, C. Brønnum-Hansen, and T. Peraro, A numerical evaluation of planar two-loop helicity amplitudes for a W-boson plus four partons, JHEP 09, 119, arXiv:1906.11862 [hep-ph] .
  * Abreu  _et al._ [2020a] S. Abreu, H. Ita, F. Moriello, B. Page, W. Tschernow, and M. Zeng, Two-Loop Integrals for Planar Five-Point One-Mass Processes, JHEP 11, 117, arXiv:2005.04195 [hep-ph] .
  * Canko  _et al._ [2020] D. D. Canko, C. G. Papadopoulos, and N. Syrrakos, Analytic representation of all planar two-loop five-point Master Integrals with one off-shell leg, (2020), arXiv:2009.13917 [hep-ph] .
  * Syrrakos [2020] N. Syrrakos, Pentagon integrals to arbitrary order in the dimensional regulator, (2020), arXiv:2012.10635 [hep-ph] .
  * Tkachov [1981] F. V. Tkachov, A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions, Phys. Lett. 100B, 65 (1981).
  * Chetyrkin and Tkachov [1981] K. G. Chetyrkin and F. V. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl. Phys. B192, 159 (1981).
  * Wang [1981] P. S. Wang, A p-adic algorithm for univariate partial fractions, in  _Proceedings of the Fourth ACM Symposium on Symbolic and Algebraic Computation_, SYMSAC ’81 (ACM, New York, NY, USA, 1981) pp. 212–217.
  * Wang  _et al._ [1982] P. S. Wang, M. J. T. Guy, and J. H. Davenport, P-adic reconstruction of rational numbers, SIGSAM Bull. 16, 2 (1982).
  * Tra [2006] _ISSAC ’06: Proceedings of the 2006 International Symposium on Symbolic and Algebraic Computation_ (ACM, New York, NY, USA, 2006) 505060.
  * Kauers [2008] M. Kauers, Fast Solvers for Dense Linear Systems, Nuclear Physics B (Proc. Suppl.) 183, 245 (2008).
  * von Manteuffel and Schabinger [2015] A. von Manteuffel and R. M. Schabinger, A novel approach to integration by parts reduction, Phys. Lett. B744, 101 (2015), arXiv:1406.4513 [hep-ph] .
  * Peraro [2016] T. Peraro, Scattering amplitudes over finite fields and multivariate functional reconstruction, JHEP 12, 030, arXiv:1608.01902 [hep-ph] .
  * Peraro [2019] T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs, JHEP 07, 031, arXiv:1905.08019 [hep-ph] .
  * Smirnov and Chuharev [2020] A. V. Smirnov and F. S. Chuharev, FIRE6: Feynman Integral REduction with Modular Arithmetic, Comput. Phys. Commun. 247, 106877 (2020), arXiv:1901.07808 [hep-ph] .
  * Klappert and Lange [2020] J. Klappert and F. Lange, Reconstructing rational functions with FireFly, Comput. Phys. Commun. 247, 106951 (2020), arXiv:1904.00009 [cs.SC] .
  * Klappert  _et al._ [2020] J. Klappert, F. Lange, P. Maierhöfer, and J. Usovitsch, Integral Reduction with Kira 2.0 and Finite Field Methods, (2020), arXiv:2008.06494 [hep-ph] .
  * Caola  _et al._ [2020] F. Caola, A. von Manteuffel, and L. Tancredi, Di-photon amplitudes in three-loop Quantum Chromodynamics, (2020), arXiv:2011.13946 [hep-ph] .
  * Badger  _et al._ [2018a] S. Badger, C. Brønnum-Hansen, H. B. Hartanto, and T. Peraro, First look at two-loop five-gluon scattering in QCD, Phys. Rev. Lett. 120, 092001 (2018a), arXiv:1712.02229 [hep-ph] .
  * Abreu  _et al._ [2018a] S. Abreu, F. Febres Cordero, H. Ita, B. Page, and M. Zeng, Planar Two-Loop Five-Gluon Amplitudes from Numerical Unitarity, Phys. Rev. D 97, 116014 (2018a), arXiv:1712.03946 [hep-ph] .
  * Badger  _et al._ [2018b] S. Badger, C. Brønnum-Hansen, T. Gehrmann, H. B. Hartanto, J. Henn, N. A. Lo Presti, and T. Peraro, Applications of integrand reduction to two-loop five-point scattering amplitudes in QCD, PoS LL2018, 006 (2018b), arXiv:1807.09709 [hep-ph] .
  * Abreu  _et al._ [2018b] S. Abreu, F. Febres Cordero, H. Ita, B. Page, and V. Sotnikov, Planar Two-Loop Five-Parton Amplitudes from Numerical Unitarity, JHEP 11, 116, arXiv:1809.09067 [hep-ph] .
  * Hodges [2013] A. Hodges, Eliminating spurious poles from gauge-theoretic amplitudes, JHEP 05, 135, arXiv:0905.1473 [hep-th] .
  * Gehrmann  _et al._ [2018] T. Gehrmann, J. Henn, and N. Lo Presti, Pentagon functions for massless planar scattering amplitudes, JHEP 10, 103, arXiv:1807.09812 [hep-ph] .
  * Chicherin and Sotnikov [2020] D. Chicherin and V. Sotnikov, Pentagon Functions for Scattering of Five Massless Particles, JHEP 12, 167, arXiv:2009.07803 [hep-ph] .
  * Gehrmann  _et al._ [2016] T. Gehrmann, J. Henn, and N. Lo Presti, Analytic form of the two-loop planar five-gluon all-plus-helicity amplitude in QCD, Phys. Rev. Lett. 116, 062001 (2016), [Erratum: Phys.Rev.Lett. 116, 189903 (2016)], arXiv:1511.05409 [hep-ph] .
  * Badger  _et al._ [2019a] S. Badger, C. Brønnum-Hansen, H. B. Hartanto, and T. Peraro, Analytic helicity amplitudes for two-loop five-gluon scattering: the single-minus case, JHEP 01, 186, arXiv:1811.11699 [hep-ph] .
  * Abreu  _et al._ [2019a] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, and B. Page, Analytic Form of Planar Two-Loop Five-Gluon Scattering Amplitudes in QCD, Phys. Rev. Lett. 122, 082002 (2019a), arXiv:1812.04586 [hep-ph] .
  * Abreu  _et al._ [2019b] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page, and V. Sotnikov, Analytic Form of the Planar Two-Loop Five-Parton Scattering Amplitudes in QCD, JHEP 05, 084, arXiv:1904.00945 [hep-ph] .
  * Abreu  _et al._ [2020b] S. Abreu, B. Page, E. Pascual, and V. Sotnikov, Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders, (2020b), arXiv:2010.15834 [hep-ph] .
  * Chawdhry  _et al._ [2020a] H. A. Chawdhry, M. Czakon, A. Mitov, and R. Poncelet, Two-loop leading-color helicity amplitudes for three-photon production at the LHC, (2020a), arXiv:2012.13553 [hep-ph] .
  * Abreu  _et al._ [2019c] S. Abreu, L. J. Dixon, E. Herrmann, B. Page, and M. Zeng, The two-loop five-point amplitude in $\mathcal{N}=4$ super-Yang-Mills theory, Phys. Rev. Lett. 122, 121603 (2019c), arXiv:1812.08941 [hep-th] .
  * Chicherin  _et al._ [2019a] D. Chicherin, T. Gehrmann, J. Henn, P. Wasser, Y. Zhang, and S. Zoia, Analytic result for a two-loop five-particle amplitude, Phys. Rev. Lett. 122, 121602 (2019a), arXiv:1812.11057 [hep-th] .
  * Chicherin  _et al._ [2019b] D. Chicherin, T. Gehrmann, J. M. Henn, P. Wasser, Y. Zhang, and S. Zoia, The two-loop five-particle amplitude in $\mathcal{N}$ = 8 supergravity, JHEP 03, 115, arXiv:1901.05932 [hep-th] .
  * Abreu  _et al._ [2019d] S. Abreu, L. J. Dixon, E. Herrmann, B. Page, and M. Zeng, The two-loop five-point amplitude in $\mathcal{N}$ = 8 supergravity, JHEP 03, 123, arXiv:1901.08563 [hep-th] .
  * Badger  _et al._ [2019b] S. Badger, D. Chicherin, T. Gehrmann, G. Heinrich, J. Henn, T. Peraro, P. Wasser, Y. Zhang, and S. Zoia, Analytic form of the full two-loop five-gluon all-plus helicity amplitude, Phys. Rev. Lett. 123, 071601 (2019b), arXiv:1905.03733 [hep-ph] .
  * Agarwal  _et al._ [2021] B. Agarwal, F. Buccioni, A. von Manteuffel, and L. Tancredi, Two-loop leading colour QCD corrections to $q\bar{q}\to\gamma\gamma g$ and $qg\to\gamma\gamma q$, (2021), arXiv:2102.01820 [hep-ph] .
  * De Laurentis and Maître [2021] G. De Laurentis and D. Maître, Two-Loop Five-Parton Leading-Colour Finite Remainders in the Spinor-Helicity Formalism, JHEP 02, 016, arXiv:2010.14525 [hep-ph] .
  * Chawdhry  _et al._ [2020b] H. A. Chawdhry, M. L. Czakon, A. Mitov, and R. Poncelet, NNLO QCD corrections to three-photon production at the LHC, JHEP 02, 057, arXiv:1911.00479 [hep-ph] .
  * Kallweit  _et al._ [2020] S. Kallweit, V. Sotnikov, and M. Wiesemann, Triphoton production at hadron colliders in NNLO QCD 10.1016/j.physletb.2020.136013 (2020), arXiv:2010.04681 [hep-ph] .
  * Catani [1998] S. Catani, The Singular behavior of QCD amplitudes at two loop order, Phys. Lett. B427, 161 (1998), arXiv:hep-ph/9802439 [hep-ph] .
  * Becher and Neubert [2009a] T. Becher and M. Neubert, On the Structure of Infrared Singularities of Gauge-Theory Amplitudes, JHEP 06, 081, [Erratum: JHEP11,024(2013)], arXiv:0903.1126 [hep-ph] .
  * Becher and Neubert [2009b] T. Becher and M. Neubert, Infrared singularities of scattering amplitudes in perturbative QCD, Phys. Rev. Lett. 102, 162001 (2009b), [Erratum: Phys. Rev. Lett.111,no.19,199905(2013)], arXiv:0901.0722 [hep-ph] .
  * Gardi and Magnea [2009] E. Gardi and L. Magnea, Factorization constraints for soft anomalous dimensions in QCD scattering amplitudes, JHEP 03, 079, arXiv:0901.1091 [hep-ph] .
  * Nogueira [1993] P. Nogueira, Automatic Feynman graph generation, J. Comput. Phys. 105, 279 (1993).
  * Larin [1993] S. A. Larin, The Renormalization of the axial anomaly in dimensional regularization, Phys. Lett. B 303, 113 (1993), arXiv:hep-ph/9302240 .
  * Moch  _et al._ [2015] S. Moch, J. A. M. Vermaseren, and A. Vogt, On $\gamma$5 in higher-order QCD calculations and the NNLO evolution of the polarized valence distribution, Phys. Lett. B 748, 432 (2015), arXiv:1506.04517 [hep-ph] .
  * Kuipers  _et al._ [2013] J. Kuipers, T. Ueda, J. A. M. Vermaseren, and J. Vollinga, FORM version 4.0, Comput. Phys. Commun. 184, 1453 (2013), arXiv:1203.6543 [cs.SC] .
  * Ruijl  _et al._ [2017] B. Ruijl, T. Ueda, and J. Vermaseren, FORM version 4.2, (2017), arXiv:1707.06453 [hep-ph] .
  * Lee [2012] R. N. Lee, Presenting LiteRed: a tool for the Loop InTEgrals REDuction, (2012), arXiv:1212.2685 [hep-ph] .
  * Laporta [2000] S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations, Int. J. Mod. Phys. A15, 5087 (2000), arXiv:hep-ph/0102033 [hep-ph] .
  * Gluza  _et al._ [2011] J. Gluza, K. Kajda, and D. A. Kosower, Towards a Basis for Planar Two-Loop Integrals, Phys. Rev. D 83, 045012 (2011), arXiv:1009.0472 [hep-th] .
  * Ita [2016] H. Ita, Two-loop Integrand Decomposition into Master Integrals and Surface Terms, Phys. Rev. D 94, 116015 (2016), arXiv:1510.05626 [hep-th] .
  * Larsen and Zhang [2016] K. J. Larsen and Y. Zhang, Integration-by-parts reductions from unitarity cuts and algebraic geometry, Phys. Rev. D 93, 041701 (2016), arXiv:1511.01071 [hep-th] .
  * Kotikov [1991] A. V. Kotikov, Differential equations method: New technique for massive Feynman diagrams calculation, Phys. Lett. B 254, 158 (1991).
  * Bern  _et al._ [1994] Z. Bern, L. J. Dixon, and D. A. Kosower, Dimensionally regulated pentagon integrals, Nucl. Phys. B 412, 751 (1994), arXiv:hep-ph/9306240 .
  * Remiddi [1997] E. Remiddi, Differential equations for Feynman graph amplitudes, Nuovo Cim. A 110, 1435 (1997), arXiv:hep-th/9711188 .
  * Gehrmann and Remiddi [2000] T. Gehrmann and E. Remiddi, Differential equations for two loop four point functions, Nucl. Phys. B 580, 485 (2000), arXiv:hep-ph/9912329 .
  * Henn [2013] J. M. Henn, Multiloop integrals in dimensional regularization made simple, Phys. Rev. Lett. 110, 251601 (2013), arXiv:1304.1806 [hep-th] .
  * Moriello [2020] F. Moriello, Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops, JHEP 01, 150, arXiv:1907.13234 [hep-ph] .
  * Goncharov [1998] A. B. Goncharov, Multiple polylogarithms, cyclotomy and modular complexes, Math. Res. Lett. 5, 497 (1998), arXiv:1105.2076 [math.AG] .
  * Remiddi and Vermaseren [2000] E. Remiddi and J. A. M. Vermaseren, Harmonic polylogarithms, Int. J. Mod. Phys. A 15, 725 (2000), arXiv:hep-ph/9905237 .
  * Goncharov [2001] A. B. Goncharov, Multiple polylogarithms and mixed Tate motives, (2001), arXiv:math/0103059 .
  * Chen [1977] K.-T. Chen, Iterated path integrals, Bull. Am. Math. Soc. 83, 831 (1977).
  * Brown [2013] F. Brown, Iterated integrals in quantum field theory, in  _6th Summer School on Geometric and Topological Methods for Quantum Field Theory_ (2013) pp. 188–240.
  * Ferguson and Bailey [1992] H. R. P. Ferguson and D. H. Bailey, A Polynomial Time, Numerically Stable Integer Relation Algorithm, RNR Technical Report RNR-91-032 (1992).
  * Boehm  _et al._ [2020] J. Boehm, M. Wittmann, Z. Wu, Y. Xu, and Y. Zhang, IBP reduction coefficients made simple, JHEP 12, 054, arXiv:2008.13194 [hep-ph] .
  * Heller and von Manteuffel [2021] M. Heller and A. von Manteuffel, MultivariateApart: Generalized Partial Fractions, (2021), arXiv:2101.08283 [cs.SC] .
  * Note [1] Since estimates of evaluation time rely on system specific parameters, we have taken a conservative value for the speed improvement.
  * Caron-Huot  _et al._ [2020] S. Caron-Huot, D. Chicherin, J. Henn, Y. Zhang, and S. Zoia, Multi-Regge Limit of the Two-Loop Five-Point Amplitudes in $\mathcal{N}=4$ Super Yang-Mills and $\mathcal{N}=8$ Supergravity, JHEP 10, 188, arXiv:2003.03120 [hep-th] .
  * Chicherin  _et al._ [2020] D. Chicherin, J. M. Henn, and G. Papathanasiou, Cluster algebras for Feynman integrals, (2020), arXiv:2012.12285 [hep-th] .
  * Hidding [2020] M. Hidding, DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions, (2020), arXiv:2006.05510 [hep-ph] .