Subleading Poles in the Numerical Unitarity Method at Two Loops

TL;DR

This paper introduces a universal algorithm to compute subleading pole contributions in two-loop scattering amplitudes, enhancing the numerical unitarity method by enabling direct calculation of complex integral coefficients.

Contribution

The authors develop a novel algorithm to extract subleading pole terms at two loops, addressing a gap in the numerical unitarity approach for scattering amplitudes.

Findings

Successfully computed two-loop four-gluon integral coefficients numerically.

Demonstrated the effectiveness of the new algorithm in handling subleading poles.

Enhanced the numerical unitarity method for two-loop amplitude calculations.

Abstract

We describe the unitarity approach for the numerical computation of two-loop integral coefficients of scattering amplitudes. It is well known that the leading propagator singularities of an amplitude's integrand are related to products of tree amplitudes. At two loops, Feynman diagrams with doubled propagators appear naturally, which lead to subleading pole contributions. In general, it is not known how these contributions can be directly expressed in terms of a product of on-shell tree amplitudes. We present a universal algorithm to extract these subleading pole terms by releasing some of the on-shell conditions. We demonstrate the new approach by numerically computing two-loop four-gluon integral coefficients.

Equations

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

$$\mathcal{A}({\ell }_l)=\sum _{\Gamma \in \Delta }\frac{1}{{\prod }_{k\in P_{\Gamma }}{\rho }_k}\sum _{i\in M_{\Gamma }\cup S_{\Gamma }}{c}_{\Gamma ,i}{m}_{\Gamma ,i}({\ell }_l),$$

$$\int\frac{d^{D}\ell_{1}d^{D}\ell_{2}}{(2\pi)^{2D}}\,\frac{m_{\Gamma,i}(\ell_{l})}{\prod_{k\in P_{\Gamma}}\rho_{k}}=\left\{\begin{array}[]{cc}I_{\Gamma,i}&\mbox{for}\quad i\in M_{\Gamma}\,,\\ 0&\mbox{for}\quad i\in S_{\Gamma}\,.\end{array}\right.$$

$$N(\Gamma ,{\ell }_l)$$

$${\ell }_l^{\Gamma }:{\ell }_l\text{mbox}with{\rho }_k=0\text{mbox}forallk\in P_{\Gamma }.$$

$${\Delta }^{\prime }\subseteq \Delta .$$

$$\underset{{\ell }_l\to {\ell }_l^{\Gamma }}{\lim }\mathcal{A}({\ell }_l)$$

$$R(\Gamma ,{\ell }_l^{\Gamma })$$

$$N(\Gamma ,{\ell }_l^{\Gamma })=R(\Gamma ,{\ell }_l^{\Gamma }).$$

$$N({\Gamma }_{\text{NM}},{\ell }_l^{{\Gamma }_{\text{NM}}})=R({\Gamma }_{\text{NM}},{\ell }_l^{{\Gamma }_{\text{NM}}})-\sum _k\frac{1}{{\rho }_k({\ell }_l^{{\Gamma }_{\text{NM}}})}N({\Gamma }_{\text{NM},k},{\ell }_l^{{\Gamma }_{\text{NM}}}),$$

$$\underset{{\ell }_l\to {\ell }_l^{{\Gamma }_d}}{\lim }\mathcal{A}({\ell }_l)$$

$$R({\Gamma }^{\prime },{\ell }_l^{{\Gamma }^{\prime }})=N({\Gamma }^{\prime },{\ell }_l^{{\Gamma }^{\prime }})+\sum _{\begin{array}{c}\Gamma \in \Delta \\ \Gamma >{\Gamma }^{\prime }\end{array}}\frac{N(\Gamma ,{\ell }_l^{{\Gamma }^{\prime }})}{{\prod }_{k\in P_{\Gamma }\setminus P_{{\Gamma }^{\prime }}}{\rho }_k({\ell }_l^{{\Gamma }^{\prime }})},$$

$$N({\Gamma }^{\prime },{\ell }_l^{{\Gamma }^{\prime }})+\sum _{\begin{array}{c}\Gamma \in \tilde{\Delta }\\ \Gamma >{\Gamma }^{\prime }\end{array}}\frac{N(\Gamma ,{\ell }_l^{{\Gamma }^{\prime }})}{{\prod }_{k\in P_{\Gamma }\setminus P_{{\Gamma }^{\prime }}}{\rho }_k({\ell }_l^{{\Gamma }^{\prime }})}=R({\Gamma }^{\prime },{\ell }_l^{{\Gamma }^{\prime }})-\sum _{\begin{array}{c}\Gamma \in \Delta \setminus \tilde{\Delta }\\ \Gamma >{\Gamma }^{\prime }\end{array}}\frac{N(\Gamma ,{\ell }_l^{{\Gamma }^{\prime }})}{{\prod }_{k\in P_{\Gamma }\setminus P_{{\Gamma }^{\prime }}}{\rho }_k({\ell }_l^{{\Gamma }^{\prime }})}.$$

$$N\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/Triangle}}},\ell^{\rm tri}\right)=R\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/Triangle}}},\ell^{\rm tri}\right)-\frac{1}{\rho}N\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/Box}}},\ell^{\rm tri}\right)\ ,$$

$$\sum_{i=1}^{m}c_{{\rm tri},i}m_{{\rm tri},i}(\ell^{\rm tri})+\frac{1}{\rho}\sum_{i=1}^{n}c_{{\rm box},i}\,m_{{\rm box},i}(\ell^{\rm tri})=R\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/Triangle}}},\ell^{\rm tri}\right)\,,$$

$$\displaystyle N\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/DoubleBoxSmall.pdf}}},\ell_{l}^{\rm c}\right)$$

$$\displaystyle N\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/BoxTriangle1Small.pdf}}},\ell_{l}^{\rm f}\right)=R\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/BoxTriangle1Small.pdf}}},\ell_{l}^{\rm f}\right)-\frac{1}{\rho_{\rm fb}}N\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/PentaTriangle1Small.pdf}}},\ell_{1}^{\rm f}\right)-\frac{1}{\rho_{\rm fc}}N\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/DoubleBoxSmall.pdf}}},\ell_{l}^{\rm f}\right)\ ,$$

$$\displaystyle N\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/BubbleBoxGeneric.pdf}}},\ell_{l}^{\rm h}\right)+$$

$$\displaystyle\quad R\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/BubbleBoxGeneric.pdf}}},\ell_{l}^{\rm h}\right)-\frac{1}{\rho_{\rm hf}}{N}\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/BoxTriangle1Small.pdf}}},\ell_{l}^{\rm h}\right)-\frac{1}{\rho_{\rm hg}}{N}\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/BoxTriangle2Small.pdf}}},\ell_{l}^{\rm h}\right)-\frac{1}{(\rho_{\rm he})^{2}}{N}\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/HexaBubbleSmall}}},\ell_{l}^{\rm h}\right)$$

$$\displaystyle\quad-\frac{1}{\rho_{\rm hf}\rho_{\rm fb}}{N}\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/PentaTriangle1Small.pdf}}},\ell_{l}^{\rm h}\right)-\frac{1}{\rho_{\rm hf}\rho_{\rm fc}}{N}\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/DoubleBoxSmall.pdf}}},\ell_{l}^{\rm h}\right)-\frac{1}{\rho_{\rm hg}\rho_{\rm gd}}{N}\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/PentaTriangle2Small.pdf}}},\ell_{l}^{\rm h}\right)\ ,$$

Full text

I Introduction

The unitarity method Unitarity constructs scattering amplitudes from their unitarity and analytic structure. It is convenient to work at the integrand level where factorization properties tie the leading coefficients of the propagator poles to products of tree amplitudes. In the last decade, numerical approaches OPP ; NumUnitarity ; BlackHat ; GKM have been developed that construct one-loop amplitudes from their propagator poles, which are given by tree amplitudes. Beyond one loop, Feynman rules naturally yield contributions with higher-order propagator powers. In order to determine the rational integrand, one thus has to obtain leading and subleading coefficients on such propagator poles. Of these, only the leading ones are directly related to a product of tree amplitudes and the subleading terms have to be obtained differently. For analytic computations this obstruction has been discussed in DPropAnalytic where the residue extraction is adjusted to pick up subleading-pole contributions. Alternatively, in analytic computations of two-loop QCD amplitudes fourgluons , the subleading-pole contributions can be tracked explicitly and evaluated or, in the case of particular helicity amplitudes, dealt with by choosing particular representations of the integrand UnitarityAnalytics2loop .

In this article we propose a numerical algorithm to extract subleading-pole contributions without resorting to analytic manipulations. The central idea is to ‘cut less’, i.e. to obtain subleading-pole contributions from their contribution to cuts which keep the respective propagators off-shell. We find that the algorithm works effectively when applied to a numerical calculation of two-loop four-gluon amplitudes. We validate our approach by comparing it to the known analytic results fourgluons , and by carrying out a number of non-trivial consistency checks.

The rest of the paper is organized as follows. In Section II we describe the organization of a calculation in the numerical unitarity method, discuss the appearance of subleading-pole terms starting at two-loops, and present our algorithm for extracting those terms. Section III contains applications in the context of one- and two-loop four-gluon amplitudes as well as a list of checks performed. Finally, in Section IV we present our conclusions and outlook.

II Setup for two-loop numerical unitarity

In this section, we review the main aspects of the numerical unitarity method at two loops. We introduce our notation, explain the appearance of subleading-pole terms at two loops and present our algorithm to extract them. Although we focus on two-loop calculations, our result can be easily generalized for computing general subleading-pole contributions in multi-loop amplitudes.

II.1 Diagrammatic decomposition

We will organize amplitudes in terms of diagrams, which can be constructed from the usual Feynman diagram decomposition in the following way: we strip Feynman diagrams of particle information and pinch all propagators through which loop momentum does not flow. We will denote the set of all diagrams constructed in this way by $\Delta$. In the remainder of this paper, a diagram $\Gamma$ denotes an element of this set, $\Gamma\in\Delta$. Each $\Gamma$ defines a set of propagator indices $P_{\Gamma}$, and we call the set of inverse propagators $\{\rho_{k}\}$ with $k\in P_{\Gamma}$ a propagator structure. As we allow for elements in $P_{\Gamma}$ to be repeated, technically we should call $P_{\Gamma}$ a multiset, but we avoid this terminology. In fig. 1 we show a generic diagram $\Gamma$ for a planar two-loop amplitude. Finally, we associate the notion of hierarchy to $\Delta$. If $\Gamma_{1}$ and $\Gamma_{2}$ are two elements of $\Delta$, such that $\Gamma_{2}$ is obtained by pinching some of the edges of $\Gamma_{1}$, then they are members of the same hierarchy. $\Gamma_{1}$ is called an ancestor of $\Gamma_{2}$, and we write $\Gamma_{1}>\Gamma_{2}$ to denote ancestry. All of the descedant’s propagators are contained in any of its ancestors, i.e. $P_{\Gamma_{2}}\subset P_{\Gamma_{1}}$. Different hierarchies in $\Delta$ are named according to their element with the fewest edges.

II.2 Master integrals and integrands

The scattering amplitude ${\cal A}$ is decomposed in the general form

[TABLE]

in terms of a set of master integrals $I_{\Gamma,i}$ and coefficient functions $c_{\Gamma,i}$. We organize the sum over master integrals according to their propagator structure labeled by $\Gamma$. By $M_{\Gamma}$ we denote the set of indices $i$ which label the different master integrals that share the same propagator structure. $M_{\Gamma}$ can be empty if no master integral is associated with diagram $\Gamma$. We work in dimensional regularization, so that the integrals and their coefficients depend on the space-time dimension $D$. In addition, the integral coefficients depend on the dimensionality $D_{s}$ of the spin space for the loop particles FDH . For simplicity we neither display the dimensional dependence nor the natural dependence on kinematic variables. Furthermore, for the present discussion it is sufficient to consider fixed values of $D$ and $D_{s}$ for which the amplitude is finite. Also, we always consider planar amplitudes, as this is sufficient for our purposes in this paper (notice that for two-loop diagrams higher propagator powers appear only in planar amplitudes).

In a numerical approach, it is helpful to analyze eq. (II.1) prior to integrating over loop momenta, i.e., to analyze the integrand ${\cal A}(\ell_{l})$. The symbol $\ell_{l}$ represents the momenta of the two loops and will be used to denote quantities defined at the integrand level. The integrand is decomposed into terms that contribute to the sum in eq. (II.1), which we call master integrands, and independent surface integrands which integrate to zero NumUnitarity2 ,

[TABLE]

where $M_{\Gamma}$ and $S_{\Gamma}$ denote the set of master integrands and surface terms associated to diagram $\Gamma$, respectively. The numerator terms $m_{\Gamma,i}(\ell_{l})$ integrate either to master integrals or to zero,

[TABLE]

Similar notation has been used for example in Badger:2016ozq . For future reference, we define the integrand numerator ${N}(\Gamma,\ell_{l})$ associated with the propagator structure $\Gamma$ by

[TABLE]

The construction of the integrand representation in eq. (II.2) has been given in NumUnitarity2 , using appropriate integration-by-parts (IBP) identities IBPGKK . The IBP relations have to be chosen sufficiently general in order to include the propagator structures already present in the integrand ${\cal A}(\ell_{l})$ of the amplitude, e.g. given by Feynman rules. In particular, given that two-loop amplitudes contain diagrams with doubled propagators, we have to consider such propagator structures as well when constructing the sets of master integrands and surface terms.

II.3 Integrand coefficients and factorization

The coefficient functions $c_{\Gamma,i}$, can be obtained by solving the linear system of equations (II.2) for generic values of the loop momentum. In the generalized unitarity approach, the system of equations is analyzed diagram by diagram. For each diagram $\Gamma$, we consider the specific values of the loop momenta $\ell_{l}^{\Gamma}$ where internal particles go on-shell,

[TABLE]

In the limit $\ell_{l}\rightarrow\ell_{l}^{\Gamma}$, both sides of eq. (II.2) diverge and the coefficients of the poles can be compared, yielding a refined system of equations. Importantly, unitarity and factorization properties of field theory amplitudes imply that the leading coefficients of the poles in ${\cal A}(\ell_{l})$ are given by products of tree amplitudes.

The required tree amplitudes can also be indexed by the diagrams $\Gamma$. Let $T_{\Gamma}$ denote the set of tree amplitudes constructed in the following way: to each $n$-point vertex of $\Gamma$, one associates an $n$-point tree amplitude evaluated on the on-shell momenta $\ell_{l}^{\Gamma}$, with matched quantum numbers on internal lines. Notice that not all diagrams $\Gamma$ give rise to a well defined product of tree amplitudes. An example of this is displayed in fig. 2. The tree amplitude associated to the four-point vertex on the top right corner of diagram (b) behaves as $1/\rho$ and thus is divergent and ill-defined in the on-shell phase space of the diagram (which sets $\rho=0$). We will denote the subset of all diagrams which give rise to a well defined product of tree amplitudes by

[TABLE]

Note that $\Delta^{\prime}$ inherits the notion of hierarchy from $\Delta$. The hierarchies in $\Delta^{\prime}$ are referred to as cut hierarchies. It is a general feature that $\Delta\neq\Delta^{\prime}$ whenever one of the diagrams $\Gamma\in\Delta$ has a propagator structure with multiple copies of a given propagator, like diagram (a) in fig. 2. In fig. 3, to the left of the dashed line, we show the sunrise hierarchy for a massless 2 $\to$ 2 amplitude. The diagrams not belonging to the corresponding cut hierarchy are drawn inside a box.

Generalized unitarity builds on the observation that, in the limit $\ell_{l}\to\ell_{l}^{\Gamma}$, we have

[TABLE]

and that in this limit $R(\Gamma,\ell_{l}^{\Gamma})$ is given as a product of trees,

[TABLE]

We stress the fact that $R(\Gamma,\ell_{l}^{\Gamma})$ is only defined on the on-shell phase space of $\Gamma$. Naively, in generalized unitarity, one expects one equation (II.7) for each diagram, such that each numerator $N(\Gamma,\ell)$, as defined in eq. (II.4), is associated to its individual on-shell limit. However, starting at two-loops, $\Delta^{\prime}\neq\Delta$ and eq. (II.8) is well defined only for the subset of diagrams $\Gamma\in\Delta^{\prime}\subset\Delta$. We are thus left with less equations to determine the coefficients in eq. (II.2). In the following subsection we outline our algorithm to overcome this issue, and in section III we apply it in the context of concrete two-loop examples.

II.4 Leading and subleading poles in generalized unitarity

The diagrams in $\Delta^{\prime}$ correspond to leading poles of the amplitude in the on-shell limit, while those in $\Delta\setminus\Delta^{\prime}$ correspond to subleading poles. We start by reviewing the standard approach to deal with leading poles.

II.4.1 Extracting leading poles

Consider for simplicity the integrand of a maximal diagram $\Gamma$, that is a configuration in which $\Gamma$ contains the maximal number of edges required for the amplitude $\cal{A}$. When working in $D$ dimensions and with external momenta defined in four dimensions, the maximum number of edges is bounded for a two-loop amplitude with $n$ external particles by ${\rm min}(n+3,11)$. Furthermore, each subloop can contain at most 6 edges (7 if a doubled propagator is present). On the on-shell phase space of $\Gamma$, we get from eq. (II.7) that

[TABLE]

Through unitarity, we can directly compute $N(\Gamma,\ell_{l}^{\Gamma})$ as a product of trees, see eq. (II.8). Using eq. (II.4), we can then extract the corresponding set $\{c_{\Gamma,i}\}$ of integrand coefficients by sampling eq. (II.9) over enough points in the on-shell phase space $\ell_{l}^{\Gamma}$. We thus obtain $N(\Gamma,\ell_{l})$ for generic $\ell_{l}$.

Consider now a next-to-maximal integrand, that is an integrand with one less propagator than a maximal one. We can still use a relation similar to eq. (II.9), taking care of subtracting contributions coming from integrands with more propagators. For concreteness, let $\Gamma_{\text{NM}}$ be a next-to-maximal diagram. We denote by $\Gamma_{\text{NM},k}$ the ancestor of $\Gamma_{\text{NM}}$ which has the same propagators as $\Gamma_{\text{NM}}$ plus an extra one, $1/\rho_{k}$. At this stage we assume that $k\notin P_{\Gamma_{\rm NM}}$ for any $k$. On the on-shell phase space of $\Gamma_{\text{NM}}$, unitarity ensures

[TABLE]

in which the inverse propagators $\rho_{k}$ are evaluated on the momenta $\ell_{l}^{\Gamma_{\text{NM}}}$. The coefficients $\{c_{\Gamma_{\text{NM}},i}\}$ in the numerator $N(\Gamma_{\text{NM}},\ell_{l})$ are determined from a linear system of equations obtained from eqs. (II.4) and (II.10).

In the absence of subleading poles, we can iterate this procedure. A systematic extraction of all integrand coefficients of a given amplitude, see eq. (II.2), can then be carried out “level by level”, from the maximal integrands to the minimal ones.

II.4.2 Extracting subleading poles

Let us now consider the case where subleading singularities are present and develop an algorithm to extract their contributions. Let $\Gamma_{p}$ and $\Gamma_{d}$ be a pair of parent-daughter diagrams, $\Gamma_{p}>\Gamma_{d}$, such that the inverse propagator $\rho_{s}$ appears both in $\Gamma_{p}$ and $\Gamma_{d}$, but to a higher power in $\Gamma_{p}$ than in $\Gamma_{d}$. At two loops, it is sufficient to assume that it is squared in $\Gamma_{p}$. According to the definition of section II.3, $\Gamma_{p}\in\Delta^{\prime}$ but $\Gamma_{d}\notin\Delta^{\prime}$.

Then, the on-shell phase space defined by the two diagrams is the same, schematically $\ell_{l}^{\Gamma_{p}}=\ell_{l}^{\Gamma_{d}}$, and we thus have both a leading (from $\Gamma_{p}$) and sub-leading pole contribution (from $\Gamma_{d}$) in the on-shell limit of the integrand,

[TABLE]

where $R(\Gamma_{d},\ell_{l}^{\Gamma_{d}})$, for $\Gamma_{d}\notin\Delta^{\prime}$, has been implicitly defined as the subleading term in the $\ell_{l}\to\ell_{l}^{\Gamma_{d}}$ limit. We stress again, that this definition applies only on the on-shell phase space of $\Gamma_{d}$.

The term $R(\Gamma_{p},\ell_{l}^{\Gamma_{p}})$ is the leading term in the $\ell_{l}\to\ell_{l}^{\Gamma_{d}}$ limit, and is obtained from a product of trees as in eq. (II.8). From it we can determine the associated numerator ${N}(\Gamma_{p},\ell_{l})$ in the standard way discussed in the previous section. In contrast, an equivalent expression for $R(\Gamma_{d},\ell_{l}^{\Gamma_{d}})$ is not known and the determination of the corresponding numerator ${N}(\Gamma_{d},\ell_{l})$ must proceed differently.

We carry on as follows. We go down the $\Delta$ hierarchy until we find a diagram that has all the same propagators as $\Gamma_{d}$ except $1/\rho_{s}$. In our case it is sufficient to consider a corresponding daughter diagram of $\Gamma_{d}$ (in turn, a granddaughter of $\Gamma_{p}$). Let $\Gamma^{\prime}$ be such a diagram, $s\notin P_{\Gamma^{\prime}}$. Assuming that the only subleading poles are associated with the propagator $s$ (that is, $\Gamma^{\prime}\in\Delta^{\prime}$ by construction), the factorization limit of the amplitude as $\ell_{l}\to\ell_{l}^{\Gamma^{\prime}}$ leads to

[TABLE]

with $R(\Gamma^{\prime},\ell_{l}^{\Gamma^{\prime}})$ given by a product of trees, see eq. (II.8). We call eq. (II.12) the cut equation, and there exists one cut equation for each element of $\Delta^{\prime}$. The sum over diagrams $\Gamma$ runs over all ancestors of $\Gamma^{\prime}$. In a standard unitarity approach and in the absence of subleading poles, all numerators $N(\Gamma,\ell_{l}^{\Gamma^{\prime}})$ for $\Gamma>\Gamma^{\prime}$ will have been determined previously from their own cut equations. We would thus use eq. (II.12) to determine the numerator $N(\Gamma^{\prime},\ell_{l}^{\Gamma^{\prime}})$, see e.g. eq. (II.10) where this is done explicitly. In the presence of subleading poles, some numerators $N(\Gamma,\ell_{l}^{\Gamma^{\prime}})$ cannot be determined from an associated cut equation, as already discussed below eq. (II.8). We thus separate the terms of the sum in eq. (II.12) into two sets: those in $\widetilde{\Delta}$ whose numerators are still unknown, and those in $\Delta\setminus\widetilde{\Delta}$ which have already been determined. We then rewrite eq. (II.12) as:

[TABLE]

In this expression, all numerator terms on the right-hand side can be extracted from the standard generalized unitarity approach outlined in section II.4.1. All numerator terms on the left-hand side are still to be determined, and we will do so for all at once. More precisely, we sample eq. (II.13) over enough points of the on-shell phase space $\ell_{l}^{\Gamma^{\prime}}$ to build a system of equations big enough to determine all coefficient functions $\{c_{\Gamma^{\prime},i}\}$ and $\{c_{\Gamma,i}\}$ for all $\Gamma>\Gamma^{\prime}$, $\Gamma\in\widetilde{\Delta}$.

We note that the algorithm that we have proposed relies only on the unitarity of the theory under consideration. In particular, it extends trivially to any loop order and is entirely process independent. Indeed, while for generic multi-loop amplitudes the structure of subleading poles is in general much richer, with for example more than one subleading term in on-shell limits like in eq. (II.11), mixed subleading poles associated to different propagators, and non-planar configurations, our algorithm still allows to find enough suitable cut equations to solve for all unknown numerators as in eq. (II.13).

III Applications to four-point gluon amplitudes

In this section we apply the algorithm introduced in the previous section to extract coefficients of ancestor diagrams on the phase space of its descendants in the context of one- and two-loop examples. The one-loop example is included as a simple illustration of our algorithm, to show that it can also be used in the absence of subleading poles. In the two-loop example, we apply it to a case with subleading poles where a standard generalized unitarity approach would not be enough. Finally, we discuss the implementation of our approach in a numerical framework and the checks that we have performed on its applicability.

III.1 Box coefficients from the triangle phase-space

We compute box and triangle coefficients from triple cuts only. The system of equations which arises is less diagonal than a standard one-loop approach, but nonetheless tractable. Consider the expression for a triangle cut of a four-point amplitude at one-loop, a specific example of the next-to-maximal case described in eq. (II.10),

[TABLE]

where $\rho$ is the inverse propagator that was pinched to obtain the triangle diagram from the box diagram. The propagator is evaluated on the on-shell momentum $\ell^{\rm tri}$. We have used a pictorial representation to show the contributing diagrams. This is the standard cut equation for the triangle and the (parent-)box cuts, however, we assume that the box coefficient functions have not yet been determined, i.e., the box diagrams belong to $\widetilde{\Delta}$ in eq. (II.13). To proceed we insert the corresponding expressions for the box and triangle numerator functions, see eq. (II.4), leading to

[TABLE]

in which the right-hand side is given by a product of trees. The triangle and box coefficients are written as $c_{{\rm tri},i}$ and $c_{{\rm box},i}$, respectively. The associated numerator insertions are denoted by $m_{{\rm tri},i}(\ell)$ and $m_{{\rm box},i}(\ell)$. The number of master and surface integrands of the triangle and box diagrams have been denoted by $m$ and $n$, respectively. Compared to the notation in eq. (II.2), $m$ is the number of elements of $M_{\rm tri}\cup S_{\rm tri}$ and $n$ the number of elements of $M_{\rm box}\cup S_{\rm box}$. All the triangle and box coefficients are then found by sampling the triangle cut over $n+m$ momenta on the on-shell phase space $\ell^{\rm tri}$ and then solving for $c_{{\rm box},i}$ and $c_{{\rm tri},i}$ by a linear regression. This one-step approach requires solving a single large linear system of equations, compared to two smaller ones when solving first for box coefficients and subsequently for the triangle coefficients.

III.2 The bubble-box hierarchy at two loops

As an example of the application of our algorithm in the presence of subleading poles, we consider a 2 $\to$ 2 amplitude in massless QCD. The maximal level diagrams have seven propagators, and the minimal diagrams are the sunrise diagrams with three propagators (see fig. 3).

Consider now the bubble-box hierarchy shown in fig. 4. All diagrams that appear above it are associated with factorization limits of (h). Six out of the seven ancestor diagrams shown have associated products of trees and their integrand coefficients can be directly extracted. Diagram (e), on the other hand, represents subleading pole contributions to the doubled-propagator diagram (a).

As an aside, we note that for all numerators in fig. 4, apart from the double-box (c) and the bubble-box (h), the integrand function space is spanned entirely by surface terms. Two master integrals are associated with (c) and one with (h). With this in hand, we solve the cut hierarchy to obtain the coefficients.

Consider first the integrands associated with the maximal diagrams (a)-(d). For example, for the numerator of (c), we have:

[TABLE]

where the on-shell momenta of diagram (c) are denoted by $\ell_{l}^{c}$. We have used a pictorial representation to denote the corresponding diagram. Analogous equations hold for the diagrams (a), (b) and (d). In practice, for each maximal diagram $\Gamma$ one generates a linear set of equations by inserting sufficiently many on-shell momentum values for $\ell_{l}^{\Gamma}$, and solves for the integrand coefficients in $N(\Gamma,\ell_{l})$.

We move then to the numerators of the two next-to-maximal diagrams, (f) and (g). These are found by their corresponding cut equations, as in eq. (II.10). For example, the numerator for the box-triangle diagram (f), fulfills

[TABLE]

where $\rho_{\rm fb}$ and $\rho_{\rm fc}$ denote the propagators that are pinched to obtain diagram (f) from (b) and diagram (f) from (c), respectively (see fig. 4). We denote the on-shell momenta associated to diagram (f) by $\ell_{l}^{\rm f}$. The integrand corresponding to diagram (g) is treated in the same manner.

Finally, we proceed to solve for the coefficients associated to diagrams (e) and (h), which involves the extraction of subleading poles. Setting $\widetilde{\Delta}=\{(e)\}$ in eq. (II.13) we obtain

[TABLE]

where the inverse propagators $\rho_{ij}$ denote the propagator pinched to obtain diagram ($i$) from diagram ($j$), with the indices ($i$) and ($j$) corresponding to pairs of daughter-parent diagrams in fig. 4. Each numerator ${N}(\Gamma,\ell_{l}^{\rm h})$ is written in terms of its basis of integrands as in eq. (II.4). At this stage, the terms on the right-hand side of eq. (III.4) are known, and we can directly solve for the coefficients associated to the integrand {N}\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/BubbleBoxGeneric.pdf}}},\ell_{l}^{\rm h}\right) as well as to the subleading-pole integrand {N}\left(\vbox{\hbox{\includegraphics[scale={0.1}]{diagrams/BubblePentagonSemiGeneric.pdf}}},\ell_{l}^{\rm h}\right) by sampling over enough on-shell loop momentum configurations.

III.3 Numerical implementation and checks

We discuss in this section the checks that we have performed on the applicability of our algorithm in a numerical unitarity framework. We have numerically reproduced the analytic results of fourgluons for all the master coefficient functions in all planar two-loop four-gluon helicity amplitudes. This was achieved through the following steps:

• •

We have implemented $D$-dimensional off-shell recursions BGrec to compute all required trees and (one- and) two-loop cuts for generic computations in $D$-dimensional numerical unitarity. Our implementation allows general values of the $D_{s}$ parameter, the dimensionality of the spin space for the loop particles. The numerical trees and one-loop cuts were cross checked against the BlackHat library BlackHat in four dimensions. We have also cross checked our two-loop cuts against results obtained using an independent implementation of the gluon amplitudes in six dimensions. Systematic checks of Ward identities and factorization limits of the trees and cuts have been performed.

• •

Two independent implementations for the construction of the $\Delta$ hierarchies and the associated subtraction structures have been produced, one based on planar configurations and another following the color decomposition of Ochirov:2016ewn . Both have been cross checked, and independently shown to correctly produce subtraction terms for diagrams in $\Delta^{\prime}$.

• •

In order to compare against known results for two-loop four-point gluon amplitudes, we have produced a set of master-surface integrand decompositions as in eq. (II.2) for all the diagrams $\Gamma$ in those amplitudes. The decompositions were produced along the lines of ref. NumUnitarity2 . We cross-checked that the surface integrands we constructed integrate to zero with the generator of IBP relations FIRE Smirnov:2014hma . We have also shown that the integrand decompositions fully span their corresponding integrand spaces by comparing them to an alternative parametrization in terms of tensor insertions tensorBasis .

• •

With all the tools described above, we were able to numerically compute all integrand coefficients that contribute to the planar two-loop four-point gluon amplitude, see fig. 3. We have validated the values of the master integral coefficients with the known analytic results fourgluons for all helicity configurations.

Most checks have been performed on a small set of phase-space points and we delay any systematic efficiency and stability checks of our implementation to future work. Nevertheless it is worth mentioning that we find that for example the master coefficient functions of the bubble-box hierarchy in fig. 4 can be extracted in less that 100$ms$ on generic phase-space points. Typically we find that they agree with numerical values obtained from the analytic expressions fourgluons to an accuracy of better than 10 digits. This study has been performed for fixed values of the dimensional parameters $D$ and $D_{s}$ using only double-precision arithmetics.

IV Conclusions

In this article we have presented an algorithm for extracting subleading-pole contributions in two-loop scattering amplitudes through numerical unitarity techniques, which are required to obtain the full amplitude. Subleading-pole integrand coefficients can be obtained by solving linear systems of equations built up from cut equations of associated descendant diagrams. The algorithm is process independent and can be naturally generalized to multi-loop amplitudes. We have performed a number of consistency checks, in which multiple subleading-pole contributions have been extracted at several levels of a two-loop cut hierarchy. In addition, we computed two-loop master coefficients through numerical unitarity, and have confirmed our results by comparing to available analytic results for four-point two-loop gluon amplitudes. Although in numerical unitarity one performs calculations in fixed dimensions $D$ and $D_{s}$, the regressions to general values of those parameters can be achieved by the observation that functionally they appear as rational or polynomial functions. We have reconstructed the full $D$ and $D_{s}$ dependence of all integral coefficients in the sunrise hierarchy in fig. 3. All these results were used in the first calculation of a two-loop amplitude in the framework of numerical unitarity 4g2loop .

A systematic study of the efficiency and stability of the numerical unitarity approach is left to future work, however, we have observed that the algorithm appears sufficiently fast and numerically stable. We hope in the future to explore the use of the numerical unitarity method to two-loop amplitudes with more than four external particles.

The algorithm we presented to deal with subleading-pole contributions requires to handle an enlarged set of linear equations when computing integral coefficients. It would be interesting to explore alternative approaches which organize the cut equations more effectively. For example, it is possible to directly associate subleading contributions to subtracted gluon amplitudes to tree amplitudes with graviton exchange Stieberger:2015kia . Further ideas include accessing the subleading poles through numerical limits, or using so-called BCJ relations Bern:2008qj to relate the subleading contributions to other color-ordered cuts.

Acknowledgments

We thank Z. Bern, A. de Freitas and D.A. Kosower for helpful discussions. We particularly thank Z. Bern for providing analytical expressions from ref. fourgluons . S.A.’s work is supported by the Juniorprofessor Program of Ministry of Science, Research and the Arts of the state of Baden-Württemberg, Germany. H.I.’s work is supported by a Marie Skłodowska-Curie Action Career-Integration Grant PCIG12-GA-2012-334228 of the European Union. The work of F.F.C., M.J. and B.P. 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. This work was performed on the bwUniCluster funded by the Ministry of Science, Research and the Arts Baden-Württemberg and the Universities of the State of Baden-Württemberg, Germany, within the framework program bwHP. The authors are grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and its partial support during the completion of this work.