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

TL;DR

This paper provides the first numerical results for two-loop helicity amplitudes involving a W-boson and four partons in QCD, using innovative finite field sampling and sector decomposition techniques.

Contribution

It introduces a novel numerical approach to evaluate complex two-loop amplitudes with an incomplete master integral basis in QCD.

Findings

First numerical results for these amplitudes

Effective use of finite field sampling and sector decomposition

Identification of master integrals with simple divergence structure

Abstract

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.

Tables

Table 1. Table 1 : Master integral topologies made up of five external legs with local numerator insertions. The topology T 𝑇 T and numerator N ​ ( k i , p i , μ i ​ j ) 𝑁 subscript 𝑘 𝑖 subscript 𝑝 𝑖 subscript 𝜇 𝑖 𝑗 N(k_{i},p_{i},\mu_{ij}) correspond to the integral definition in Eq. ( 40 ). Numerator building blocks Ψ Ψ \Psi , Φ Φ \Phi and Ω Ω \Omega are defined in Eq. ( 4 ). δ 1 subscript 𝛿 1 \delta_{1} is the order at which the expansion of the master integral starts, while δ 2 subscript 𝛿 2 \delta_{2} the highest order in ϵ italic-ϵ \epsilon needed from the master integral for the amplitude evaluation.

$T$	$N(k_i,p_i,{\mu }_{ij})$	${\delta }_1$	${\delta }_2$	Crossing
	${\Omega }_{2;4|56}{\mu }_{11}$	$ϵ$	1	$\left(1\leftrightarrow 4,2\leftrightarrow 3\right)$
	${\Omega }_{2;4|56}^{\ast }{\mu }_{11}$	$ϵ$	1
	${\Psi }_{1;1\left|2\right|3}{\Omega }_{2;4|56}$	1	1
	${\Phi }_{2;2|3}{\mu }_{11}$	$ϵ$	1	—
	${\Phi }_{2;2|3}^{\ast }{\mu }_{11}$	$ϵ$	1
	${\Psi }_{1;4\left|56\right|1}{\Phi }_{2;2|3}$	1	1
	${\mu }_{12}$	$ϵ$	1	$\left(1\leftrightarrow 4,2\leftrightarrow 3\right)$
	${\Omega }_{1;4|56}{\Phi }_{2;2|3}$	1	1
	${\Omega }_{1;4|56}{\Phi }_{2;2|3}^{\ast }$	1	1
	${\mu }_{11}$	$ϵ$	1	$\left(1\leftrightarrow 4,2\leftrightarrow 3\right)$
	1	${ϵ}^{-2}$	$ϵ$	—
	${\Phi }_{1;2|3}$	1	$ϵ$
	${\Phi }_{1;2|3}^{\ast }$	1	$ϵ$
	${\Phi }_{1;2|3}{\Omega }_{2;1|56}$	1	$ϵ$
	${\Phi }_{1;2|3}{\Omega }_{2;1|56}^{\ast }$	1	$ϵ$
	${\Phi }_{1;2|3}^{\ast }{\Omega }_{2;1|56}$	1	$ϵ$
	${\Psi }_{1;2\left|3\right|4}$	${ϵ}^{-1}$	$ϵ$	$\left(1\leftrightarrow 4,2\leftrightarrow 3\right)$
	${\Psi }_{1;2\left|3\right|4}^{\ast }$	${ϵ}^{-1}$	$ϵ$
	${\Psi }_{1;4\left|56\right|1}$	${ϵ}^{-1}$	$ϵ$	—
	${\Psi }_{1;4\left|56\right|1}^{\ast }$	${ϵ}^{-1}$	$ϵ$

Table 2. Table 2 : Master integral topologies made up of four external legs with local numerator insertions. The topology T 𝑇 T and numerator N ​ ( k i , p i , μ i ​ j ) 𝑁 subscript 𝑘 𝑖 subscript 𝑝 𝑖 subscript 𝜇 𝑖 𝑗 N(k_{i},p_{i},\mu_{ij}) correspond to the integral definition in Eq. ( 40 ). Numerator building blocks Ψ Ψ \Psi , Φ Φ \Phi and Ω Ω \Omega are defined in Eq. ( 4 ). δ 1 subscript 𝛿 1 \delta_{1} is the order at which the expansion of the master integral starts, while δ 2 subscript 𝛿 2 \delta_{2} the highest order in ϵ italic-ϵ \epsilon needed from the master integral for the amplitude evaluation.

$T$	$N(k_i,p_i,{\mu }_{ij})$	${\delta }_1$	${\delta }_2$	Crossing
	${\Omega }_{1;1|23}{\Omega }_{2;4|56}$	1	1	$\left(1\leftrightarrow 4,2\leftrightarrow 3\right)$
	${\Omega }_{1;1|23}{\Omega }_{2;4|56}^{\ast }$	1	1
	${\Omega }_{1;4|56}{\Omega }_{2;3|12}$	1	1	$\left(1\leftrightarrow 4,2\leftrightarrow 3\right)$
	${\Omega }_{1;4|56}{\Omega }_{2;3|12}^{\ast }$	1	1
	${\Omega }_{1;3|12}$	1	1	$\left(1\leftrightarrow 4,2\leftrightarrow 3\right)$
	${\Omega }_{1;1|23}$	1	1	$\left(1\leftrightarrow 4\right)$ $\left(2\leftrightarrow 5,3\leftrightarrow 6\right)$ $\left(1\leftrightarrow 4,2\leftrightarrow 5,3\leftrightarrow 6\right)$

Table 3. Table 3 : The numerical results of A ^ ( 2 ) , [ i ] ​ ( 1 q , 2 g , 3 g , 4 q ¯ ′ , 5 ν ¯ , 6 ℓ ) superscript ^ 𝐴 2 delimited-[] 𝑖 subscript 1 𝑞 subscript 2 𝑔 subscript 3 𝑔 subscript 4 superscript ¯ 𝑞 ′ subscript 5 ¯ 𝜈 subscript 6 ℓ \hat{A}^{(2),[i]}\left(1_{q},2_{g},3_{g},4_{\bar{q}^{\prime}},5_{\bar{\nu}},6_{\ell}\right) using kinematic point in Eq. ( 44 ) for each ( d s − 2 ) subscript 𝑑 𝑠 2 (d_{s}-2) component.

$qgg{\overline{q}}^{\prime }\overline{\nu }\ell$	${ϵ}^{-4}$	${ϵ}^{-3}$	${ϵ}^{-2}$	${ϵ}^{-1}$	${ϵ}^0$
${\hat{A}}_{-++++-}^{\left(2\right),\left[0\right]}$	4.50000	-2.38581(9)	17.0943(2)	69.09(3)	-137.1(3)
${\hat{A}}_{-++++-}^{\left(2\right),\left[1\right]}$	0.00000	-0.62498(3)	-147.7288(3)	-347.221(9)	863.63(8)
${\hat{A}}_{-++++-}^{\left(2\right),\left[2\right]}$	0.00000	0.00000	-0.031254(7)	-3.72543(2)	-65.7967(5)
$qgg{\overline{q}}^{\prime }\overline{\nu }\ell$	${ϵ}^{-4}$	${ϵ}^{-3}$	${ϵ}^{-2}$	${ϵ}^{-1}$	${ϵ}^0$
${\hat{A}}_{-+-++-}^{\left(2\right),\left[0\right]}$	4.50000	-2.3858(1)	-16.3282(4)	-1.397(3)	70.05(4)
${\hat{A}}_{-+-++-}^{\left(2\right),\left[1\right]}$	0.00000	-0.624996(6)	0.7603(2)	4.457(1)	1.04(1)
${\hat{A}}_{-+-++-}^{\left(2\right),\left[2\right]}$	0.00000	0.00000	-0.0312498(2)	-0.025800(5)	0.72620(5)
$qgg{\overline{q}}^{\prime }\overline{\nu }\ell$	${ϵ}^{-4}$	${ϵ}^{-3}$	${ϵ}^{-2}$	${ϵ}^{-1}$	${ϵ}^0$
${\hat{A}}_{--+++-}^{\left(2\right),\left[0\right]}$	4.50000	-2.38579(5)	-22.18117(7)	-16.113(6)	90.06(4)
${\hat{A}}_{--+++-}^{\left(2\right),\left[1\right]}$	0.00000	-0.625000(1)	1.131987(9)	5.7364(2)	-2.1289(6)
${\hat{A}}_{--+++-}^{\left(2\right),\left[2\right]}$	0.00000	0.00000	-0.0312502(3)	0.005162(2)	1.21279(4)

Table 4. Table 4 : The numerical results of A ^ ( 2 ) , [ i ] ​ ( 1 q , 2 Q ¯ , 3 Q , 4 q ¯ ′ , 5 ν ¯ , 6 ℓ ) superscript ^ 𝐴 2 delimited-[] 𝑖 subscript 1 𝑞 subscript 2 ¯ 𝑄 subscript 3 𝑄 subscript 4 superscript ¯ 𝑞 ′ subscript 5 ¯ 𝜈 subscript 6 ℓ \hat{A}^{(2),[i]}\left(1_{q},2_{\bar{Q}},3_{Q},4_{\bar{q}^{\prime}},5_{\bar{\nu}},6_{\ell}\right) using kinematic point in Eq. ( 44 ) for each ( d s − 2 ) subscript 𝑑 𝑠 2 (d_{s}-2) component.

$q\overline{Q}Q{\overline{q}}^{\prime }\overline{\nu }\ell$	${ϵ}^{-4}$	${ϵ}^{-3}$	${ϵ}^{-2}$	${ϵ}^{-1}$	${ϵ}^0$
${\hat{A}}_{-+-++-}^{\left(2\right),\left[0\right]}$	2.00000	-6.00283(9)	-12.7724(2)	31.869(6)	158.89(6)
${\hat{A}}_{-+-++-}^{\left(2\right),\left[1\right]}$	0.00000	-0.583333(1)	0.96122(6)	5.2453(4)	2.853(3)
${\hat{A}}_{-+-++-}^{\left(2\right),\left[2\right]}$	0.00000	0.00000	-0.055555	-0.240170(1)	-0.25365(2)
$q\overline{Q}Q{\overline{q}}^{\prime }\overline{\nu }\ell$	${ϵ}^{-4}$	${ϵ}^{-3}$	${ϵ}^{-2}$	${ϵ}^{-1}$	${ϵ}^0$
${\hat{A}}_{--+++-}^{\left(2\right),\left[0\right]}$	2.00000	-6.00282(8)	-18.0013(1)	34.592(7)	222.52(6)
${\hat{A}}_{--+++-}^{\left(2\right),\left[1\right]}$	0.00000	-0.583334(1)	2.059832(8)	5.4211(2)	-13.5049(5)
${\hat{A}}_{--+++-}^{\left(2\right),\left[2\right]}$	0.00000	0.00000	-0.055555	-0.081689(2)	1.10832(2)

Table 5. Table 5 : The numerical comparison of A ^ ( 2 ) ​ ( 1 q , 2 g , 3 g , 4 q ¯ ′ , 5 ν ¯ , 6 ℓ ) superscript ^ 𝐴 2 subscript 1 𝑞 subscript 2 𝑔 subscript 3 𝑔 subscript 4 superscript ¯ 𝑞 ′ subscript 5 ¯ 𝜈 subscript 6 ℓ \hat{A}^{(2)}\left(1_{q},2_{g},3_{g},4_{\bar{q}^{\prime}},5_{\bar{\nu}},6_{\ell}\right) with the universal pole structure 𝒫 ( 2 ) superscript 𝒫 2 {\mathcal{P}}^{(2)} defined in Eq. ( 24 ), using kinematic point of Eq. ( 44 ), in the HV scheme.

$qgg{\overline{q}}^{\prime }\overline{\nu }\ell$	${ϵ}^{-4}$	${ϵ}^{-3}$	${ϵ}^{-2}$	${ϵ}^{-1}$	${ϵ}^0$
${\hat{A}}_{-++++-}^{\left(2\right)}$	4.50000	-3.63577(3)	-277.2182(7)	-344.56(1)	2051.1(2)
$P_{-++++-}^{\left(2\right)}$	4.5	-3.63576	-277.2186	-344.569(6)	—
${\hat{A}}_{-+-++-}^{\left(2\right)}$	4.50000	-3.63581(9)	-13.6826(2)	6.143(5)	66.21(7)
$P_{-+-++-}^{\left(2\right)}$	4.5	-3.63576	-13.6824	6.145(1)	—
${\hat{A}}_{--+++-}^{\left(2\right)}$	4.50000	-3.63579(5)	-18.79219(7)	-6.633(6)	79.02(4)
$P_{--+++-}^{\left(2\right)}$	4.5	-3.63576	-18.79212	-6.6303(5)	—

Table 6. Table 6 : The numerical comparison of A ^ ( 2 ) ​ ( 1 q , 2 Q ¯ , 3 Q , 4 q ¯ ′ , 5 ν ¯ , 6 ℓ ) superscript ^ 𝐴 2 subscript 1 𝑞 subscript 2 ¯ 𝑄 subscript 3 𝑄 subscript 4 superscript ¯ 𝑞 ′ subscript 5 ¯ 𝜈 subscript 6 ℓ \hat{A}^{(2)}\left(1_{q},2_{\bar{Q}},3_{Q},4_{\bar{q}^{\prime}},5_{\bar{\nu}},6_{\ell}\right) with the universal pole structure 𝒫 ( 2 ) superscript 𝒫 2 {\mathcal{P}}^{(2)} defined in Eq. ( 24 ), using kinematic point of Eq. ( 44 ), in the HV scheme.

$q\overline{Q}Q{\overline{q}}^{\prime }\overline{\nu }\ell$	${ϵ}^{-4}$	${ϵ}^{-3}$	${ϵ}^{-2}$	${ϵ}^{-1}$	${ϵ}^0$
${\hat{A}}_{-+-++-}^{\left(2\right)}$	2.00000	-7.16949(9)	-9.9055(2)	39.922(6)	154.79(7)
$P_{--+++-}^{\left(2\right)}$	2	-7.16944	-9.9054	39.9245(8)	—
${\hat{A}}_{--+++-}^{\left(2\right)}$	2.00000	-7.16948(8)	-12.9371(1)	41.432(8)	189.53(6)
$P_{--+++-}^{\left(2\right)}$	2	-7.16944	-12.9370	41.4353(6)	—

Equations

$${\mathcal{A}}_n^{(2)}(\{p\})=\sum _c{\mathcal{C}}_c{A}_{n,c}^{(2)}(\{p\}),$$

$$A_n^{(2)}(\{p\})=\int \prod _{i=1}^2\frac{d^d{k}_i}{i{\pi }^{d/2}{e}^{-ϵ{\gamma }_E}}\sum _T\frac{N_T(d_s,\{k\},\{p\})}{{\prod }_{\alpha \in T}{D}_{\alpha }(\{k\},\{p\})},$$

$$N_{T}(d_{s},\{k\},\{p\})=N_{T}\big{(}d_{s},k_{i}.k_{j},\mu_{ij},k_{i}.q_{j},q_{i}.q_{j},\bar{u}(p_{i})f(k,q)u(p_{j})\big{)},$$

$$k_i={\bar{k}}_i+{\tilde{k}}_i.$$

$$A_n^{(2),h}(\{p\})=\int \prod _{i=1}^2\frac{d^d{k}_i}{i{\pi }^{d/2}{e}^{-ϵ{\gamma }_E}}\sum _T\frac{N_T^h(d_s,\{k\},\{p\})}{{\prod }_{\alpha \in T}{D}_{\alpha }(\{k\},\{p\})},$$

$$N^{h}_{T}(d_{s},\{k\},\{p\})=N^{h}_{T}\big{(}d_{s},x_{i},k_{i}.k_{j},\bar{k}_{i}.p_{j},\mu_{ij},\langle p_{a}|\bar{k}_{i}|p_{b}],\langle p_{a}|\bar{k}_{i}|\bar{k}_{j}|p_{b}\rangle,[p_{a}|\bar{k}_{i}|\bar{k}_{j}|p_{b}]\big{)}.$$

$$A_n^{(2),h}(\{p\})=\int \prod _{i=1}^2\frac{d^d{k}_i}{i{\pi }^{d/2}{e}^{-ϵ{\gamma }_E}}\sum _T\frac{{\Delta }_T^h(d_s,\{k\},\{p\})}{{\prod }_{\alpha \in T}{D}_{\alpha }(\{k\},\{p\})},$$

$$G_{a_1{a}_2{a}_3{a}_4{a}_5{a}_6{a}_7{a}_8{a}_9{a}_{10}{a}_{11}}$$

$$\times \frac{1}{k_1^{2a_1}}\frac{1}{(k_1-p_1{)}^{2a_2}}\frac{1}{(k_1-p_1-p_2{)}^{2a_3}}\frac{1}{(k_1+p_4+p_{56}{)}^{2a_4}}$$

$$\times \frac{1}{k_2^{2a_5}}\frac{1}{(k_2-p_{56}{)}^{2a_6}}\frac{1}{(k_2-p_4-p_{56}{)}^{2a_7}}\frac{1}{(k_1+k_2{)}^{2a_8}}$$

$$\times \frac{1}{(k_1+p_{56}{)}^{2a_9}}\frac{1}{(k_2+p_1{)}^{2a_{10}}}\frac{1}{(k_2+p_1+p_2{)}^{2a_{11}}},$$

$$\displaystyle\Delta^{h}\bigg{(}\parbox{45.52458pt}{\includegraphics[width=45.52458pt,trim=0.0pt 0.0pt 25.09373pt 0.0pt]{figs/twomassdbox}}\bigg{)}=\sum c^{h}_{(1,a_{2},1,1,1,1,1,1,a_{9},a_{10},a_{11})}$$

$$\times$$

$$-4$$

$$-6$$

$${\bar{k}}_i^{\mu }$$

$$⟨p_a|{\bar{k}}_i|p_b]$$

$${\mu }_{ij}$$

$$\displaystyle\Delta^{h}\bigg{(}\parbox{45.52458pt}{\includegraphics[width=45.52458pt,trim=0.0pt 0.0pt 25.09373pt 0.0pt]{figs/twomassdbox}}\bigg{)}+\frac{\Delta^{h}\bigg{(}\parbox{42.67912pt}{\includegraphics[width=42.67912pt,trim=0.0pt 0.0pt 40.15pt 0.0pt]{figs/T431_1}}\bigg{)}}{(k_{1}-p_{1})^{2}}=N^{h}\bigg{(}\parbox{45.52458pt}{\includegraphics[width=45.52458pt,trim=0.0pt 0.0pt 25.09373pt 0.0pt]{figs/twomassdbox}}\bigg{)}+\frac{N^{h}\bigg{(}\parbox{42.67912pt}{\includegraphics[width=42.67912pt,trim=0.0pt 0.0pt 40.15pt 0.0pt]{figs/T431_1}}\bigg{)}}{(k_{1}-p_{1})^{2}},$$

$$A_n^{(2),h}(\{p\})=\sum _{\mathbf{a}}{c}_{\mathbf{a}}^h(x_i)G_{\mathbf{a}},$$

$$A_n^{(2),h}(\{p\})=\sum _k{c}_k^{\mathrm{IBP},h}(x_i,ϵ)J_k(\{p\},ϵ).$$

$${\mathcal{A}}^{(L)}(1_q,2_{\bar{Q}},3_Q,4_{{\bar{q}}^{\prime }},5_{\bar{\nu }},6_{\ell })$$

$${\mathcal{A}}^{(L)}(1_q,2_g,3_g,4_{{\bar{q}}^{\prime }},5_{\bar{\nu }},6_{\ell })$$

$$Z=\left(\begin{array}{cccccc}1& 0& y_1& y_2& y_3& y_4\\ 0& 1& 1& 1& 1& 1\\ 0& 0& 0& \frac{x_5}{x_2}& x_6& 1\\ 0& 0& 1& 1& x_7& 1-\frac{x_8}{x_5}\end{array}\right),$$

$$x_1$$

$$x_5$$

$${\mathcal{P}}^{(1)}$$

$${\mathcal{P}}^{(2)}$$

$$I_2(ϵ)=-\frac{1}{2}{I}_1(ϵ)\left[I_1(ϵ)+\frac{{\beta }_0}{ϵ}\right]+\frac{N(ϵ)}{N(2ϵ)}\left[\frac{{\beta }_0}{2ϵ}+\frac{{\gamma }_1^{\mathrm{cusp}}}{8}\right]I_1(2ϵ)+H^{(2)}(ϵ),$$

$$I_1^{q\bar{Q}Q{\bar{q}}^{\prime }\bar{\nu }\ell }(ϵ)$$

$$I_1^{qgg{\bar{q}}^{\prime }\bar{\nu }\ell }(ϵ)$$

$$H_{q\bar{Q}Q{\bar{q}}^{\prime }\bar{\nu }\ell }^{(2)}(ϵ)$$

$$H_{qgg{\bar{q}}^{\prime }\bar{\nu }\ell }^{(2)}(ϵ)$$

$${\beta }_0=$$

$${\beta }_1=$$

$${\gamma }_0^g=$$

Full text

††preprint: IPPP/19/56, ZU-TH 33/19

1 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\to 2$ scattering and N3LO for $2\to 1$ processes111For recent summaries of the state-of-the-art see Badger:2016bpw ; Bendavid:2018nar ..

The need to match experimental precision has led to increasing efforts from the theoretical community to develop new techniques for $2\to 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 vonManteuffel:2014ixa ; Peraro:2016wsq ; Peraro:2019svx have recently produced the first analytic results for five-parton amplitudes in the leading colour approximation. Combined with the recently computed analytic master integrals Papadopoulos:2015jft ; Gehrmann:2018yef using the canonical basis approach to differential equations, a form suitable for combination with the unresolved contributions to the cross section has been obtained Gehrmann:2015bfy ; Badger:2018enw ; Abreu:2018zmy ; Abreu:2019odu .

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\to W+j$ was among the first $2\to 2$ process computed at NLO Giele:1993dj . The amplitudes for $pp\to W+2j$ were computed using the recently developed on-shell unitarity method Bern:1994zx ; Bern:1994cg ; Bern:1996ka ; Bern:1997sc and were implemented into MCFM to provide differential cross-section predictions Campbell:2002tg . NLO results for associated $W$-boson production with three or more jets are accessible through automation and the use of the generalised unitarity method Ellis:2008qc ; Ellis:2009zw ; Berger:2008sj ; Berger:2009ep ; Berger:2010zx ; Bern:2013gka .

NNLO corrections to $pp\to 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) Tkachov:1981wb ; Chetyrkin:1981qh reduction, which in recent times has involved following Laporta’s algorithm Laporta:2001dd , 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 Kotikov:1990kg ; Gehrmann:1999as ; Henn:2013pwa has been employed to find analytic expressions, most recently for the complete set of planar Gehrmann:2015bfy ; Gehrmann:2018yef and non-planar integrals Abreu:2018rcw ; Boehm:2018fpv ; Abreu:2018aqd ; Chicherin:2018mue ; Chicherin:2018old . For the case, in which the amplitudes considered here fall, only one of the three planar families has been evaluated Papadopoulos:2015jft . 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 Gluza:2010ws ; vonManteuffel:2014ixa ; Larsen:2015ped ; Ita:2015tya ; Kosower:2018obg and public implementations Lee:2012cn ; Smirnov:2014hma ; Smirnov:2019qkx ; vonManteuffel:2012np ; Maierhoefer:2017hyi . Applications to five-particle problems have been possible though yielded large IBP reduction tables Boels:2018nrr ; Chawdhry:2018awn . 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. Badger:2018enw ; Peraro:2019svx , 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 Ossola:2006us 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 Mastrolia:2011pr ; Mastrolia:2012an ; Zhang:2012ce ; Badger:2012dp ; Mastrolia:2012wf ; Mastrolia:2013kca ; Mastrolia:2016dhn and methods to employ unitarity cuts Bern:1994zx ; Bern:1994cg to build amplitudes by directly incorporating IBP decomposition have been established Kosower:2011ty ; CaronHuot:2012ab ; Ita:2015tya ; Abreu:2017xsl ; Abreu:2017idw .

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 Badger:2017jhb ; Abreu:2017hqn ; Badger:2018gip ; Abreu:2018jgq . These algorithms have been generalised to allow for a full reconstruction of the coefficients of the pentagon functions classified in Gehrmann:2018yef leading to an analytic form of the single-minus helicity amplitudes Badger:2018enw and the complete leading colour five-parton helicity amplitudes within the numerical unitarity framework Abreu:2018zmy ; Abreu:2019odu . The success of computations in the planar sector has shifted focus to the non-planar sector of massless two-loop five-point amplitudes with a series of new results in super-symmetric Yang-Mills Abreu:2018aqd ; Chicherin:2018yne and gravity Abreu:2019rpt ; Chicherin:2019xeg as well as in the all-plus sector of QCD Badger:2019djh .

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 Peraro:2019svx , 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 Binoth:2000ps ; Smirnov:2015mct ; Borowka:2017idc . Analytic results for the following classes of master integral are available: one of the three families of the off-shell five-point pentagon-box Papadopoulos:2015jft and four-point functions with one Gehrmann:2000zt and two off-shell Gehrmann:2015ora ; vonManteuffel:2015msa ; Henn:2014lfa ; Caola:2014lpa 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 ArkaniHamed:2010kv ; ArkaniHamed:2010gh with simplified divergence structure and therefore easier numerical evaluation. We consider both the $q\bar{Q}Q\bar{q}^{\prime}\bar{\nu}\ell$ and $qgg\bar{q}^{\prime}\bar{\nu}\ell$ 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.

2 Calculational framework

The framework described in this section is a modification of the numerical algorithm for two-loop amplitudes presented in Badger:2018enw 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 Nogueira:1991ex 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

[TABLE]

where $\mathcal{A}^{(2)}_{n}(\{p\})$ is the two-loop colour-dressed $n$-point amplitude, $A^{(2)}_{n}(\{p\})$ is the two-loop colour-stripped $n$-point amplitude and $\mathcal{C}_{c}$ is the corresponding colour factor. The colour-stripped amplitude is made up of numerator functions, $N_{T}(\{k\},\{p\})$, and a set of loop propagator denominators, $D_{\alpha}(\{k\},\{p\})$, for each diagram topology $T$

[TABLE]

where $k_{i}$ is the loop momenta, $d=4-2\epsilon$ is the space-time dimension and $d_{s}=g^{\mu}_{\mu}$ is the spin dimension. The loop amplitude in t’Hooft-Veltman (HV) scheme tHooft:1972tcz can be obtained by setting $d_{s}=d$, while the Four-Dimensional-Helicity (FDH) scheme Bern:2002zk 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 Kuipers:2012rf ; Ruijl:2017dtg and the Spinney library Cullen:2010jv . In general, the explicit functional dependence of the numerator function at this point is given by

[TABLE]

where $q_{i}=\{p_{i},\varepsilon_{i}\}$ ($\varepsilon_{i}$ is the polarisation vector of the external vector boson) and $\bar{u}(p_{i})f(k,q)u(p_{j})$ is a spinor string made up of slashed momenta ($\not{q}_{i}$ and $\not{k}_{i}$). The $d$-dimensional loop momenta can be decomposed into a four-dimensional part and an extra-dimensional part

[TABLE]

Due to rotational invariance in the extra dimensions, $\tilde{k}_{i}$ can only appear in the numerator function as $\mu_{ij}=-\tilde{k}_{i}\cdot\tilde{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}$ Hodges:2009hk . This allows us to express the spinor products of external momenta $(\langle ij\rangle,[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

[TABLE]

where the explicit functional dependence on the helicity-dependent numerator function $N^{h}_{T}(\{k\},\{p\})$ is

[TABLE]

In processing the algebraic expressions in Eqs. (1), (2) and (5), we have used in-house Form Kuipers:2012rf ; Ruijl:2017dtg and Mathematica scripts.

In this form the helicity-dependent numerator functions in Eq. (5) must be re-expressed in terms of integral families that can later be reduced using IBP equations. To achieve this we apply an integrand reduction algorithm to obtain

[TABLE]

where $\Delta$ is the irreducible numerator for an independent topology $T$. In order to determine $\Delta$, 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

[TABLE]

where $p_{ij\cdots k}=p_{i}+p_{j}+\cdots+p_{k}$. Up to cyclic permutations of the external legs, all integrals appearing in Eq. (7) can be written in the form of Eq. (8). We follow the conventions used in Badger:2018enw where negative exponents, $a_{i}<0$, correspond to the ISPs of the irreducible numerator, $\Delta_{T}^{h}$. 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

[TABLE]

where the figure represents the topology $T$. The bounds on the exponents are

[TABLE]

The helicity-dependent coefficients are functions of the spin dimension, $d_{s}$, and the external kinematics, $c^{h}=c^{h}(d_{s},\{p\})$. 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 momenta

[TABLE]

The coefficients of the spanning vectors, $p_{j}$, are functions of the inverse propagators and ISPs, $a_{ij}=a_{ij}(D_{\alpha},\text{ISPs})$. $a_{ij}$ can be determined by solving a linear system of equations constructed by contracting Eq. (13) with the spanning vectors, $p_{j}$. All variables in the numerators Eq. (6) can then be expressed in terms of these coefficients. For example

[TABLE]

The last relation is obtained by squaring Eq. (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 Baikov:1996iu , 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

[TABLE]

which is valid only when both sides are evaluated on the hepta-cut for the two-mass double box, i.e. $D_{\alpha}=0$ for $\alpha=1,3,\dots,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. (8),

[TABLE]

where we sum over tuples $\mathbf{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_{\mathbf{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. Bourjaily:2017wjl , 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 Laporta:2001dd with the aid of LiteRed Lee:2012cn , and solved numerically over finite fields within the FiniteFlow framework Peraro:2019svx . We can finally write the helicity amplitudes in the master integral basis $J_{k}$

[TABLE]

3 Planar two-loop $W$ plus four parton scattering

The number of Feynman diagrams contributing to $q\bar{Q}Q\bar{q}^{\prime}\bar{\nu}\ell$ and $qgg\bar{q}^{\prime}\bar{\nu}\ell$ processes at leading colour are 210 and 603, respectively, and the leading colour partial amplitudes are extracted according to,

[TABLE]

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 note that the vector boson only couples to the quark line connecting $q$ and $\bar{q}^{\prime}$ and does not couple to the equal flavour quark pair $Q,\bar{Q}$.

We choose a rational parametrisation of the massless $2\to 4$ kinematics using the momentum twistor parametrisation Hodges:2009hk

[TABLE]

where $y_{i}=\sum_{j=1}^{i}\prod_{k=1}^{j}\frac{1}{x_{k}}$ and

[TABLE]

We stress that while we generate a complete parametrisation for the $2\to 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. (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, $\epsilon$. The $\epsilon$-expanded partial amplitude contains a divergent part, manifested by the poles in $\epsilon$, 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 Catani:1998bh ; Becher:2009qa ; Becher:2009cu ; Gardi:2009qi .

The pole structure of the unrenormalised amplitude in the HV scheme at one and two loops is given by

[TABLE]

where $\hat{A}^{(1)}$ is the unrenormalised one-loop amplitude normalised to the tree-level amplitude. The $I_{2}(\epsilon)$ operator is defined by

[TABLE]

while the $I_{1}(\epsilon)$ operators for the $W+4$ parton process at leading colour are

[TABLE]

where $N(\epsilon)={e^{\epsilon\gamma_{E}}}/{\Gamma(1-\epsilon)}$ and

[TABLE]

Note that the $H_{q\bar{Q}Q\bar{q}^{\prime}\bar{\nu}\ell}^{(2)}(\epsilon)$ and $H_{qgg\bar{q}^{\prime}\bar{\nu}\ell}^{(2)}(\epsilon)$ functions are given in the leading colour limit. The $\beta$ function coefficients and anomalous dimensions without the contribution from closed fermion loops $N_{f}$ are

[TABLE]

where $C_{A}=N_{c}$, $C_{F}=(N_{c}^{2}-1)/(2N_{c})$.

With the pole structures available, we can check that the divergent part of the two-loop 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. (24) is expanded to $\mathcal{O}(\epsilon^{0})$. The analytic form of the finite remainder is in general much simpler than the finite part, as demonstrated in Badger:2018enw ; Abreu:2018zmy ; Abreu:2019odu ; Badger:2019djh .

4 Local master integrals

The master integrals appearing in the two-loop leading colour $q\bar{Q}Q\bar{q}^{\prime}\bar{\nu}\ell$ and $qgg\bar{q}^{\prime}\bar{\nu}\ell$ amplitudes, that are not known analytically, are evaluated numerically using sector decomposition Binoth:2000ps ; Smirnov:2015mct ; Borowka:2017idc . 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 ArkaniHamed:2010kv ; ArkaniHamed:2010gh 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 Panzer:2014gra ; vonManteuffel:2014qoa . 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

[TABLE]

or

[TABLE]

labelled according to Eq. (8). In Eq. (38), we use a notation for the integral with a numerator insertion $N(k_{i},p_{i},\mu_{ij})$

[TABLE]

where $T$ is the diagram topology that we will specify by drawing it, and $D_{\alpha}$ is a set of massless propagator denominator for a given topology $T$. The Laurent expansion for those master integrals starts at $\mathcal{O}(\epsilon^{-4})$ for $G_{11111111000}$ and $\mathcal{O}(\epsilon^{-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 Badger:2016ozq 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

[TABLE]

where $\operatorname{\rm tr}_{\pm}(ijkl)=\frac{1}{2}\operatorname{\rm tr}((1\pm\gamma_{5})\not{p}_{i}\not{p}_{j}\not{p}_{k}\not{p}_{l})$. The first two integrals in Eq. (41) evaluate to $\mathcal{O}(\epsilon)$ 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 $\epsilon$, therefore, the master integrals need to be expanded to $\mathcal{O}(1)$. To evaluate the last integral in Eq. (41) using sector decomposition method, we first need to write the numerator insertion in terms of scalar products and momentum twistor variables ($k_{i}\cdot k_{j}$, $k_{i}\cdot p_{j}$, $x_{i}$) using loop momentum decomposition given in Eqs. (13) - (15)

[TABLE]

The integral on the RHS can be directly evaluated using pySecDec Borowka:2017idc 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

[TABLE]

We also include in Tables 1 and 2 the order at which the expansion of the master integral starts, $\mathcal{O}(\delta_{1})$, the highest order in $\epsilon$ needed from the master integral for the amplitude evaluation, $\mathcal{O}(\delta_{2})$, as well as the list of possible permutations/crossings to obtain the full list of master integral topologies with local numerator insertions. $\delta_{1}>\delta_{2}$ indicates that an integral does not contribute to the finite part of the amplitude, while $\delta_{1}=\delta_{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 $\epsilon$ expansion of the integral is required.

5 Numerical results

We select a random Euclidean phase-space point by choosing rational values in the momentum twistor parametrisation from Eq. (21),

[TABLE]

The numerical results of the leading colour partial amplitudes are obtained by evaluating the master integrals in three different ways:

We make use of the master integral solutions that are known analytically and readily available for evaluation in the Euclidean kinematics Gehrmann:2000zt ; Papadopoulos:2015jft . 2. 2.

Master integrals defined with the local numerator insertions, shown in Tables 1 and 2, are numerically evaluated using pySecDec Borowka:2017idc . 3. 3.

The remaining master integrals, which contain either one-loop squared topologies or 5 propagator or fewer two-loop topologies, are evaluated directly using Fiesta Smirnov:2015mct and pySecDec Borowka:2017idc .222The one-loop integrals, required to $\mathcal{O}(\epsilon^{2})$, are evaluated numerically. Since these terms do not cause any stability issues there was no need to consider optimisation with analytic expressions.

Note that there is a class of four-point master integrals with two off shell legs, that are not covered in Papadopoulos:2015jft , but available in Henn:2014lfa ; Gehrmann:2015ora ; vonManteuffel:2015msa , 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 numerically333 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\to 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.. 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 $q\bar{Q}Q\bar{q}^{\prime}\bar{\nu}\ell$ and $qgg\bar{q}^{\prime}\bar{\nu}\ell$ helicity amplitudes, normalised to the tree level amplitude

[TABLE]

with helicities $\lambda_{i}$. We can further split the amplitude into components of $d_{s}-2$

[TABLE]

In Tables 3 and 4 we display numerical evaluations of the helicity amplitudes for $qgg\bar{q}^{\prime}\bar{\nu}\ell$ and $q\bar{Q}Q\bar{q}^{\prime}\bar{\nu}\ell$ channels, respectively, using the kinematic point given in Eq. (44) 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. (24) in the HV scheme. To obtain the two-loop pole structure of Eq. (24) up to the single pole, we need to compute the one-loop amplitude up to $\mathcal{O}(\epsilon)$ for both the $qgg\bar{q}^{\prime}\bar{\nu}\ell$ and $q\bar{Q}Q\bar{q}^{\prime}\bar{\nu}\ell$ processes. The one-loop 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 $\mathcal{O}(\epsilon)$ part of the two-mass easy box and one-mass pentagon integrals are evaluated numerically using $\textsc{Fiesta}/\textsc{pySecDec}$, while the rest are obtained from available analytic expressions Ellis:2007qk ; Papadopoulos:2014lla . Therefore, the numerical values quoted for the poles in Tables 5 and 6 are exact up to $\mathcal{O}(\epsilon^{-2})$. We assess the uncertainty of the $\mathcal{O}(\epsilon^{-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 $\mu_{ij}$. The integrals containing $\mu_{ij}$ are first written as dimension-shifted integrals and all integrals appearing in the amplitude are evaluated directly using $\textsc{Fiesta}/\textsc{pySecDec}$. This approach is similar to the method we employed for the numerical evaluation of the planar two-loop five gluon amplitude Badger:2017jhb . The results obtained using this method are in perfect agreement with the results reported in this paper.

6 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 Chicherin:2017dob ; Abreu:2018aqd ; Abreu:2019rpt ; Chicherin:2018yne ; Chicherin:2018mue ; Chicherin:2019xeg ; Chicherin:2018old , 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 Badger:2018enw ; Abreu:2018zmy ; Abreu:2019odu .

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.

Acknowledgements.

We thank Johannes Henn, Thomas Gehrmann, Andreas von Manteuffel for helpful discussions. We also thank Johannes Schlenk for his assistance with pySecDec. SB is supported by an STFC Rutherford Fellowship ST/L004925/1 and CBH and HBH received partial support from Rutherford Grant ST/M004104/1. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 772099. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement 746223.

Appendix A Kinematic invariants at the benchmark phase space point

From the momentum twistor parametrisation in Eq. (44) it is possible to evaluate all kinematic quantities in the external momenta. Since it may be useful to have reference values in momentum space we list here six two-particle and three three-particle invariants,

[TABLE]

In four dimensions only eight of these are independent due to the vanishing of the Gram determinant of the five independent momenta.