Analytic Form of the Planar Two-Loop Five-Parton Scattering Amplitudes in QCD

TL;DR

This paper derives the complete analytic expressions for all leading-color two-loop five-parton helicity amplitudes in QCD, facilitating precise calculations of three-jet production at hadron colliders.

Contribution

It introduces a novel method combining finite field evaluations and variable choices to efficiently reconstruct and simplify two-loop five-parton amplitudes analytically.

Findings

Analytic expressions for all leading-color two-loop five-parton amplitudes obtained.

Amplitudes are simplified to a compact form using multivariate partial fractions.

Results enable next-to-next-to-leading order QCD corrections for three-jet production.

Abstract

We present the analytic form of all leading-color two-loop five-parton helicity amplitudes in QCD. The results are analytically reconstructed from exact numerical evaluations over finite fields. Combining a judicious choice of variables with a new approach to the treatment of particle states in $D$ dimensions for the numerical evaluation of amplitudes, we obtain the analytic expressions with a modest computational effort. Their systematic simplification using multivariate partial-fraction decomposition leads to a particularly compact form. Our results provide all two-loop amplitudes required for the calculation of next-to-next-to-leading order QCD corrections to the production of three jets at hadron colliders in the leading-color approximation.

Figures

Tables

Table 1. Table 1: Color-ordered Feynman rules for vertices with scalar particles explicitly introduced by dimensional reduction of gluons.

	$\frac{i}{\sqrt{2}}{\left(p_2-p_1\right)}^{\mu }{g}_{\left[D_s-6\right]}^{\alpha \beta }$
	$-\frac{i}{2}{g}_{\left[6\right]}^{\mu \nu }{g}_{\left[D_s-6\right]}^{\alpha \beta }$
	$i{g}_{\left[6\right]}^{\mu \nu }{g}_{\left[D_s-6\right]}^{\alpha \beta }$

Equations

$$M=\sum _n{v}_n{M}_n,$$

$$A(q,\bar{q},g,\dots ,g)$$

$$A(q,\bar{q},Q,\bar{Q},g,\dots ,g)$$

$$\displaystyle\begin{split}A(1_{g},2_{g},3_{g},4_{g},5_{g})\big{|}_{\textrm{leading color}}=&\sum_{\sigma\in S_{5}/Z_{5}}\operatorname{Tr}\left(T^{a_{\sigma(1)}}T^{a_{\sigma(2)}}T^{a_{\sigma(3)}}T^{a_{\sigma(4)}}T^{a_{\sigma(5)}}\right)\\ &\times\mathcal{A}({\sigma(1)}_{g},{\sigma(2)}_{g},{\sigma(3)}_{g},{\sigma(4)}_{g},{\sigma(5)}_{g})\,,\end{split}$$

$$\displaystyle\begin{split}A(1_{q},2_{\bar{q}},3_{g},4_{g},5_{g})\big{|}_{\textrm{leading color}}=&\sum_{\sigma\in S_{3}}\left(T^{a_{\sigma(3)}}T^{a_{\sigma(4)}}T^{a_{\sigma(5)}}\right)^{\;\bar{\imath}_{2}}_{i_{1}}\\ &\times\mathcal{A}(1_{q},2_{\bar{q}},\sigma(3)_{g},\sigma(4)_{g},\sigma(5)_{g})\,,\end{split}$$

$$\displaystyle\begin{split}A(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{g})\big{|}_{\textrm{leading color}}=&\,(T^{a_{5}})^{\;\bar{\imath}_{2}}_{i_{3}}\delta^{\;\bar{\imath}_{4}}_{i_{1}}\;\mathcal{A}(1_{q},2_{\bar{q}},5_{g},3_{Q},4_{\bar{Q}})\,\,+\\ &\,(T^{a_{5}})^{\;\bar{\imath}_{4}}_{i_{1}}\delta^{\;\bar{\imath}_{2}}_{i_{3}}\;\mathcal{A}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{g})\,,\end{split}$$

$$\mathcal{A}=g_0^3\left({\mathcal{A}}^{(0)}+\frac{{\alpha }_0}{4\pi }{N}_c{\mathcal{A}}^{(1)}+{\left(\frac{{\alpha }_0}{4\pi }\right)}^2{N_c}^2{\mathcal{A}}^{(2)}+\mathcal{O}({\alpha }_0^3)\right),$$

$$\begin{array}{rl}{\mathcal{A}}^{(1)}& ={\mathcal{A}}^{(1)[N_f^0]}+\frac{N_f}{N_c}{\mathcal{A}}^{(1)[N_f^1]},\\ {\mathcal{A}}^{(2)}& ={\mathcal{A}}^{(2)[N_f^0]}+\frac{N_f}{N_c}{\mathcal{A}}^{(2)[N_f^1]}+{\left(\frac{N_f}{N_c}\right)}^2{\mathcal{A}}^{(2)[N_f^2]}.\end{array}$$

$${\alpha }_0{\mu }_0^{2ϵ}{S}_{ϵ}={\alpha }_s{\mu }^{2ϵ}\left(1-\frac{{\beta }_0}{ϵ}\frac{{\alpha }_s}{4\pi }+\left(\frac{{\beta }_0^2}{{ϵ}^2}-\frac{{\beta }_1}{ϵ}\right){\left(\frac{{\alpha }_s}{4\pi }\right)}^2+\mathcal{O}\left({\alpha }_s^3\right)\right),$$

$${\mathcal{A}}_R=S_{ϵ}^{-\frac{3}{2}}{g}_s^3\left({\mathcal{A}}_R^{(0)}+\frac{{\alpha }_s}{4\pi }{N}_c{\mathcal{A}}_R^{(1)}+{\left(\frac{{\alpha }_s}{4\pi }\right)}^2{N_c}^2{\mathcal{A}}_R^{(2)}+\mathcal{O}({\alpha }_s^3)\right),$$

$$\begin{array}{rl}{\mathcal{A}}_R^{(1)}& ={\mathbf{I}}_{[n]}^{(1)}(ϵ){\mathcal{A}}_R^{(0)}+\mathcal{O}({ϵ}^0),\\ {\mathcal{A}}_R^{(2)}& ={\mathbf{I}}_{[n]}^{(2)}(ϵ){\mathcal{A}}_R^{(0)}+{\mathbf{I}}_{[n]}^{(1)}(ϵ){\mathcal{A}}_R^{(1)}+\mathcal{O}({ϵ}^0).\end{array}$$

$${\mathcal{R}}^{(2)}={\mathcal{A}}_R^{(2)}-{\mathbf{I}}_{[n]}^{(1)}{\mathcal{A}}_R^{(1)}-{\mathbf{I}}_{[n]}^{(2)}{\mathcal{A}}_R^{(0)}+\mathcal{O}(ϵ).$$

$${\mathcal{A}}^{(2)}({\ell }_l)=\sum _{\Gamma \in \Delta }\sum _{i\in M_{\Gamma }\cup S_{\Gamma }}{c}_{\Gamma ,i}\frac{m_{\Gamma ,i}({\ell }_l)}{{\prod }_{j\in P_{\Gamma }}{\rho }_j},$$

$$\sum _{\mathrm{states}}\prod _{i\in T_{\Gamma }}{\mathcal{A}}_i^{(0)}({\ell }_l^{\Gamma })=\sum _{\begin{array}{c}{\Gamma }^{\prime }\ge \Gamma ,\\ i\in M_{{\Gamma }^{\prime }}\cup S_{{\Gamma }^{\prime }}\end{array}}\frac{c_{{\Gamma }^{\prime },i}{m}_{{\Gamma }^{\prime },i}({\ell }_l^{\Gamma })}{{\prod }_{j\in (P_{{\Gamma }^{\prime }}\setminus P_{\Gamma })}{\rho }_j({\ell }_l^{\Gamma })},$$

$${\mathcal{A}}^{(2)}=\sum _{\Gamma \in \Delta }\sum _{i\in M_{\Gamma }}{c}_{\Gamma ,i}{\mathcal{I}}_{\Gamma ,i},$$

$${\mathcal{I}}_{\Gamma ,i}=\int d^D{\ell }_l\frac{m_{\Gamma ,i}({\ell }_l)}{{\prod }_{j\in P_{\Gamma }}{\rho }_j}.$$

$${\mathcal{A}}^{(2)}=\sum _{i\in B}\sum _{k=-4}^0{ϵ}^k{d}_{k,i}{h}_i,$$

$${\mathcal{R}}^{(2)}=\sum _{i\in B}{r}_i{h}_i.$$

$$\vec{s}=\{s_{12},s_{23},s_{34},s_{45},s_{51}\},$$

$${\mathrm{tr}}_5=4i{\epsilon }^{\mu \nu \rho \sigma }{p}_{1\mu }{p}_{2\nu }{p}_{3\rho }{p}_{4\sigma },$$

$$r_i=\frac{n_i}{W^{{\vec{q}}_i}}.$$

$$A_{\alpha }=\{i|i\in A\text{mbox}with{\alpha }_i\ne 0\}.$$

$$W^{\vec{\alpha }}<W^{\vec{\beta }}\leftrightarrow \vec{\alpha }<\vec{\beta }.$$

$${\ell }_1=({\ell }_{1[4]},{\vec{\mu }}_1),{\ell }_2=({\ell }_{2[4]},{\vec{\mu }}_2).$$

$${\vec{\mu }}_1(t)=\frac{1}{2}\left(\begin{array}{c}t+\frac{{\mu }_{11}}{t}\\ t-\frac{{\mu }_{11}}{t}\end{array}\right),{\vec{\mu }}_2(t)=\frac{{\mu }_{12}}{{\mu }_{11}}{\vec{\mu }}_1(t)-\frac{r}{{\mu }_{11}}\frac{1}{2}\left(\begin{array}{c}t-\frac{{\mu }_{11}}{t}\\ t+\frac{{\mu }_{11}}{t}\end{array}\right),$$

$$A(D_s)=\sum _{i=0}^N{\mathcal{K}}_i{D}_s^i,$$

$$A(D_s)=\sum _{i=0}^N{\tilde{\mathcal{K}}}_i(D_s-D_0{)}^i,$$

$$g_{[D_s]}^{\mu \nu }=g_{[D_s-6]}^{\mu \nu }+g_{[6]}^{\mu \nu },g_{[D_s-6]}^{\mu \nu }{g}_{\mu \nu [6]}^{}=0,$$

$$(\gamma_{[D_{s}]}^{\mu})_{a\kappa}^{\,b\lambda}=\left\{\begin{array}[]{ll}\left(\gamma_{[6]}^{\mu}\right)_{a}^{\;b}\,\delta_{\kappa}^{\lambda}\,,&\quad 0\leq\mu\leq 5\,,\\ &\\ \left(\gamma^{\star}_{[6]}\right)_{a}^{\;b}\left(\gamma_{[D_{s}-6]}^{(\mu-6)}\right)_{\kappa}^{\;\lambda}\,,&\quad\mu\geq 6\,,\end{array}\right.$$

$$\mathrm{Tr}\left(\prod _{{\mu }_i\in \mathcal{G}}{\gamma }_{[D_s]}^{{\mu }_i}\right)=\mathrm{Tr}\left(\prod _{{\mu }_i\in \tilde{\mathcal{G}}}{\gamma }_{[D_s-6]}^{{\mu }_i}\right)\cdot \mathrm{Tr}\left(\prod _{{\mu }_i\in \mathcal{G}}{\gamma }_{[6]}^{\mathfrak{I}({\mu }_i)}\right),$$

$$\mathfrak{I}:\mu \to \{\begin{array}{ll}\mu ,& 0\le \mu \le 5,\\ \star ,& \mu \ge 6.\end{array}$$

$$\displaystyle\hskip 64.58313pt=\hskip 25.83325pt\vbox{\hbox{\includegraphics[width=77.49976pt]{figures/Ds-example.pdf}}}\quad+$$

Full text

††preprint: CP3-19-13, FR-PHENO-2019-002, IPhT-19/024

1 Introduction

Scattering amplitudes in QCD provide the basic building blocks for hadron collider phenomenology. Integrated over phase space, they yield theoretical predictions that can be compared to experimental measurements like those performed at the CERN Large Hadron Collider (LHC). These theoretical predictions can be systematically improved through inclusion of higher-order QCD corrections, which require higher-loop scattering amplitudes. The analytic computation of these multi-loop amplitudes is a non-trivial task.

Experience at one loop has shown that the direct numerical computation of scattering amplitudes can avoid difficulties encountered in analytic approaches. However, while numerical evaluation is sufficient for phenomenological applications, compact analytic results are still very useful. Indeed, the numerical evaluation of analytic expressions is generally straightforward to set up and oftentimes more stable and efficient than the purely numerical approach. Furthermore, explicit formulae allow for a detailed analysis of the analytic properties of amplitudes in order to learn about higher-order perturbation theory.

In this article we present the analytic expressions of the two-loop five-parton QCD amplitudes in the leading-color approximation. Recently, significant progress has been made regarding the computation of two-loop multi-particle amplitudes. In the case of five-point QCD amplitudes, the first one to be studied was the leading-color five-gluon amplitude with all helicities positive, initially evaluated numerically Badger:2013gxa and afterwards presented analytically Gehrmann:2015bfy ; Dunbar:2016aux . Building on this result, the all-plus two-loop six- and seven-gluon amplitudes were obtained Dunbar:2016gjb ; Dunbar:2017nfy . By now all leading-color two-loop five-parton amplitudes (i.e., with external gluons and/or massless quarks) have been computed numerically Badger:2017jhb ; Abreu:2017hqn ; Badger:2018gip ; Abreu:2018jgq . Very recently, numerical algorithms were combined with functional reconstruction techniques Peraro:2016wsq resulting in analytic expressions for the planar two-loop five-gluon single-minus helicity amplitude Badger:2018enw and for all two-loop five-gluon helicity amplitudes Abreu:2018zmy . These calculations rely on the availability of the planar two-loop five-point master integrals, which have been given in refs. Papadopoulos:2015jft ; Gehrmann:2018yef . Progress with full-color five-gluon amplitudes Badger:2015lda is gaining momentum with recent results towards the computation of non-planar master integrals Abreu:2018aqd ; Chicherin:2018old . In more conventional approaches, integration-by-parts relations were obtained Boels:2018nrr ; Chawdhry:2018awn which can be used to compute the same type of two-loop five-point amplitudes.

Here we apply a numerical variant Ossola:2006us ; Ellis:2007br ; Giele:2008ve ; Berger:2008sj of the unitarity method Bern:1994zx ; Bern:1994cg ; Bern:1997sc ; Britto:2004nc which has recently been extended to two loops Ita:2015tya ; Abreu:2017idw ; Abreu:2017xsl . Furthermore, we take advantage of computations with exact kinematics through the usage of finite-field arithmetic vonManteuffel:2014ixa ; Peraro:2016wsq in order to functionally reconstruct Peraro:2016wsq the multivariate rational coefficients of a basis of special functions Gehrmann:2018yef . We have already shown that this strategy can be employed to compute amplitudes of relevance to LHC phenomenology, by providing the analytic expressions of all planar two-loop five-gluon amplitudes Abreu:2018zmy . In this article we further improve on the latter work by producing compact analytic results for all leading-color two-loop five-parton amplitudes in QCD. This includes the five-parton processes with five gluons, two quarks and three gluons, and four quarks and one gluon, with zero, one and two light-quark loops. These are all the two-loop amplitudes required for the computation of the leading-color next-to-next-to-leading order (NNLO) QCD corrections to three-jet production at hadron colliders.

A number of new developments allows to obtain these results. First, we set up an efficient method to determine the dependence on the dimension $D_{s}$ associated to the particle states circulating in the loop. We extend the approach of Anger:2018ove and remove a bottleneck of dimensional reconstruction Giele:2008ve ; Ellis:2008ir ; Boughezal:2011br in fermion amplitudes Abreu:2018jgq by analytically precomputing part of the dependence on the dimensional regulator.

Second, we modify the multivariate reconstruction algorithm that we employed in ref. Abreu:2018zmy in order to use an ansatz in terms of Mandelstam variables instead of twistor variables Hodges:2009hk . This is achieved by exploiting the analytic structure of amplitudes with five external massless partons in order to obtain target functions which are rational functions of Mandelstam variables. This has a dramatic impact in reducing the degree of the multivariate polynomials that we reconstruct, allowing us to obtain all analytic expressions with a modest computational effort.

Finally, in order to employ a single finite field for rational reconstruction, we perform a systematic analysis of the analytic structure of the reconstructed functions which takes advantage of two major simplification procedures. As a by-product, the procedure results in rather compact expressions for all helicity amplitudes, which we provide as ancillary files. We first analyze the dimension of the function space spanned by all the pentagon-function coefficients on each amplitude. It turns out that this dimension is an order of magnitude smaller (on average) than the number of pentagon functions. Then we simplify each member of the basis of this space with a multivariate partial-fraction decomposition. All taken into account, the final expressions we present are reduced in byte size by up to two orders of magnitude, with final results that can easily be handled for future analytic and numerical evaluations.

This article is organized as follows. Section 2 is devoted to the description of the numerical calculation of the multi-parton amplitudes. There we precisely define the objects we compute, including the two-loop helicity amplitudes and the finite remainder functions in dimensional regularization. In section 3 our method of handling the $D_{s}$ dependence is presented, and section 4 describes the analytic reconstruction and the simplification of the results by means of multivariate partial fractioning. We discuss the analytic results in section 5 before concluding in section 6. Additional information related to the infrared structure of the amplitudes, Feynman rules and a rationalization of the momentum-space variables is presented in three appendices.

2 Calculation of Multi-Parton Planar Amplitudes

2.1 Multi-Parton Helicity Amplitudes

In this work we compute the five-parton two-loop QCD helicity amplitudes in the leading-color approximation. More precisely, we keep the leading terms in the formal limit of a large number of colors ${N_{c}}$, and scale the number of light flavors ${N_{f}}$ while keeping the ratio ${N_{f}}/{N_{c}}$ fixed.

We evaluate amplitudes in dimensional regularization. Care must be taken when applying dimensional regularization in a numerical approach with external fermions Abreu:2018jgq . In such an approach amplitudes are computed with integer dimensional particle representations, while dimensionally-regulated amplitudes are analytically continued to non-integer dimensions $D=4-2\epsilon$. We will follow the same procedure described in detail in ref. Abreu:2018jgq , which is related to the idea of using suitable projection operators Glover:2004si allowing one to work with Lorentz-invariant objects in the main computational steps. A scattering amplitude $M$ can thus be written as

[TABLE]

where the $M_{n}$ are Lorentz scalars and the $v_{n}$ provide a basis for the $(D_{s}-4)$-dimensional spinor structures. The basis $\{v_{n}\}$ is process dependent, and the loop order at which a given spinor structure begins to contribute depends on the regularization scheme. For concreteness, we work in the ’t Hooft-Veltman (HV) scheme and use the bases of $\{v_{n}\}$ considered in ref. Abreu:2018jgq . For NNLO phenomenology, in cases where two-loop virtual corrections are interfered with a tree amplitude, it is sufficient to compute $M_{0}$ through

[TABLE]

where on the right-hand-side of the equations we compute the traces to project onto the structures $v_{0}$ in eq. (1). This means that the relevant contributions are those in which the open $(D_{s}-4)$-dimensional spinor indices are traced over on each quark line, see figure 1.

While the amplitudes in eq. (2) are easily defined and evaluated analytically, it is important to find an efficient way to compute them in a numerical setup. We discuss our approach in section 3 where we give details of our implementation for the numerical evaluation of amplitudes with fermions. This definition of helicity amplitudes with quarks is consistent with that of ref. Glover:2004si and the prescription given in ref. Anger:2018ove .

When considering amplitudes with two identical quark lines, contributions where index contractions lead to a single trace as in figure 2 should also be considered. Nevertheless, at the level of the finite remainder, amplitudes with identical quarks can be obtained by antisymmetrizing distinct-flavor expressions DeFreitas:2004kmi ; Abreu:2018jgq obtained from eq. (2b). Thus, our results are sufficient for NNLO QCD phenomenological studies of processes involving identical quarks.

The gluon and fermion amplitudes (2) can be decomposed in terms of color structures. We denote the fundamental generators of the $SU({N_{c}})$-group by $(T^{a})^{\;\bar{\jmath}}_{i}$, where the adjoint index $a$ runs over ${N_{c}}^{2}-1$ values and the (anti-) fundamental indices $i$ ($\bar{\jmath}$) run over ${N_{c}}$ values. The fundamental generators are normalized as $\operatorname{Tr}(T^{a}T^{b})=\delta^{ab}$. One can then consider the color decomposition of each process as

[TABLE]

where $S_{n}$ denotes all permutations of $n$ indices and $S_{n}/Z_{n}$ all non-cyclic permutations of $n$ indices. We write the particle type explicitly as a subscript, and all remaining properties of each particle (momentum, helicity, etc.) are implicit.

The $\mathcal{A}$ in eqs. 3, 4 and 5 are called partial amplitudes. For each of the amplitudes considered, all partial amplitudes can be related to the others by exchanging external legs, so only one partial amplitude is independent. These are expanded in a perturbative expansion,

[TABLE]

where $\alpha_{0}=g_{0}^{2}/(4\pi)$ is the bare QCD coupling and $\mathcal{A}^{(k)}$ denotes a $k$-loop partial amplitude. Each $\mathcal{A}^{(k)}$ can be further expanded as a series in powers of $N_{f}/N_{c}$,

[TABLE]

We compute the coefficients $\mathcal{A}^{(k)[N_{f}^{l}]}$, with $0\leq l\leq k\leq 2$, where only planar diagrams contribute as we work in the leading-color approximation. In figs. 3, 4 and 5 we give representative diagrams for each of these contributions.

2.2 Finite Remainder

The bare scattering amplitudes defined in eq. (7) have divergences of ultraviolet and infrared origin. Both can be predicted from lower-loop amplitudes. It is convenient to remove this redundant information and define a finite remainder that contains the genuine two-loop information. There is no unique way to define the remainder, so we now discuss our conventions, with more details given in appendix A.

The renormalized amplitudes can be obtained from their bare counterparts by replacing in eq. (6) the bare QCD coupling $\alpha_{0}$ by the renormalized coupling $\alpha_{s}$ in $D=4-2\epsilon$ dimensions. The two couplings are related by

[TABLE]

which we can use to define the perturbative expansion of the renormalized amplitude,

[TABLE]

where $S_{\epsilon}=(4\pi)^{\epsilon}e^{-\epsilon\gamma_{E}}\,$ and $\alpha_{s}=g_{s}^{2}/(4\pi)$. The $\beta_{i}$ are the coefficients in the perturbative expansion of the QCD $\beta$-function, which we give explicitly in appendix A. Here, $\mu_{0}^{2}$ is the scale introduced in dimensional regularization to keep the coupling dimensionless in the QCD Lagrangian, and $\mu^{2}$ is the renormalization scale. In the following, we set $\mu^{2}=\mu_{0}^{2}=1$ (with arbitrary dimensions).

The renormalized amplitudes $\mathcal{A}^{(k)}_{R}$ have only infrared divergences, which can be determined from lower-loop functions and well known universal factors Catani:1998bh ; Sterman:2002qn ; Becher:2009cu ; Gardi:2009qi . More precisely, we have

[TABLE]

The ${\bf I}^{(k)}_{[n]}$ are operators in color space that become diagonal in the leading-color approximation. They contain some process-specific components, and explicit expressions are given in appendix A.

Using eqs. 8, 9 and 10, we can predict the poles of the two-loop amplitudes we wish to compute. Alternatively, we can use them to define a finite remainder $\mathcal{R}^{(2)}$, according to

[TABLE]

In this expression, we extend the expansion of ${\bf I}_{[n]}^{(1)}\mathcal{A}_{R}^{(1)}$ and ${\bf I}_{[n]}^{(2)}\mathcal{A}_{R}^{(0)}$ to also include terms of order $\epsilon^{0}$. This subtracts non-trivial contributions from the finite term of $\mathcal{A}_{R}^{(2)}$ that are related to the lower-loop amplitudes. In this paper, we directly compute analytic expressions for the remainders $\mathcal{R}^{(2)}$, from which one can then recover the full two-loop bare amplitude $\mathcal{A}^{(2)}$, see eq. (61) for the explicit relation.

2.3 Numerical Unitarity

The first step in the analytic reconstruction of the remainders defined in eq. (11) is to evaluate the bare amplitudes numerically. To this end, we employ numerical unitarity Ita:2015tya ; Abreu:2017xsl ; Abreu:2017idw in order to evaluate the coefficients of a decomposition of the amplitude into a linear combination of master integrals. First, we consider a decomposition of the integrand of the amplitude $\mathcal{A}^{(2)}(\ell_{l})$ in terms of master integrands and surface terms (we use $\ell_{l}$ to denote the loop momenta). The former correspond to master integrals and the latter integrate to zero. Specifically, we have

[TABLE]

where $\Delta$ is the set of all propagator structures $\Gamma$, $P_{\Gamma}$ the set of all inverse propagators $\rho_{j}$ in $\Gamma$, and $M_{\Gamma}$ and $S_{\Gamma}$ denote the corresponding sets of master integrands and surface terms.

To compute the decomposition of the amplitude in terms of master integrals we must evaluate the coefficients $c_{\Gamma,i}$, with $i\in M_{\Gamma}$. These are rational functions of spinor components of the external momenta and the dimensional regulator $\epsilon$. We determine them using the standard approach in numerical unitarity. We first build a system of linear equations through sampling of on-shell values of the loop momenta, i.e., where $\ell_{l}\to\ell_{l}^{\Gamma}$ with $\ell_{l}^{\Gamma}$ such that $\rho_{j}(\ell_{l}^{\Gamma})=0$ for $j\in P_{\Gamma}$. The leading contribution to eq. (12) in this limit factorizes into products of tree amplitudes

[TABLE]

where we label the set of tree amplitudes associated with the vertices in the diagram corresponding to $\Gamma$ by $T_{\Gamma}$, and the sum on the right-hand side runs over the propagator structures $\Gamma^{\prime}$ such that $P_{\Gamma}\subseteq P_{\Gamma^{\prime}}$. The sum on the left-hand side runs over the (scheme dependent) physical states of each internal line of $\Gamma$. At two loops there are also subleading contributions in the limit $\ell_{l}\to\ell_{l}^{\Gamma}$ which can easily be dealt with even though no factorization theorem is known Abreu:2017idw (for a recent related study on these contributions see Baumeister:2019rmh ). The coefficients $c_{\Gamma,i}$ can then be obtained at a given phase-space point by solving the linear system in eq. (13) numerically. The numerical evaluations are performed using finite-field arithmetic, allowing to efficiently obtain exact results for rational phase-space points, circumventing problems of numerical instabilities. This is key for the task of functional reconstruction.

2.4 Pentagon-Function Decomposition

Once the coefficients $c_{\Gamma,i}$ in eq. (12) are computed, we obtain a decomposition of the amplitude into master integrals,

[TABLE]

where the $\mathcal{I}_{\Gamma,i}$ are the master integrals,

[TABLE]

For planar five-parton amplitudes, a basis of master integrals has been computed Papadopoulos:2015jft ; Gehrmann:2018yef . The master integrals evaluate to linear combinations of so-called multiple polylogarithms (MPLs), which can be numerically evaluated using available programs (e.g. Vollinga:2004sn ). Once numerical values for the coefficients $c_{\Gamma,i}$ have been computed (using for example the approach described in the previous subsection) one can obtain numerical values for the amplitudes.

The MPLs are a class of special functions with only logarithmic singularities that can be equipped with algebraic structures that allow one to algorithmically find relations between them Goncharov:2010jf ; Duhr:2011zq ; Duhr:2012fh . One can then construct a basis for the space of MPLs relevant for five-parton scattering amplitudes, and this was achieved in Gehrmann:2018yef where the so-called pentagon functions were introduced. Further, MPLs are equipped with a notion of weight, which can be used to organize the pentagon functions. Indeed, there are no relations between pentagon functions of different weight so we can separate the space of pentagon functions into subspaces of different weights. For two-loop amplitudes, we need functions of at most weight 4. In the following, we denote the pentagon functions by $\{h_{i}\}_{i\in B}$ and the associated set of labels $B$.

After expansion in epsilon, the amplitude can thus be expressed in terms of pentagon functions

[TABLE]

where we use the fact that the poles in $\epsilon$ of two-loop amplitudes are at most $\mathcal{O}(\epsilon^{-4})$ and the $d_{k,i}$ are rational functions of the external data. The motivation for this decomposition is that order by order in epsilon the master integrals satisfy more relations than integration-by-parts relations and these are manifested by the pentagon function decomposition. The same basis $B$ can be used to express the terms ${\bf I}_{[n]}^{(1)}\mathcal{A}_{R}^{(1)}$ and ${\bf I}_{[n]}^{(2)}\mathcal{A}_{R}^{(0)}$ in eq. (11), and we can thus decompose the remainder in terms of pentagon functions:

[TABLE]

Again, we suppress the dependence of the algebraic coefficient functions $r_{i}$ and the pentagon functions $h_{i}$ on the external data. The remainders, and consequently the coefficient functions $r_{i}$, depend on particle types and helicities as well as the number of flavors ${N_{f}}$.

For convenience we use the convention that the set of pentagon functions $\{h_{i}\}_{i\in B}$ not only includes genuine functions, such as $\log(-s_{12})$, but also constants with weight, such as $\pi^{2}$, that correspond to the pentagon functions evaluated at specific points. These are boundary conditions that are specific to the kinematic region where the amplitude is evaluated. In this paper we specialize our discussion to the Euclidean region. In ref. Gehrmann:2018yef , the pentagon functions have been continued to all relevant physical regions, and it is thus rather straightforward to extend our results to any physical region.

2.5 Analytic Properties of Coefficient Functions

The coefficient functions $\{r_{i}\}_{i\in B}$ introduced in eq. (17) have a universal pole structure Abreu:2018zmy which we find to be expressible in terms of the so-called planar alphabet $A$ of the pentagon functions Gehrmann:2018yef . The planar alphabet $A$ is given in terms of 26 letters $\{W_{i}\}_{i\in A}$ which in turn are rational functions of the independent Mandelstam variables $s_{ij}=(p_{i}+p_{j})^{2}$,

[TABLE]

and the parity-odd contraction of four momenta

[TABLE]

with the Levi-Civita symbol $\varepsilon_{\mu\nu\rho\sigma}$ and the convention $\varepsilon_{0123}=1$. The square of ${\rm tr}_{5}$ gives the five-point Gram determinant $\Delta_{5}={\rm tr}_{5}^{2}$. The letters can be further grouped into parity even and odd letters, $A^{+}$ and $A^{-}$ respectively, with $A=A^{+}\cup A^{-}$. (Parity here refers to the transformation properties of the letters under a parity transformation in momentum space.) While ${\rm tr}_{5}$ is parity odd, all Mandelstam variables $\vec{s}$ are parity even. The functions $r_{i}$ are rational functions in the variables ${\vec{s}}$ and ${\rm tr}_{5}$, whose denominators are monomials in the letters,

[TABLE]

The vector of exponents $\vec{q}_{i}=\{q_{i,1},...q_{i,26}\}$ differs between the coefficient functions, and for a given coefficient function $r_{i}$ not all letters contribute. The pattern of contributing letters is linked to the pentagon functions and is helicity dependent Abreu:2018zmy . We will often refer to the set of independent factors in a monomial such as $W^{\vec{\alpha}}$ as $A_{\alpha}$,

[TABLE]

Finally, we will have to specify the notion of a monomial ordering on the exponent vectors,

[TABLE]

We will often use the lexicographic monomial ordering, which amounts to comparing the size of the leading entries in the vectors and, if they are equal, comparing the following entries etc. We refer to text books such as ref. cox2013ideals for further information concerning the concepts of monomial ordering in the context of polynomial-division algorithms.

2.6 Amplitude Evaluation

First we summarize the tools employed to obtain the results we present. The propagator structures $\Delta$ in eq. (12) are obtained by generating colored cut diagrams with QGRAF Nogueira:1991ex followed by color decomposition performed in Mathematica according to Ochirov:2016ewn ; FermionColour . We carry out the decomposition into master and surface integrands as described in refs. Abreu:2017xsl ; Abreu:2017hqn ; Abreu:2018jgq ; Abreu:2018zmy , where we used SINGULAR DGPS to obtain the unitarity-compatible surface terms. The master integral coefficients in eq. (12) are evaluated over finite fields,111With cardinalities of order $\mathcal{O}(2^{31})$. employing our C++ framework for multi-loop numerical unitarity. We use Givaro Givaro for basic finite-fields arithmetic, and improve the multiplication speed with a custom implementation of a Barrett reduction Barrett1987 ; HoevenLQ14 . The $D_{s}$-dimensional tree-level amplitude products appearing on the left-hand-side of eq. (13) are evaluated through off-shell recursion Berends:1987me , and the corresponding linear system of equations are solved by PLU factorization and back substitution.

Furthermore, we document some technical aspects of the computations. The transformation from master integral coefficients in eq. (14) into the coefficients of pentagon functions in eq. (17) is accomplished with our C++ representation of the master integrals in terms of pentagon functions Gehrmann:2018yef . We refine the treatment of square roots appearing from solving the quadratic on-shell conditions on the loop momenta compared to the previous implementation described in ref. Abreu:2018jgq . As before, we rotate the $D$-dimensional components of the loop momenta into a six-dimensional subspace,

[TABLE]

We first use the algorithm of Abreu:2017hqn to solve the on-shell conditions such that the $\mu_{ij}$ are rational in the input parameters. In order to represent these momenta in a 6-dimensional embedding we proceed as follows. Without loss of generality, we choose to use an alternating metric signature $(+,-,+,-,\ldots)$, and parametrize the two-dimensional loop-momentum components $\vec{\mu}_{1}$ and $\vec{\mu}_{2}$ as follows:

[TABLE]

where $t$ is a free dimensionful parameter that leaves the scalar products $r=\sqrt{\mu_{12}^{2}-\mu_{11}\mu_{22}}$ and $\mu_{ij}=\mu_{i}^{1}\mu_{j}^{1}-\mu_{i}^{2}\mu_{j}^{2}$ invariant. We perform numerical computations with the external kinematic data $\{p_{i}\}_{i=1,5}$ taking values in a finite field. While we do not require loop momenta to take values in the same number field, eq. (24) guarantees that their components take values in an algebra generated by the basis $\{1,r\}$ over the same number field. The above parametrization is an improvement compared to the basis of four elements employed in ref. Abreu:2018jgq .

3 $D_{s}$ Dependence from Dimensional Reduction

The amplitudes $A$ defined in eq. (2) are polynomials in $D_{s}$,

[TABLE]

where $N$, the maximal power of $D_{s}$, varies depending on the process, the loop order, and the choice of tensor structure in eq. (1). For the amplitudes considered in this paper $N\leq 2$. In this section we will suppress all arguments of $A$ and only keep track of the dependence on $D_{s}$. In a numerical framework, $A(D_{s})$ can only be evaluated for integer $D_{s}$ values for which the particle states are well defined. To be able to set $D_{s}=4-2\epsilon$ in the HV scheme, the knowledge of the coefficients $\mathcal{K}_{i}$ is required. One way to obtain them is to reconstruct the polynomial (25) from a sample of $(N+1)$ integer values of $D_{s}$. This procedure is known as dimensional reconstruction Giele:2008ve and has previously been applied in Ellis:2008ir ; Boughezal:2011br ; Abreu:2017xsl ; Abreu:2017hqn .

While being generic and straightforward to implement, this approach has drawbacks which become particularly evident in amplitudes with fermions. The dimension of the spinor representation scales exponentially (as $2^{D_{s}/2}$) with $D_{s}$, as opposed to the linear scaling of the vector representation. Furthermore, the external spinor states with definite helicity can be embedded consistently only for even values of $D_{s}$, which pushes the sample values higher compared to the case of vector particles in the loops. Beyond the obvious detrimental effect on the numerical complexity, the dimensionality of the spinor representation determines the number of terms entering the evaluation of the traces to obtain the helicity amplitudes through eq. (2). For the case of amplitudes with multiple external quark pairs this makes the computation of traces unnecessarily time consuming.

These considerations motivate the search for more efficient alternatives to dimensional reconstruction. Here we employ one such alternative, based on the idea of dimensional reduction, which has recently been presented in ref. Anger:2018ove and already applied to the computation of one-loop amplitudes in ref. Anger:2017glm . In the remainder of this section we give a brief overview of this method and refer the reader to ref. Anger:2018ove for more technical details.

We start by rearranging eq. (25) in the following way:

[TABLE]

where $D_{0}$ is some base dimension, and the coefficients $\tilde{\mathcal{K}}_{i}$ can be obtained by a linear transformation of the coefficients $\mathcal{K}_{i}$ in eq. (25). It turns out that, given a suitable choice of $D_{0}$, the dependence of $A(D_{s})$ on degrees of freedom higher than $D_{0}$ can be captured in a kinematic-independent way. This observation allows one to analytically separate this dependence, and thus evaluate each coefficient $\tilde{\mathcal{K}}_{i}$ directly. Furthermore, these evaluations are then performed in the base dimension $D_{0}$, resulting in spinor representations of much lower dimensionality compared to those encountered in the framework of dimensional reconstruction.

For reasons that will become clear shortly, one chooses the base dimension $D_{0}$ to be the minimal dimension which allows to embed all loop-momentum components without introducing new relations. For two-loop amplitudes we have $D_{0}=6$, and we shall specialize to this case henceforth. We write the metric tensor as a direct sum,

[TABLE]

and the gamma matrices as a direct product (see e.g. Collins:1984xc ; Kreuzer:susylectures ),

[TABLE]

where $\gamma^{\star}_{[6]}$ is a six-dimensional analogue of $\gamma_{5}$ in four dimensions, i.e. $\{\gamma^{\star}_{[6]},\gamma_{[6]}^{\mu}\}=0$ for all $\mu\in\{0,5\}$, and $(\gamma^{\star}_{[6]})^{2}=1$. The product representation allows us to factorize any chain of gamma matrices into a product of $6$- and $D_{s}-6$-dimensional gamma-matrix chains. Then, using the fact that the trace of a direct product of two matrices is the product of their traces, we can split the traces required to obtain the coefficients of tensor structures in eq. (1) as follows:

[TABLE]

where the product on the left-hand side is over a sequence $\mathcal{G}$ of $D_{s}$-dimensional Lorentz indices, $\tilde{\mathcal{G}}=\{\mu_{i}\in\mathcal{G}~{}|~{}\mu_{i}\geq 6\}$, and the map $\mathfrak{I}$ is defined as

[TABLE]

The traces of $\prod_{\mu_{i}\in\tilde{\mathcal{G}}}\gamma^{\mu_{i}}_{[D_{s}-6]}$ can be evaluated analytically using well-known Clifford algebra identities, which produce sums of products of $g^{\mu_{i}\mu_{j}}_{[D_{s}-6]}$. The crucial observation is that the only object to be contracted with the indices beyond $(D_{s}-6)$ is $g^{\mu\nu}_{[D_{s}-6]}$. This is ensured by our choice of the base dimension. These indices then always contribute terms of the form $g^{\mu}_{[D_{s}-6]\mu}=(D_{s}-6)$, generating contributions to the coefficients $\tilde{\mathcal{K}}_{i}$ with $i>0$ in eq. (26). At this point, all degrees of freedom beyond $(D_{s}-6)$ are traded for polynomials in $(D_{s}-6)$ with integer factors, and the coefficients $\tilde{\mathcal{K}}_{i}$ are expressed in terms of six-dimensional objects only.

From a Feynman diagrammatic perspective, the contributions to the coefficients of the polynomial in $(D_{s}-6)$ can be represented by introducing a scalar particle. The Feynman rules associated to this particle can be readily derived from dimensional reduction of the original Feynman rules of the theory Bern:2002zk ; Badger:2013gxa ; Giele:2008ve ; Anger:2018ove by applying the relations (27) and (28). We list them in the appendix B for convenience. To illustrate this procedure we now give an example. Consider a decomposition of a Feynman diagram with $D_{s}$-dimensional particles on the left-hand side of figure 6. The four non-vanishing contributions after evaluating partial traces and contracting all $(D_{s}-6)$-dimensional indices are shown on the right-hand side, where the scalars are introduced to represent what remains of these contractions.

We would like to conclude this section with some remarks. First, the remaining non-trivial traces of eq. (29), i.e. those of $\prod_{k}\gamma^{\mu_{k}}_{[6]}$, cannot be simplified generically as the indices appear contracted with the loop momenta at the integrand level. We evaluate them by the direct summation over a specially constructed set of external states (in an analogous way to the sums performed in Abreu:2018jgq for dimensional reconstruction). Secondly, in the absence of fermions in the loops, this method coincides with the so called six-dimensional formalism employed in refs. Badger:2013gxa ; Badger:2017jhb , and can thus be viewed as an extension thereof to amplitudes with fermions. Finally, we note that the method presented in this section can be straightforwardly generalized to higher number of loops by adjusting the base dimension $D_{0}$, as well as to the extraction of coefficients of different tensor structures in eq. (1).

4 The Analytic Structure of Five-Parton Amplitudes

We apply three methods to reduce the complexity of the remainder functions and facilitate their reconstruction from numerical samples. First, we introduce a particular ansatz for the coefficients $r_{i}$ of eq. (17) which simplifies the reconstruction procedure by decreasing the required number of sample points. Next, we comment on how to remove redundancies in the large number of rational functions in a scattering amplitude. Finally, we reduce the difficulty of rationally reconstructing the finite-field result by introducing an algorithm to uniquely obtain a multivariate partial-fraction decomposition of the analytic expressions.

4.1 Analytic Reconstruction

In order to perform the reconstruction of the five-parton amplitudes from samples over finite fields, we employ a modification of the algorithm of Peraro:2016wsq . Our aim is to reconstruct the functions $r_{i}$ in terms of $\{\vec{s},{\rm tr}_{5}\}$, with $\vec{s}$ the five independent Mandelstam invariants defined in (18) and the ${\rm tr}_{5}$ defined in (19), while using a suitable rational parametrization of the phase space of five massless particles.

Let us start by introducing such a parametrization of the external kinematics. We keep four of the Mandelstam invariants $\{s_{12},s_{23},s_{45},s_{51}\}$ and introduce an auxiliary variable $x$. More explicitly, we consider the Mandelstam invariants as functions of this new set of variables,

[TABLE]

and, rescaling the momenta such that $s_{12}=1$, $s_{34}$ is defined in terms of $x$ as

[TABLE]

This parametrization follows from the twistor matrix given in appendix C and so naturally rationalizes ${\rm tr}_{5}$. For any (generic) fixed value of the invariants $\vec{s}$ there are two corresponding values of $\{x,\bar{x}\}$, corresponding to the two solutions to the quadratic equation in (32). This manifests the known fact that it is not possible to rationally parametrize the phase space in terms of Mandelstam variables. These two different points in phase space are in fact parity conjugates, and correspond to opposite signs of ${\rm tr}_{5}$ (which is the square root of the discriminant of the quadratic equation in (32)). Under this transformation,

[TABLE]

Having established the set of variables that rationalize the phase space of five massless particles, we next discuss our ansatz for the rational coefficients $r_{i}$ and the sampling procedure to fix the parameters in the ansatz. We decompose each function $r_{i}$ into parity odd and even parts,

[TABLE]

As ${\rm tr}_{5}{}^{2}$ is manifestly polynomial in the invariants, it is clear that no higher powers in ${\rm tr}_{5}$ are required. The motivation for choosing such a representation is two-fold. First, a generic expression of this form has a higher polynomial degree when expressed in terms of the twistor variables, as the relation is non-linear. Consequently, as we shall observe, expressing the finite remainders in this form allows for a more efficient analytic reconstruction. Second, the Mandelstam variables make manifest the dependence on the external momenta, and therefore also any possible symmetries related to exchanges of external legs.

While the functions $r_{i}^{\pm}$ are rational in the Mandelstam invariants $\vec{s}$, their definition through eq. (34) only allows them to be evaluated via the functions $r_{i}$ which, due to the presence of ${\rm tr}_{5}$, are not rational in the invariants themselves. For this reason, we first isolate the even and odd parts $r_{i}^{\pm}$ and reconstruct them in terms of the Mandelstam variables $\vec{s}$. As the parity conjugate points $\{s_{23},s_{45},s_{51},x\}$ and $\{s_{23},s_{45},s_{51},\bar{x}\}$ correspond to the same value of the invariants $\vec{s}$, we can use this pair of points to evaluate the $r_{i}^{\pm}$ on a single point $\vec{s}$,

[TABLE]

A further simplification to the reconstruction procedure is achieved by conjecturing that the ansatz for the denominators given in eq. (20) extends to the $r_{i}^{\pm}$, with the extra requirement that no odd letters appear in the denominator. That is, we conjecture that

[TABLE]

where the $n_{i}^{\pm}(\vec{s}\,)$ are polynomials in the invariants and the contributing $W_{j}(\vec{s}\,)$ are polynomials in the $\vec{s}$ only (and not ${\rm tr}_{5}$). To test this ansatz and determine the exponent vector $\vec{q}_{i}$, we consider the functions $r_{i}^{\pm}$ on a ‘univariate slice’ where the variables depend on a single parameter $t$, $\{s_{23}(t),s_{45}(t),s_{51}(t),x(t)\}$. We then reconstruct the univariate functions $r_{i}^{\pm}(t)$ using Thiele’s method and match the denominators with a monomial in the letters $W_{i}(\vec{s}\,[t])$ in order to obtain the exponent vector $\vec{q}_{i}$. In order to recover this information from the univariate slice, we must choose a curve on which all letters $W_{i}(\vec{s}\,[t])$ are distinct functions of $t$. Furthermore, if we require that all invariants are linear in $t$, this procedure will also tell us the total degree of the numerator polynomials $n_{i}^{\pm}(\vec{s}\,)$ in the Mandelstam variables. Such a curve is given by, e.g.,

[TABLE]

where the variables $c_{i}$ and $d_{i}$ take fixed, arbitrary values.222It is much more transparent in the parametrization of Abreu:2018zmy as to why $s_{34}(t)$ is then linear in $t$. With this procedure, we confirmed the ansatz of eq. (37) and thus determined all the denominators of the rational functions $r_{i}^{\pm}$.

The problem of reconstructing the $r_{i}^{\pm}$ is now reduced to the reconstruction of the polynomials $n^{\pm}_{i}$. We perform the analytic reconstruction of $n_{i}^{\pm}(\vec{s})$ directly in terms of the Mandelstam variables. This is achieved with a modified form of the recursive Newton method Peraro:2016wsq . If we consider a univariate polynomial in the variable $z$, the method makes an ansatz which is adapted to the choice of sample points $\{z_{1},\ldots,z_{n}\}$. Specifically, the polynomial is written as a linear combination of the basis polynomials

[TABLE]

As the $p_{j}(z_{k})$ vanish for $k<j$, when solving for the coefficients of this basis the linear system is triangular by construction. Once one establishes the coefficients of these basis elements, it is trivial to rewrite the polynomial in terms of monomials in $z$. In the multivariate case, one can apply this strategy recursively by singling out one variable and letting the coefficients be polynomials in the remaining variables Peraro:2016wsq .

In our approach we also use the observation that it is not strictly required that the arguments of the $p_{j}$ be the input parameters. In fact, we can choose them to be a function, in our case $s_{34}=s_{34}(s_{23},s_{45},s_{51},x)$. Specifically, in the first step of the recursion, we write our numerator polynomial as

[TABLE]

where

[TABLE]

the value $R$ is obtained from the highest degree term $\sim t^{R}$ found in the univariate slice, and the $n_{i,j}^{\pm}$ are polynomials. To be able to evaluate the $n_{i,j}^{\pm}$ at a point in the recursive approach, we evaluate the function $n_{i}^{\pm}$ at $R+1$ values of $x$, and solve the associated triangular linear system.

To illustrate the advantage of the approach outlined above compared with the strategy used in Abreu:2018zmy , we show in figure 7 the polynomial-degree drop for the two most complex remainders, separating the different ${N_{f}}$ contributions. This observation is common to all five-parton remainders at $N_{f}^{0}$ and $N_{f}^{1}$: we observe a drop of the polynomial degree in the numerators $n_{i}^{\pm}(\vec{s})$ by 40%-50% when written in terms of Mandelstam variables as compared to the twistor variables. The efficiency of the analytic reconstruction algorithm scales as a binomial $\left(\begin{subarray}{c}n+R\\ n\end{subarray}\right)$ in the number of variables $n$ and the polynomial degree $R$. Such a drop in the polynomial degree thus implies an important improvement in efficiency. The $N_{f}^{2}$ remainders are trivial to reconstruct in both approaches.

Finally, as an indication for the complexity of the analytic reconstruction we show the numerator degrees for the $N_{f}^{0}$ contribution of all the planar five-parton remainders in figure 8. The figure displays the expected reduction in complexity of the remainders for external fermions compared to gluons.

4.2 Basis of Rational Coefficients

It was noted in Abreu:2018aqd ; Chicherin:2018yne ; Abreu:2019rpt ; Chicherin:2019xeg that the space of independent rational functions appearing in the $\mathcal{N}=4$ super Yang-Mills (SYM) and $\mathcal{N}=8$ super gravity two-loop five-point amplitudes is much smaller than the space of possible independent pentagon functions. This motivates us to apply the same analysis here, both in order to allow performing minimal reconstruction work and as a way to obtain simpler analytic expressions. These linear dependencies between the rational functions allow us to express the remainder of each amplitude as

[TABLE]

where $r_{i}$ is a vector of rational coefficients and $K$ denotes the dimension of the space they span. The index $j$ runs over a subspace $\hat{B}$ of the space $B$ of all pentagon functions introduced in eq. (17), i.e., $\hat{B}\subseteq B$. The constant matrix $M$, the rational functions $r_{i}$ and the index sets $K$ and $\hat{B}$ are helicity-amplitude dependent.

The basis of rational functions $r_{i}$ is not unique and we choose it by taking the linearly independent subset with the lowest total polynomial degree. The associated matrix $M$ can be computed numerically in the finite field used for the analytic reconstruction, using as many evaluations as its rank. In order to rationally reconstruct $M$ we only require one finite field of cardinality $\mathcal{O}(2^{31})$, apart from two cases where we combine information from two finite fields (see more details in section 5).

4.3 Partial Fractions

Having chosen a basis of rational coefficients, we study the $\{r_{i}^{\pm}\}_{i\in K}$ in a finite field, obtained with the strategy described in section (4.1). The next step is to lift the finite-field result to the rational numbers, i.e., to rationally reconstruct the result. The correctness of the rational reconstruction can be verified by comparison against a single evaluation on a different finite field. In general, the success of the rational reconstruction is guaranteed if the numerator and denominator of the original rational number satisfy a bound dependent on the finite-field cardinality Wang:1981:PAU:800206.806398 ; vonManteuffel:2014ixa ; Peraro:2016wsq , i.e., rational numbers with numerators and denominators with smaller absolute values are more easily reconstructed from their finite-field images. As such, it is beneficial if we can find a form of the $r_{i}^{\pm}$ where the rational numbers involved have ‘simple’ numerators and denominators. To this end, we will employ a multivariate partial-fraction decomposition of the $r_{i}^{\pm}$, which will also have the side effect of simplifying their analytic form as a whole. Indeed, the partial-fraction decomposition manifests the singularities of the coefficients, which are controlled by physical properties such as the factorization of the amplitude in specific limits.

We will employ the partial-fraction decomposition introduced by Leĭnartas leinartas1978factorization ; raichev2012leinartas . This decomposition has recently been used for bringing differential equations into canonical form Meyer:2016slj , as well as for the analytic reduction of loop integrands Zhang:2012ce ; Mastrolia:2012an . A Leĭnartas decomposition of our basis functions is given by

[TABLE]

That is, the $r_{i}^{\pm}$ are written as a linear combination of rational functions which we label by the exponents $\alpha$ of the denominators. The set of all $\vec{\alpha}$ for a given $r_{i}^{\pm}$ is denoted by $D_{i}$ (we combine the exponents for the even and odd coefficients). The $\{n_{i,\vec{\alpha}}^{\pm}\}_{\vec{\alpha}\in D_{i}}$ are the numerators associated with the rational function labelled by $\vec{\alpha}$, which are polynomials in the $\{\vec{s}\}$ variables. We construct the Leĭnartas decomposition (42) coefficient-by-coefficient, such that the exponents in $\vec{\alpha}\in D_{i}$ are non-vanishing only for the letters appearing in the starting denominator $W^{\vec{q}_{i}}$, see eq. (37), which we recall are polynomials in the $\vec{s}$. The particular algorithm which we employ to uniquely specify the Leĭnartas decomposition is related to the techniques of Smirnov:2005ky .

A Leĭnartas decomposition (42) is characterized by exponents of the denominator factors. All denominators $W^{\vec{\alpha}}$ in a Leĭnartas decomposition are ‘algebraically independent’, in that there exists no non-zero polynomial $P(w)$, such that it vanishes upon inserting the denominator factors,

[TABLE]

where $A_{\alpha}$ denotes the set of letters in $W^{\vec{\alpha}}$, see eq. (21). Let us make a few comments about such a decomposition. First, individual terms in eq. (42) may have factors raised to higher degree than when expressed over a common denominator as in eq. (37). Second, no term can have more denominator factors than the number of variables $\vec{s}$ as these factors would necessarily be algebraically dependent. Also, we note that the Leĭnartas decomposition is different to an iterated univariate partial-fraction decomposition in that it maintains the original denominator structures of the input function, i.e., in our case, of eq. (37). Finally, in ref. raichev2012leinartas it is noted that such a decomposition is in general not unique, due to how the algebraic dependence relations (43) are resolved. Let us illustrate this point by considering

[TABLE]

which clearly vanishes upon using the definitions of the letters, $W_{3}=s_{34},W_{4}=s_{45}$ and $W_{6}=s_{34}+s_{45}$. The dependence relation implies a relation between multiple terms with different denominator factors,

[TABLE]

which can be used to remove either of the three denominator factors on the right-hand side. This arbitrariness in the implementation of the Leĭnartas decomposition will in general lead to expressions that are not as compact as possible if inconsistent choices are made across different terms being decomposed. In the following we describe an approach that produces a unique Leĭnartas decomposition by providing a global prescription to expand algebraic dependencies of the denominators of a given coefficient $r_{i}^{\pm}$.

Our approach makes use of multivariate polynomial division and Gröbner basis techniques. In summary, we reinterpret a rational function in a polynomial way in order to employ division by a Gröbner basis to find a canonical form of a polynomial which is subject to a series of constraints. The set of constraints among the variables form a generating set of an ideal and if one divides a polynomial in these variables by a Gröbner basis of the ideal, the remainder is unique and ‘minimal’ with respect to the relations. In the rest of this section we describe these steps in more detail. We refer the reader to ref. cox2013ideals for a pedagogical discussion of all the required concepts in polynomial-division algorithms.

Our starting point is a rational function $r_{i}^{\pm}$ of the form (37), and we introduce an auxiliary variable for each denominator factor $W_{j}$ for $j\in A_{q_{i}}$, where we recall $A_{q_{i}}$ denotes the set of labels of letters appearing in the denominator of $r_{i}^{\pm}$. More precisely we introduce variables $Q_{j}$ to which we associate the constraint polynomial $C_{j}(Q_{j},\vec{s}\,)$,

[TABLE]

Setting all $C_{j}=0$ (i.e., working in the equivalence class of the ideal generated by $\{C_{j}\}_{j\in A_{q_{i}}}$) imposes the constraint that multiplication by $Q_{j}$ is equivalent to division by $W_{j}$. As such, subject to this constraint we can express eq. (37) in a polynomial fashion as

[TABLE]

where the notation ’$\sim$’ emphasises that the relation holds modulo the ideal $\{C_{j}\}_{j\in A_{q_{i}}}$. In order to uniquely implement the constraints $C_{j}=0$, we first divide the right-hand side of eq. (47) by a Gröbner basis of the constraint polynomials $\{C_{j}\}_{j\in A_{q_{i}}}$.333It is important to note that Gröbner basis division works over any field. Hence, the algorithm can be be applied before or after rational reconstruction. After the division, the original polynomial is rewritten as a linear combination of the $C_{j}$ and a unique remainder, which is a polynomial in the $Q_{j}$ and the $\vec{s}$. We then impose $C_{j}=0$, leaving this remainder as the result of the partial-fraction decomposition. Finally, we replace the $Q_{j}$ by $1/W_{j}$ and recover an expression of the form of eq. (42). We note that the Gröbner basis of the ideal $\{C_{j}\}_{j\in A_{q_{i}}}$ depends on a monomial ordering of the variables $\{Q_{j}\}_{j\in A_{q_{i}}}$ and $\{\vec{s}\}$, which in turn affects the final expression. Effectively, by choosing different orderings one can specify which of the denominator factors $\{W_{j}\}_{j\in A_{q_{i}}}$ one prefers in the final expression. Once this choice has been made, the result of the procedure is unique.

Let us now comment on why this simple polynomial division technique achieves a minimal Leĭnartas decomposition. First, it is important to recall that the remainder of the division by a Gröbner basis is minimal in the sense that it cannot be written in terms of a member of the ideal and a new remainder without that remainder having higher polynomial degree. That is, subject to the choice of monomial ordering, the remainder has the lowest possible polynomial degree. With this in mind, it is clear why the division by the constraint polynomials (46) will maximally cancel numerator and denominator: since the second term of the constraint polynomial is $1$, which is the term with the lowest possible degree (in any monomial orderings), all possible cancellations between numerator and denominator will occur, as that ensures they will result in the lowest polynomial degree for the remainder.

Another important question is how this division results in a set of denominator structures with no algebraic dependencies, see eq. (43), and how the relations are guaranteed to be resolved in a consistent way, see the discussion around eq. (45). The fact that algebraic dependencies are implemented follows from a key property of division by a Gröbner basis, which is that it will also apply non-trivial relations implied by the original set of relations. As already argued in raichev2012leinartas , algebraic dependencies induce such non-trivial relations. Indeed, if the set of denominator factors of a given term

[TABLE]

are algebraically dependent, then by definition there exists a non-zero polynomial $P(w)$ such that P\big{(}\{W_{i}\}_{i\in A_{\alpha}}(\vec{s})\big{)}=0, which is trivially a member of the ideal of constraints. We associate the same monomial ordering to the $w$ variables as we do to the $Q_{j}$ and write the polynomial $P(w)$ in the form

[TABLE]

That is, we distinguish the term of lowest monomial degree $\vec{\beta}$ (we refer to the notation introduced in eq. (22)). Let us now multiply the $Q^{\vec{\alpha}}$ in eq. (48) by this polynomial and $Q^{\vec{\beta}}$. We get

[TABLE]

where we introduced the function $\theta(\vec{\alpha})$ which yields a vector with the negative entries removed. This shows that relations such as that of eq. (45) are now also in the ideal, and we do not have to construct them explicitly. Furthermore, given that by construction $\vec{\gamma}>\vec{\beta}$, we have

[TABLE]

That is, for each term in the sum over $\gamma$, the $Q$-dependent part of the monomial is always of lower degree than the original denominator structure $Q^{\vec{\alpha}}$. For any monomial ordering where the $Q_{j}$ are sorted before the $\vec{s}$, the term $c_{\vec{\beta}}\,Q^{\vec{\alpha}}$ is thus leading and will therefore be reduced by the polynomial division algorithm, which produces remainders of lower degree. In summary, the algebraic dependence relations are accounted for in the constraint ideal, and will be removed by Gröbner basis division. If, for instance, one chooses a lexicographic order (with the $Q_{j}$’s sorted before the $\vec{s}$ variables), this will decrease the power of the $Q_{j}$ for the lowest $j$, at the price of potentially increasing the power of $Q_{j}$ for higher $j$. This argument can be applied recursively until all algebraic dependencies have been resolved. Because the resolution of the dependencies is based on the polynomial ordering and implemented by the Gröbner basis division, it will be consistently implemented across all denominator factors of a given $r_{i}^{\pm}$ and give a unique Leĭnartas decomposition.

We close by noting that this procedure is both simple to implement as well as very efficient, easily handling the complex rational functions present in the problem. Furthermore, an important consequence of this decomposition is that the rational numbers in the expression are more likely to be small. As discussed at the beginning of this section, this makes the result in a single finite field of cardinality $\mathcal{O}(2^{31})$ sufficient to determine them.

5 Results

In this section we take a closer look at the analytic results we derived for a basis set of two-loop five-parton amplitudes. We list the checks we have performed and describe the format of the ancillary files which are used to distribute the results. We also comment on the structure of the remainders.

For each five-parton process we choose a basis set of helicity amplitudes, such that all other choices of helicities and color structures from eqs. 3, 4 and 5 can be obtained by a combination of parity transformation, charge conjugation, and permutations of external momenta. Since we reconstruct two-loop remainders, in order to assemble two-loop amplitudes we also need expressions for the one-loop amplitudes through order $\epsilon^{2}$, see eq. (10). Hence, we have computed analytic expressions for all the corresponding leading-color one-loop amplitudes to all orders in $\epsilon$ in the HV scheme (extending corresponding results of Bern:1993mq ; Kunszt:1994nq ; Bern:1994fz ) and checked that they reproduce the numerical one-loop results presented in Abreu:2018jgq up to $\mathcal{O}(\epsilon^{2})$. The calculation of these amplitudes was also performed by analytically reconstructing the expressions from numerical data.

In order to validate the expressions of the finite two-loop remainders, we have performed a number of checks. First, to check the correctness of the reconstruction procedure, we have compared our analytic results to an independent numerical evaluation of the amplitude on rational phase-space points and found agreement. Second, we have used the analytic results for one-loop amplitudes and two-loop remainders to reproduce all available numerical targets Badger:2013gxa ; Badger:2017jhb ; Gehrmann:2015bfy ; Abreu:2017hqn ; Badger:2018gip ; Abreu:2018jgq ; Dunbar:2016aux . Finally, the analytic results presented here for five-gluon amplitudes at $N_{f}^{0}$ agree with the previously computed expressions Gehrmann:2015bfy ; Dunbar:2016aux ; Badger:2018enw ; Abreu:2018zmy .

We present our results in the form of ancillary files which can be downloaded from the source files of this work on the arXiv server. We include the following files, in Mathematica-readable format:

• •

Analytic expressions for the basis set of one-loop five-parton amplitudes expressed in terms of one-loop master integrals, valid to all orders in $\epsilon$;

• •

Expressions for all one-loop master integrals written in terms of pentagon functions through order $\epsilon^{2}$;

• •

The basis set of two-loop five-parton finite remainders, which are the main result of this work. For each remainder and at each power of $N_{f}$, we give a list of independent coefficient functions $r_{i}$, a matrix $M_{ij}$ and a list of pentagon functions $h_{j}$, consistent with the decomposition of eq. (41);

• •

Scripts that assemble two-loop amplitudes from the remainders. In particular, the Mathematica script TwoLoopAmplitudesNumerical.m uses the analytic expressions to reproduce the numeric values of refs. Abreu:2018jgq ; Badger:2018gip , as detailed in the included README.md file.

Our expressions can be easily adapted to perform numerical integration over phase space. Indeed, they are very compact (with a compressed size of 1.6 Mb), and ready for automated algorithms for optimized evaluation like those included in FORM Kuipers:2012rf ; Ruijl:2014spa . It is interesting to note that whilst there are around 400 independent pentagon functions, for all amplitudes the number of independent rational structures is much lower. A further simplification of the coefficients can be obtained by combining different powers of $N_{f}$ in a decomposition similar to that of eq. (41). Indeed, in figure 9 we observe an overlap of the spaces of independent coefficients functions between distinct ${N_{f}}$ contributions for a given helicity assignment (these are expected from the cancellations present in supersymmetric amplitudes). For the most complicated amplitude, $(g^{-},g^{+},g^{-},g^{+},g^{+})$, the dimension of the combined set of coefficient functions is only 95. We have not presented the amplitudes in this way in the expressions we provide in order to give easier access to the different powers of $N_{f}$. While our expressions are specialized to the Euclidean phase-space region, our coefficients can be used to extract the required information to cover all regions of phase space with minimal work. We leave this to a future publication, where we will explore the numerical evaluation of the amplitudes in more detail.

We end this section by summarizing the computational resources used to obtain our results. The evaluation of the multivariate polynomials $n_{i}^{\pm}(\vec{s})$ in eq. (37), which is performed in a single finite field for all processes, is the most computationally intense. Every other step of the computation, including the extraction of denominator factors, the construction of the matrix $M_{ij}$, 444We note that, for all processes considered, the matrix $M_{ij}$ was computed with only the information contained in the $n_{i}^{\pm}(\vec{s})$, except for the $N_{f}^{0}$ remainder of the processes $(q^{+},{\bar{q}}^{-},Q^{-},{\bar{Q}}^{+},g^{+})$ and $(q^{+},{\bar{q}}^{-},g^{+},g^{-},g^{+})$ where, in order to rationally reconstruct its entries, the numerical computation of the remainder over 100 extra phase-space points in a second finite field was required. performing the partial fractioning of the expressions and rationally reconstructing the numerical coefficients obtained in the previous step, can be performed quickly on a modern laptop computer. The most demanding remainder reconstruction was that of the $N_{f}^{1}$ contribution to the $(g^{-},g^{+},g^{-},g^{+},g^{+})$ process. In total, 94,696 phase-space points (coming in parity conjugate pairs, as discussed in section 4.1) were necessary for its reconstruction. On average (considering all contributions to the reconstruction procedure) the evaluation time per phase-space point amounts to 4.5 minutes.555This timing information is measured on an Intel Xeon E5-2670v3 CPU while running the maximum amount of threads it allows. The total computational resources required to evaluate all our results are relatively modest, and can easily be obtained on a midsize computer cluster.

6 Conclusion

We have presented the leading-color five-parton two-loop scattering amplitudes in analytic form for the first time. These are provided in a set of ancillary files, where we give compact analytic expressions for five-gluon amplitudes as well as amplitudes with two and four external quarks, including in all cases the contributions of closed light-quark loops. These results have been obtained employing a functional-reconstruction approach Peraro:2016wsq to promote numerical unitarity Ita:2015tya ; Abreu:2017xsl ; Abreu:2017idw finite-field evaluations of the amplitude to analytical expressions.

This approach has recently been applied for the planar five-gluon two-loop amplitudes in ref. Abreu:2018zmy and here we develop it further in two main directions. First, we have improved the handling of fermions in the numerical approach. We use dimensional reduction to analytically obtain part of the dependence on the dimensional-regularization parameter by introducing scalar particles. This leads to significant efficiency improvements as compared to a numerical dimensional-reconstruction approach Giele:2008ve ; Ellis:2008ir ; Boughezal:2011br as applied earlier in refs. Abreu:2017xsl ; Abreu:2017hqn ; Abreu:2018jgq . Enhancing the dimensional reconstruction by dimensional reduction is not a new idea and has been applied to gluon amplitudes in the original work Giele:2008ve . Here we apply this improvement Anger:2018ove starting from the five-point amplitudes including external fermions Abreu:2018jgq . Second, we have improved the analytic-reconstruction algorithm from numerical samples. We reconstruct directly in terms of Mandelstam variables as opposed to momentum-twistor variables, which requires considerably fewer evaluations. Furthermore, we simplify the reconstruction procedure by focusing on a linearly independent set of rational functions. Moreover, we have developed a multivariate partial-fraction algorithm in order to perform rational reconstruction with input from only a single finite field of cardinality $\mathcal{O}(2^{31})$. As a by-product of these techniques, we find compact analytic forms and an interestingly small set of independent functions. With these methods, a sample with a size comparable to that required for numerical Monte-Carlo integration can be used to produce analytic expressions.

The computational method we present is robust and efficient. Here, we computed compact analytic expressions for QCD amplitudes depending on five kinematic scales, which had long been a bottleneck. We expect that important scattering amplitudes in the Standard Model depending on a higher number of scales are within the reach of our approach. With the recent progress in non-planar master integral computations Abreu:2018aqd ; Chicherin:2018old , it would also be interesting to explore the amplitudes beyond the leading-color approximation.

The five-parton amplitudes are an important ingredient for obtaining precision phenomenology for the three-jet production process at next-to-next-to-leading order QCD at the LHC in the coming years. In particular, a complete set of compact analytic expressions, such as the ones presented here, will be important for an efficient and numerically stable evaluation of the amplitudes. On a shorter time scale, we believe that our results will be valuable for exploring the analytic properties of scattering amplitudes in QCD.

Acknowledgments

We thank L. J. Dixon, D. A. Kosower and M. Zeng for many useful discussions. We thank C. Duhr for the use of his Mathematica package PolyLogTools. The work of S.A. is supported by the Fonds de la Recherche Scientifique–FNRS, Belgium. The work of J.D., F.F.C. and V.S. is supported by the Alexander von Humboldt Foundation, in the framework of the Sofja Kovalevskaja Award 2014, endowed by the German Federal Ministry of Education and Research. The work of B.P. is supported by the French Agence Nationale pour la Recherche, under grant ANR–17–CE31–0001–01. The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 39/963-1 FUGG.

Appendix A Infrared Structure of Two-Loop Five-Parton Amplitudes

In this appendix we give more details on our definition of the remainders we compute in this paper. In section 2.2 we stated that divergences of renormalized two-loop amplitudes obey a universal structure,

[TABLE]

The renormalized amplitudes can be written in terms of the bare amplitudes as

[TABLE]

with

[TABLE]

For amplitudes in the leading-color approximation the operators $\mathbf{I}^{(1)}_{[n]}$ and $\mathbf{I}^{(2)}_{[n]}$ are diagonal in color space. The operator $\mathbf{I}^{(1)}_{[n]}$ is given by

[TABLE]

with $s_{i,j}=(p_{i}+p_{j})^{2}$ and the indices defined cyclically. The index $a_{i}$ denotes a type of particle with momentum $p_{i}$, i.e., $a_{i}\in\{g,q,\bar{q},Q,\bar{Q}\}$. The symbols $\gamma_{a,b}$ are symmetric under the exchange of indices, $\gamma_{a,b}=\gamma_{b,a}$, and given by:

[TABLE]

The operator ${\bf I}^{(2)}_{[n]}$ is

[TABLE]

where

[TABLE]

and ${\bf H}_{[n]}(\epsilon)$ is a diagonal operator at leading color that depends on the number of external quarks and gluons in the process,

[TABLE]

with

[TABLE]

The two-loop bare amplitude ${\mathcal{A}}^{(2)}$ can be obtained from the remainder defined in eq. (11) using

[TABLE]

Appendix B Dimensionally Reduced Feynman Rules

In the table 1 we list the color-ordered Feynman rules for vertices involving the scalar particles introduced in section 3. The $(D_{s}-6)$-dimensional part of these rules can be fully contracted in each Feynman diagram yielding kinematic-independent factors.

Appendix C Rational Phase-Space Parametrization

We give a twistor parametrization Hodges:2009hk used to rationalize the external on-shell momenta $\{p_{i}\}_{i=1,5}$ with $p_{i}^{2}=0$. This parametrization yields a momentum point with the kinematic invariants $(s_{12},s_{23},s_{34},s_{45},s_{51})$ from the input variables $\{s_{23},s_{45},s_{51},x\}$ and $s_{12}=1$ as given in the main text in eq. (32).

To each external momentum $p_{i}$, we associate the spinors $\lambda_{i}$ and $\tilde{\lambda}_{i}$, from which we can compute the associated momenta through

[TABLE]

The spinors are then parametrized through the momentum-twistor matrix,

[TABLE]

Conjugate spinors $\tilde{\lambda}_{i}$ can then be obtained through

[TABLE]

where $\langle i,j\rangle=\det(\{\lambda_{i},\lambda_{j}\})$ is the usual spinor-bracket computed by taking the determinant of the associated sub-matrix of (63). The Mandelstam invariants $s_{ij}$ and ${\rm tr}_{5}$ required in the main text are given in terms of spinors by,

[TABLE]

using $[i,j]=\det(\{\tilde{\lambda}_{j},\tilde{\lambda}_{i}\})$. For further details, see e.g. section 2 of ref. Bourjaily:2010wh .