# Analytic form of planar two-loop five-gluon scattering amplitudes in QCD

S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page

► To cite this version:

S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page. Analytic form of planar two-loop five-gluon scattering amplitudes in QCD. Physical Review Letters, 2019, 122 (8), pp.082002. 10.1103/Phys-RevLett.122.082002. hal-01965343

HAL Id: hal-01965343

https://hal.science/hal-01965343v1

Submitted on 6 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Analytic Form of Planar Two-Loop Five-Gluon Scattering Amplitudes in QCD

S. Abreu,¹ J. Dormans,² F. Febres Cordero,²,³ H. Ita,² and B. Page⁴

¹Center for Cosmology, Particle Physics and Phenomenology (CP3),
Université Catholique de Louvain, 1348 Louvain-La-Neuve, Belgium

²Physikalisches Institut, Albert-Ludwigs-Universität Freiburg, D-79104 Freiburg, Germany

³Physics Department, Florida State University Tallahassee, Florida 32306, USA

⁴Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay, F-91191 Gif-sur-Yvette cedex, France

<figure><img src="image_2.png" /></figure>

(Received 17 December 2018; published 26 February 2019)

We present the analytic form of the planar two-loop five-gluon scattering amplitudes in QCD for a complete set of independent helicity configurations of external gluons. These include the first analytic results for five-point two-loop amplitudes relevant for the computation of next-to-next-to-leading-order QCD corrections at hadron colliders. The results were obtained by reconstructing analytic expressions from numerical evaluations. The complexity of the computation is reduced by exploiting physical and analytical properties of the amplitudes, employing a minimal basis of so-called pentagon functions that have recently been classified.

DOI: 10.1103/PhysRevLett.122.082002

Scattering amplitudes in QCD are beautiful mathemati-
cal objects that condense the physics of particle collisions.
They play a central role in providing theoretical predictions
for experimental observables, such as those measured in
particle collisions at the Large Hadron Collider at CERN,
making higher-order quantum corrections highly desirable.
The higher the corrections one considers, the more intricate
the amplitudes become, needing to accommodate richer
unitarity and factorization properties. Nevertheless, the
analytic properties are so constraining that final results
appear very compact in spite of the large computational
effort required to obtain them. In this Letter, we exploit
these properties to directly access the compact analytic
results from numerical evaluations.

Analytic expressions for two-loop four-gluon amplitudes
have been known for more than a decade [1,2]. At higher
multiplicity, analytic results have been obtained for all-plus
helicity amplitudes up to seven points [3–8] and up to five
points for the single-minus helicity configuration [9],
reaching a complexity where traditional methods start to
break down. Conventional approaches to analytic calcu-
lations have seen considerable progress [10,11] but involve
large intermediate expressions, the size of which is highly
sensitive to details of the approach and obscures the
simplicity of the result. To sidestep this issue and make
the most of recent progress on five-point two-loop integral

calculations [12–15], two groups have used numerical
approaches to compute five-parton two-loop QCD scatter-
ing amplitudes [16–19]. Numerical approaches are more
resilient to the swelling of intermediate expressions, and
can be improved with finite-field calculations to provide
exact results [20]. While numerical results are sufficient for
Monte Carlo integration, analytic results increase the
evaluation efficiency and facilitate the study of the math-
ematical structure of amplitudes. It has been pointed out
that the analytic form of rational functions in scattering
amplitudes can be reconstructed from numerical samples
[20,21] as already demonstrated for the two-loop four-point
amplitudes [22] and recently for the single-minus five-
gluon amplitude [9]. Combined with the better under-
standing of the analytic properties of the functions appear-
ing in loop amplitudes that was developed over the last few
years, the reconstruction of analytic expressions from
numerical samples circumvents the prohibitively large
intermediate analytic steps and directly targets the simpler
final results.

In this Letter, we use these techniques to compute a
complete set of analytic two-loop five-gluon amplitudes in
the limit of a large number of colors N_c. We apply the
extension of the numerical variant [23–26] of the unitarity
method [27–29] to multiloop computations [22,30,31] to
obtain numerical samples on finite fields based on
Refs. [17,19]. From these samples we obtain analytic
expressions for the coefficients in a decomposition of
the amplitude in terms of so-called pentagon functions
[13]. The efficiency of the reconstruction is increased by
finding convenient bases in the space of pentagon functions
and by an a priori knowledge of the coefficients of the
denominators.

Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI. Funded by SCOAP³.
We present for the first time analytic results for the two-loop five-gluon amplitudes known as MHV amplitudes, which are an important ingredient in the next-to-next-to-leading-order (NNLO) QCD corrections for three-jet production at hadron colliders. The results will serve as a testing ground for exploring the structure of two-loop scattering amplitudes in quantum field theory.

**Scattering amplitudes.**—The main results of this Letter are analytic expressions for two-loop five-gluon scattering amplitudes, in the leading-color approximation of pure Yang-Mills theory. We write the perturbative expansion of a bare five-gluon helicity amplitude as

$$ \mathcal{A}(\{p_i, h_i\}_{i=1,...,5})|_{\text{leading color}} = \sum_{\sigma \in S_5/Z_5} \text{Tr}[T^{a_{\sigma(1)}} T^{a_{\sigma(2)}} T^{a_{\sigma(3)}} T^{a_{\sigma(4)}} T^{a_{\sigma(5)}}] \times g_0^3(\mathcal{A}^{(0)} + \lambda \mathcal{A}^{(1)} + \lambda^2 \mathcal{A}^{(2)} + \mathcal{O}(\lambda^3)), \quad (1) $$

where $\lambda = N_c g_0^2/(4\pi)^2$ and $g_0$ is the bare QCD coupling. The set $S_5/Z_5$ denotes all noncyclic permutations of five indices. The planar amplitudes $\mathcal{A}^{(k)}$ are functions of the momenta $p_{\sigma(i)}$ and helicities $h_{\sigma(i)}$. In the leading-color approximation there is a single color structure and it is sufficient to specify a helicity assignment for the ordered set of legs, which we specify by subscripts

$$ \mathcal{A}_{h_1 h_2 h_3 h_4 h_5}^{(k)} $$

when required, assuming an ordered set of momentum assignments $\{p_1, p_2, p_3, p_4, p_5\}$. In this section we describe our approach to the calculation of the $\mathcal{A}^{(k)}$ in dimensional regularization with $D = 4 - 2\epsilon$. Although we focus on the leading-color contributions, the approach we use is fully generic and applicable beyond this approximation.

Our calculation of the amplitudes $\mathcal{A}^{(k)}$ is performed in the framework of two-loop numerical unitarity [17,22,30,31]. The starting point is the general parametrization of the integrand of the amplitude [30]

$$ \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}, \quad (2) $$

with $M_\Gamma$ being a set of master integrands, $S_\Gamma$ being a set of surface terms, $P_\Gamma$ as the set of propagators $\rho_j$ being associated with propagator structure $\Gamma$, and $\ell_l$ as the set of loop momenta. The coefficients $c_{\Gamma, i}$ are determined by exploiting the factorization properties of $\mathcal{A}(\ell_l)$ in loop-momenta limits $\ell_l \to \ell_l^\Gamma$ where propagators go on-shell [32]:

$$ \sum_{\text{states}} \prod_{i \in T_{\Gamma}} \mathcal{A}_{i}^{\text{tree}}(\ell_{l}^{\Gamma}) = \sum_{\substack{\Gamma' \ge \Gamma_{\Gamma} \\ i \in M_{\Gamma'} \cup S_{\Gamma'}}} \frac{c_{\Gamma', i} m_{\Gamma', i}(\ell_{l}^{\Gamma})}{\prod_{j \in (P_{\Gamma'} \setminus P_{\Gamma})} \rho_{j}(\ell_{l}^{\Gamma})}. \quad (3) $$

The set $T_\Gamma$ labels all tree amplitudes corresponding to the vertices in the diagram $\Gamma$. The sum runs over the $\Gamma'$ with more propagators. The state sum runs over the (scheme dependent) gluon-helicity states of each internal line of $\Gamma$. Solving the linear system in Eq. (3) directly yields the coefficients of master integrals and no further integral reduction is needed. In numerical unitarity, the linear system [Eq. (3)] is constructed and solved numerically. Using finite-field arithmetic removes any issues related to loss of precision in numerical operations, at the price of minor modifications to the standard numerical unitarity approach. We describe our algorithm in Ref. [17].

**Numerical reconstruction.**—We aim to compute analytic expressions for five-gluon amplitudes by reconstructing them from numerical samples. For a suitably chosen parametrization of phase space, such as momentum twistors [33], the coefficients $c_{\Gamma, i}$ of Eq. (2) are rational functions. For concreteness, we use the parametrization [21]

$$ \begin{aligned} s_{12} &= x_4, & s_{23} &= x_2 x_4, & s_{34} &= x_4 \left( \frac{(1+x_1)x_2}{x_0} + x_1(x_3-1) \right), \\ s_{45} &= x_3 x_4, & s_{51} &= x_1 x_4 (x_0 - x_2 + x_3), \\ \text{tr}_5 &= 4i\epsilon(p_1, p_2, p_3, p_4) \\ &= x_4^2 \left( x_2(1+2x_1) + x_0 x_1 (x_3-1) - \frac{x_2(1+x_1)(x_2-x_3)}{x_0} \right), \end{aligned} \quad (4) $$

where $s_{ij} = (p_i + p_j)^2$ with the indices defined cyclically. The difficulty of analytic reconstruction is determined by the complexity of the functions under study. More precisely, this means that one should attempt to reconstruct rational functions with low total degree of numerator and denominator polynomials. With this aim in mind, we first exploit a series of known physical and analytical properties of the amplitudes, and only apply the reconstruction algorithms to simpler objects that we cannot further constrain. In this section we summarize our approach.

The divergence structure of scattering amplitudes is governed by known universal functions [34–37]. It is thus sufficient to compute the so-called finite remainder [34], defined by subtracting contributions that are related to tree and one-loop amplitudes from a two-loop amplitude. There is freedom in how the remainders are defined, so we now give our definitions. For helicity amplitudes which vanish at tree-level, $\mathcal{A}_{\pm++++}^{(k)}$, we use

$$ \mathcal{R}_{\pm++++}^{(2)} = \bar{\mathcal{A}}_{\pm++++}^{(2)} + S_{\epsilon} \bar{\mathcal{A}}_{\pm++++}^{(1)} \sum_{i=1}^{5} \frac{(-s_{i,i+1})^{-\epsilon}}{\epsilon^2} + \mathcal{O}(\epsilon), \quad (5) $$
where $S_\epsilon = (4\pi)^\epsilon e^{-\epsilon\gamma_\epsilon}$. The $\bar{\mathcal{A}}^{(k)}$ denote amplitudes normalized to remove any ambiguity related to overall phases. In the case of amplitudes that vanish at tree level, we normalize to the leading order in $\epsilon$ of the (finite) one-loop amplitude. For the MHV amplitudes, $\mathcal{A}_{-\mp\pm\pm}^{(k)}$, which we normalize to the corresponding tree amplitude, we define

$$ \mathcal{R}_{-\mp\pm\pm}^{(2)} = \bar{\mathcal{A}}_{-\mp\pm\pm}^{(2)} - \left(\frac{5\tilde{\beta}_0}{2\epsilon} + \mathbf{I}^{(1)}\right) S_\epsilon \bar{\mathcal{A}}_{-\mp\pm\pm}^{(1)} \\ + \left(\frac{15\tilde{\beta}_0^2}{8\epsilon^2} + \frac{3}{2\epsilon}(\tilde{\beta}_0\mathbf{I}^{(1)} - \tilde{\beta}_1) - \mathbf{I}^{(2)}\right) S_\epsilon^2 + \mathcal{O}(\epsilon), \qquad (6) $$

where $\tilde{\beta}_i$ are the coefficients of the QCD $\beta$-function divided by $N_c^{i+1}$. $\mathbf{I}^{(1)}$ and $\mathbf{I}^{(2)}$ are the standard Catani operators at leading color, and Eq. (6) is obtained by writing Catani's prediction for the poles of a bare amplitude. Precise expressions for the operators in our conventions can be found in Appendix B of Ref. [19]. We note that in both Eqs. (5) and (6) we require one-loop amplitudes expanded up to order $\epsilon^2$.

The amplitudes we are concerned with in this Letter can be written as linear combinations of special functions, the so-called multiple polylogarithms (MPLs). The coefficients in this linear combination are functions of the external data and of the dimensional regulator $\epsilon$, just as in the decomposition of Eq. (2). It is by now well understood that there are relations between different MPLs, and newly developed mathematical tools allow us to find such relations in an algorithmic way [38–40]. It is then possible to construct a basis for the space of special functions required for planar five-point massless amplitudes up to two-loops, and this was done in Ref. [13] where a set of pentagon functions was defined. Amplitudes can then be written as

$$ \bar{\mathcal{A}}^{(2)} = \sum_{i \in B} \sum_{k=-4}^{0} \epsilon^k \tilde{c}_{k,i}(\vec{x}) h_i(\vec{x}) + \mathcal{O}(\epsilon), \qquad (7) $$

where $\vec{x} = \{x_0, x_1, x_2, x_3, x_4\}$ and the $\tilde{c}_{k,i}(\vec{x})$ are rational functions of the twistor variables. B denotes the basis of pentagon functions $h_i$. There is a notion of (transcendental) weight associated with MPLs, and thus with pentagon functions, that allows us to organize the basis B. For instance, there is a single element of weight 0, the trivial function 1. At weight one, there are five elements, $\ln(-s_{i,i+1})$. At weight two, there are products of weight-one functions, irreducible weight-two functions, and a new constant, $\pi^2$. The same pattern holds at weight three and four, the highest weight required for two-loop amplitudes at order $\epsilon^0$. To make the simplifications that one expects to find in the remainders explicit, we write both one- and two-loop amplitudes in terms of pentagon functions so that the remainders themselves are written as combinations of pentagon functions:

$$ \mathcal{R}^{(2)} = \sum_{i \in B} r_i(\vec{x}) h_i(\vec{x}). \qquad (8) $$

The coefficients $r_i(\vec{x})$ have lower total degrees compared with the $\tilde{c}_{k,i}(\vec{x})$ of Eq. (7), but to further increase the efficiency of the reconstruction we find it useful to implement a series of improvements by investigating the structure of the coefficients on generic “univariate slices.” Such slices are defined by parametrizing the twistor variables $\vec{x}$ by a single variable $t$, $x_i = a_i + b_i t$, with the $a_i$ and $b_i$ chosen randomly in the finite field (for high enough cardinality, this ensures the $x_i$ do not satisfy simple relations leading to artificial simplifications). The coefficients $r_i(\vec{x})$ are themselves univariate rational functions of $t$, $r_i[\vec{x}(t)] = r_i(t)$. Importantly, on such univariate slices the degrees of the numerator and denominator of the rational functions $r_i(t)$ correspond to the total degrees of $r_i(\vec{x})$ in the $x_j$. We use these slices to probe the complexity of the functions and find simplifications. First, we classify the pole structures of the coefficient functions $r_i(\vec{x})$. A similar classification has been exploited for the computation of one-loop QCD amplitudes in the past [26]. On physical grounds we expect the pole structure to be determined by the so-called alphabet of the pentagon functions. The alphabet determines the points in phase space where the pentagon functions (or their discontinuities) have logarithmic singularities, and they provide a natural ansatz for the poles of the coefficients. The five-point planar alphabet can be written in terms of 26 letters $W_i$, see Ref. [13], which we rewrite in terms of twistor variables. This gives a set A of 26 independent letters $A = \{w_i(\vec{x})\}$. Indeed, we find that the denominators do factorize into products of letters,

$$ r_i(\vec{x}) = \frac{n_i(\vec{x})}{\prod_{j \in A} w_j(\vec{x})^{q_{ij}}}. \qquad (9) $$

The integer exponents $q_{ij}$ are determined by matching this ansatz on univariate slices. The computation of the $r_i(\vec{x})$ is then reduced to the much simpler multivariate polynomial reconstruction of the numerators $n_i(\vec{x})$. Second, it is expected that cancellations between different basis functions in B take place in exceptional kinematic configurations, which implies relations between their coefficients. This motivates a search for a different basis of pentagon functions with coefficients of a lower total degree. To find this new (helicity-dependent) basis, we construct linear combinations of coefficients

$$ \sum_{i \in B} a_{i,k} r_i(\vec{x}) = \frac{N_k(\vec{x}, a_{i,k})}{\prod_{j \in A} w_j(\vec{x})^{q_{kj}}}, \qquad (10) $$

and solve for phase space independent $a_{i,k}$ such that the numerators $N_k(\vec{x}, a_{i,k})$ factorize a subset of the $w_j \in A$. This can be performed on univariate slices by only accepting solutions which are invariant over multiple slices. The
<table><caption>TABLE I. Each $t^n/t^d$ denotes the total degree of numerator ($n$) and denominator ($d$) of the most complex coefficient for each helicity amplitude in the decomposition of Eq. (7) (second column) or Eq. (8) (third column). The fourth column lists the highest polynomial we reconstruct. The final column lists the number of letters $w_j(\vec{x})$ that contribute in the denominator of Eq. (9).</caption><thead><tr><th>Helicity</th><th>$\tilde{c}_{k,i}(t)$</th><th>$r_i(t)$</th><th>$n'_i(t)$</th><th>$w_j$s in denominator</th></tr></thead><tbody><tr><td>+++++</td><td>$t^{34}/t^{28}$</td><td>$t^{10}/t^4$</td><td>$t^{10}$</td><td>3</td></tr><tr><td>-++++</td><td>$t^{50}/t^{42}$</td><td>$t^{35}/t^{28}$</td><td>$t^{35}$</td><td>14</td></tr><tr><td>--+++</td><td>$t^{70}/t^{65}$</td><td>$t^{50}/t^{45}$</td><td>$t^{40}$</td><td>17</td></tr><tr><td>---+</td><td>$t^{84}/t^{82}$</td><td>$t^{68}/t^{66}$</td><td>$t^{53}$</td><td>20</td></tr></tbody></table>

matrix $a_{i,k}$ allows us to change to a new basis $B'$ in the space of special functions, in which remainders can be decomposed as in Eq. (8), with coefficients $r'_i(\vec{x})$ whose numerators $n'_i(\vec{x})$ are polynomials of lower total degree than those of Eq. (9).

We believe that a deeper understanding of these two improvements might be fruitfully exploited in the future to tackle more complex computations. In particular, it would be interesting to better understand how the branch structure of pentagon functions is related to the poles of the coefficients, and if this suggests a superior basis of pentagon functions with simpler coefficients.

**Implementation.**—Let us now discuss how we applied the techniques of the previous section to perform the analytic reconstruction. We first compute numerical values for the coefficients in a decomposition of one- and two-loop amplitudes in terms of master integrals, using our C++ implementation of numerical unitarity on a finite field [17,19]. We then introduce expressions for the master integrals in terms of pentagon functions and compute the remainders in Eq. (5) or Eq. (6). After changing the basis of pentagon functions and multiplying by the predetermined denominator structures, we arrive at an evaluation of the numerator polynomials $n'_i(\vec{x})$. Table I summarizes the impact of the improvements, the difficulty in the reconstruction of each amplitude, and the number of different letters appearing in the denominators of each helicity. Given the simplicity of $R_{\pm++++}$, we did not implement any basis change for these helicities. We observe a drastic increase in complexity from top to bottom, as expected from one-loop amplitudes [41]. Comparing the second and third columns we find that the coefficients in Eq. (8) are indeed simpler than those in Eq. (7). We also show that not all letters contribute to the denominator in Eq. (9), even for the most complicated helicity.

The reconstruction of the multivariate polynomials is done with an in-house implementation of the algorithm presented in Ref. [21]. For a detailed description we refer to the original article. The algorithm [21] correlates the sampled phase space to the function it reconstructs,

learning about simplifications as it progresses. For this reason, it is not trivial to parallelize. In our implementation this is addressed by anticipating a suitable superset of the points for a set of functions. This implementation is well adapted to evaluations on modern computer clusters.

In order to demonstrate the efficiency of our approach, we now summarize the computational requirements for the numerical reconstruction of the presented amplitudes. The dense nature of the algorithm implies that the complexity grows as the number of terms in a polynomial of total degree $R$ in $n$ variables, $(R+n/n)$. The functions we reconstruct are dimensionless and thus only depend on the four variables $x_0, x_1, x_2$, and $x_3$, see Eq. (4). For example, for a generic polynomial with $R=53$, we would need around 400 000 evaluations. In practice, we observe that the polynomials we reconstruct are not completely generic and we require fewer evaluations. For the most complicated remainder, $R_{-++++}$, we require 237 098 phase-space points for the reconstruction of all the functions $n'_i(x)$. For $R_{-++++}$, we require 85 979 phase-space points. The average time for the numerical computation of the remainder of an MHV amplitude in a finite-field is 4 minutes per phase-space point. The reconstruction of the remainders $R_{\pm++++}$ is simpler. We perform the computation on a finite field of cardinality $O(2^{31})$. For the MHV amplitudes we also evaluate on a second finite field and apply the Chinese remainder theorem. With this, we are able to rationally reconstruct the result.

**Results.**—We have validated the approach described above by computing the one-loop amplitudes to all orders in $\epsilon$ (extending the results of Ref. [41]). We then used it to obtain the analytic form of all planar five-gluon two-loop amplitudes in the leading-color approximation, in the 't Hooft-Veltman regularization scheme. Our results are included in a set of ancillary files [42]. The files contain analytic expressions for the remainders of four helicity configurations—$A_{\pm++++}$, see Eq. (5), and $A_{-\mp++++}$, see Eq. (6)—from which any other helicity configuration can be obtained by symmetry or permutation of the momenta. These remainders are written as a linear combination of pentagon functions, similar to Eq. (8) but with the adjusted bases. This allows us to use the library provided in Ref. [13] to evaluate the remainders. We also include analytic expressions for the one-loop amplitudes, which are required to obtain the two-loop amplitudes from their remainder. These are written as a linear combination of master integrals in the one-loop basis of Ref. [25] and valid to all orders in $\epsilon$. We provide expressions for the expansion of the one-loop master integrals through order $\epsilon^2$, written in terms of pentagon functions, so that they can easily be combined with the expressions for the remainders. Finally, we include a script that assembles all components to evaluate a two-loop amplitude.

The results we present are valid in the Euclidean region. Since they are written in terms of pentagon functions [13]
there is minimal work to extend the results to all kinematic regions. We checked our expressions against recently computed two-loop amplitudes [3,5,7,9,16,17,19]. In the ancillary files we include numerical values for the pentagon functions at the phase-space point of Ref. [19] and a script that numerically evaluates the amplitudes at that point.

**Conclusions.**—We have presented the analytic form of a complete set of planar two-loop five-gluon amplitudes in the leading-color approximation. The amplitudes were obtained by reconstructing the analytic expressions from numerical samples on finite fields, computed in the framework of two-loop numerical unitarity. This allows us to make the most of the resilience of numerical calculations to handle intermediate steps and to only target the analytic form of final expressions, which is constrained by various physical properties. By focusing on finite remainders, we reduced the complexity of the objects to reconstruct. Efficiency of the reconstruction was further enhanced by an *a priori* determination of the denominators and changes of basis that considerably reduce the total degrees of the numerators of the coefficients. Through this process, we reduced the calculation of the coefficients to the reconstruction of a multivariate polynomial of relatively low total degree, rendering the most complex MHV amplitudes easily reconstructible on a modern computer cluster.

Our results include the first computation of the analytic form of the five-gluon MHV amplitudes at two loops. They are an important contribution to the NNLO QCD corrections to three-jet production at hadron colliders. While in principle the numerical evaluation of amplitudes is sufficient, the efficiency requirements for phase-space integration over the final states are high due to helicity and color sums. The analytic expressions that we provide will help to control both the evaluation times and the numerical stability of future phenomenological studies.

We expect our computational approach to greatly contribute to formal developments in the study of scattering amplitudes in quantum field theory, and to the new era of precision QCD in high-energy physics.

We thank V. Sotnikov 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. and F. F. C. 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 No. ANR-17-CE31-0001-01. The authors acknowledge support by the state of Baden-Württemberg through bwHPC. S. A. and B. P. thank the Galileo Galilei Institute for Theoretical Physics for the hospitality and

the INFN for partial support during the final stages of this work.

[1] E. W. N. Glover, C. Oleari, and M. E. Tejeda-Yeomans, Nucl. Phys. **B605**, 467 (2001).

[2] Z. Bern, A. De Freitas, and L. J. Dixon, J. High Energy Phys. 03 (2002) 018.

[3] S. Badger, H. Frellesvig, and Y. Zhang, J. High Energy Phys. 12 (2013) 045.

[4] S. Badger, G. Mogull, A. Ochirov, and D. O'Connell, J. High Energy Phys. 10 (2015) 064.

[5] T. Gehrmann, J. M. Henn, and N. A. Lo Presti, Phys. Rev. Lett. **116**, 062001 (2016); **116**, 189903(E) (2016).

[6] S. Badger, G. Mogull, and T. Peraro, J. High Energy Phys. 08 (2016) 063.

[7] D. C. Dunbar and W. B. Perkins, Phys. Rev. D **93**, 085029 (2016).

[8] D. C. Dunbar and W. B. Perkins, Phys. Rev. Lett. **117**, 061602 (2016).

[9] S. Badger, C. Brønnum-Hansen, H. B. Hartanto, and T. Peraro, J. High Energy Phys. 01 (2019) 186.

[10] R. H. Boels, Q. Jin, and H. Luo, arXiv:1802.06761.

[11] H. A. Chawdhry, M. A. Lim, and A. Mitov, arXiv: 1805.09182.

[12] C. G. Papadopoulos, D. Tommasini, and C. Wever, J. High Energy Phys. 04 (2016) 078.

[13] T. Gehrmann, J. M. Henn, and N. A. Lo Presti, J. High Energy Phys. 10 (2018) 103.

[14] S. Abreu, B. Page, and M. Zeng, J. High Energy Phys. 01 (2019) 006.

[15] D. Chicherin, T. Gehrmann, J. M. Henn, N. A. Lo Presti, V. Mitev, and P. Wasser, arXiv:1809.06240.

[16] S. Badger, C. Brønnum-Hansen, H. B. Hartanto, and T. Peraro, Phys. Rev. Lett. **120**, 092001 (2018).

[17] S. Abreu, F. Febres Cordero, H. Ita, B. Page, and M. Zeng, Phys. Rev. D **97**, 116014 (2018).

[18] S. Badger, C. Brønnum-Hansen, T. Gehrmann, H. B. Hartanto, J. Henn, N. A. Lo Presti, and T. Peraro, Proc. Sci., LL2018 (2018) 006.

[19] S. Abreu, F. Febres Cordero, H. Ita, B. Page, and V. Sotnikov, J. High Energy Phys. 11 (2018) 116.

[20] A. von Manteuffel and R. M. Schabinger, Phys. Lett. B **744**, 101 (2015).

[21] T. Peraro, J. High Energy Phys. 12 (2016) 030.

[22] S. Abreu, F. Febres Cordero, H. Ita, M. Jaquier, B. Page, and M. Zeng, Phys. Rev. Lett. **119**, 142001 (2017).

[23] G. Ossola, C. G. Papadopoulos, and R. Pittau, Nucl. Phys. **B763**, 147 (2007).

[24] R. K. Ellis, W. T. Giele, and Z. Kunszt, J. High Energy Phys. 03 (2008) 003.

[25] W. T. Giele, Z. Kunszt, and K. Melnikov, J. High Energy Phys. 04 (2008) 049.

[26] C. F. Berger, Z. Bern, L. J. Dixon, F. Febres Cordero, D. Forde, H. Ita, D. A. Kosower, and D. Maitre, Phys. Rev. D **78**, 036003 (2008).

[27] Z. Bern, L. J. Dixon, D. C. Dunbar, and D. A. Kosower, Nucl. Phys. **B425**, 217 (1994).

[28] Z. Bern, L. J. Dixon, D. C. Dunbar, and D. A. Kosower, Nucl. Phys. **B435**, 59 (1995).
[29] R. Britto, F. Cachazo, and B. Feng, Nucl. Phys. **B725**, 275 (2005).

[30] H. Ita, Phys. Rev. D **94**, 116015 (2016).

[31] S. Abreu, F. Febres Cordero, H. Ita, M. Jaquier, and B. Page, Phys. Rev. D **95**, 096011 (2017).

[32] At two loops, there are subleading poles in this limit, which we handle with the approach described in Ref. [31].

[33] A. Hodges, J. High Energy Phys. 05 (2013) 135.

[34] S. Catani, Phys. Lett. B **427**, 161 (1998).

[35] G. F. Sterman and M. E. Tejeda-Yeomans, Phys. Lett. B **552**, 48 (2003).

[36] T. Becher and M. Neubert, Phys. Rev. Lett. **102**, 162001 (2009); **111**, 199905(E) (2013).

[37] E. Gardi and L. Magnea, J. High Energy Phys. 03 (2009) 079.

[38] A. B. Goncharov, M. Spradlin, C. Vergu, and A. Volovich, Phys. Rev. Lett. **105**, 151605 (2010).

[39] C. Duhr, H. Gangl, and J. R. Rhodes, J. High Energy Phys. 10 (2012) 075.

[40] C. Duhr, J. High Energy Phys. 08 (2012) 043.

[41] Z. Bern, L. J. Dixon, and D. A. Kosower, Phys. Rev. Lett. **70**, 2677 (1993).

[42] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.122.082002, for the analytic expressions.