Two-Loop Four-Gluon Amplitudes with the Numerical Unitarity Method

TL;DR

This paper introduces a numerical unitarity method to compute two-loop four-gluon amplitudes, enabling automated, multi-scale calculations relevant for collider physics.

Contribution

It presents the first numerical approach for two-loop amplitudes using unitarity, demonstrating its universality and potential for automation in collider phenomenology.

Findings

Successful numerical computation of two-loop four-gluon amplitudes

Analysis of the method's numerical stability and properties

Reconstruction of analytic amplitude forms

Abstract

We present the first numerical computation of two-loop amplitudes based on the unitarity method. As a proof of principle, we compute the four-gluon process. We discuss the new method, analyze its numerical properties and apply it to reconstruct the analytic form of the amplitudes. The numerical method is universal, and can be automated to provide multi-scale two-loop computations for phenomenologically relevant signatures at hadron colliders.

Equations

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

$$\mathcal{A}$$

$$\sum _{\mathrm{states}}\prod _{i\in T_{\Gamma }}{\mathcal{A}}_i^{\mathrm{tree}}({\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 })},$$

$$0=\int \prod _{l=1,2}{d}^D{\ell }_l\frac{\partial }{\partial {\ell }_j^{\nu }}\left[\frac{u_j^{\nu }}{{\prod }_{k\in P_{\Gamma }}{\rho }_k}\right],$$

$$u_i^{\nu }\frac{\partial }{\partial {\ell }_i^{\nu }}{\rho }_j=f_j{\rho }_j,$$

$${\ell }_l=\sum _{j\in B_l^p}{v}_l^j{r}^{lj}+\sum _{j\in B_l^t}{v}_l^j{\alpha }^{lj}+\sum _{i\in B^{ct}}{n}^i{\alpha }^{li}+\sum _{i\in B^{ϵ}}{n}^i{\mu }_l^i,$$

$$r^{li}=-\frac{1}{2}({\rho }_{li}-(q_{li}{)}^2-{\rho }_{l(i-1)}+(q_{l(i-1)}{)}^2),$$

$${\mu }_{ll}\equiv ({\mu }_l{)}^2$$

$$u=f_i{\rho }_i\frac{\partial }{\partial {\rho }_i}+u_j\frac{\partial }{\partial {\alpha }^j}+\sum _{l,l^{\prime }=1,2}{f}_{l^{\prime }}^l{\mu }_l^k\frac{\partial }{\partial {\mu }_{l^{\prime }}^k},$$

$$\begin{array}{rl}\bar{u}({\mu }_{11})& =2{\mu }_{11}{f}_1^1+2{\mu }_{12}{f}_1^2,\\ \bar{u}({\mu }_{22})& =2{\mu }_{22}{f}_2^2+2{\mu }_{12}{f}_2^1,\\ \bar{u}({\mu }_{12})& ={\mu }_{12}(f_1^1+f_2^2)+{\mu }_{11}{f}_2^1+{\mu }_{22}{f}_1^2.\end{array}$$

$$\begin{array}{rl}& {\alpha }^{lj}{\alpha }^{l^{\prime }j}-b_1\frac{{\mu }_{l{l}^{\prime }}}{ϵ},l,l^{\prime }\in \{1,2\},\\ & ({\alpha }^{lj}{)}^2({\alpha }^{l^{\prime }j}{)}^2-b_1\frac{{\mu }_{ll}}{ϵ}({\alpha }^{l^{\prime }j}{)}^2-b_2\frac{{\mu }_{l{l}^{\prime }}}{ϵ}{\alpha }^{lj}{\alpha }^{l^{\prime }j},l\ne l^{\prime },\end{array}$$

$$m_{\Gamma ,u}({\ell }_l)$$

Full text

The experiments at the Large Hadron Collider (LHC) at CERN are entering a new phase in which observables will be studied with relative errors of the order of a few percent. New discoveries through precision measurements require an equal or better control over the theoretical uncertainties of predictions from the Standard Model of particle physics. A central bottleneck to obtaining such predictions is the complexity of computing quantum corrections. We demonstrate a new, automatable algorithm for computing two-loop corrections in QCD based on the established unitarity method. In particular we focus on the numerical variant of the method which has proven valuable for dealing with multi-scale problems. In addition to the flexibility of this approach, the geometric nature of our method simplifies intermediate computational steps and promises good numerical behavior.

As a proof of principle of the new method we recompute the two-loop gluon-gluon scattering amplitudes 4gluonGOT ; 4gluonBFD . This process exposes much of the complexity of two-loop scattering amplitudes in QCD, including their universal infrared and ultraviolet behavior. The expressions in ref. 4gluonBFD have been obtained using the analytic variant of the unitarity method Unitarity which constructs scattering amplitudes from their unitarity and analytic properties. This method has been applied to a number of one-loop computations for the LHC, and is the method of choice in formal research on scattering amplitudes in supersymmetric theories. In QCD, at two-loop level, the lack of (super)symmetry and the need for an infrared and ultraviolet regulator pose further challenges due to the appearance of one-loop sub divergences. Nevertheless, analytic computations of five- and six-gluon amplitudes  AllPlusTwoLoop have recently become available, albeit for constrained helicity configurations. Numerical variants of the unitarity approach at one-loop level OPP ; NumUnitarity ; DNumUnitarity ; BlackHat have by now provided a large number of phenomenologically relevant predictions for the LHC. The flexibility of the numerical unitarity method in combination with its good numerical stability motivate us to extend it to multi-loop amplitudes. The algorithm which we put forward generalizes the one-loop approach in a non-trivial way, as the two-loop variant requires additional geometric input IntDec . Being automatable and less susceptible to the complexity of analytic multi-scale computations, our approach has the potential to mirror the successes found at one-loop level. In addition, analytic expressions can be efficiently reconstructed from a numerical algorithm Peraro:2016wsq .

In this letter, we present the first numerical computation of two-loop QCD amplitudes with the unitarity method. We focus on the leading-color contributions to the two-loop four-gluon amplitudes which we validate by comparing with known results 4gluonBFD . First, we set up the equations necessary for a hierarchical extraction of an amplitude’s integrand DProp . Second, we decompose integrands of massless four-point amplitudes into master integrals and surface terms, extending the results of IntDec . Third, we describe our numerical implementation and the linear algebra techniques employed to compute coefficients for fixed values of the dimensional regulator. Subsequently, we reconstruct the full regulator dependence and show explicit numerical results. Finally, we discuss numerical reconstruction of the analytic amplitude and give our conclusions.

Numerical Unitarity Method. We apply a variant of the unitarity method suitable for analytic and numerical computations which generalizes one-loop methods to higher loop orders. For more details of our approach we refer the reader to ref. DProp . We start with an ansatz for the integrand ${\cal A}(\ell_{l})$ of a two-loop amplitude IntDec ,

[TABLE]

where $\Delta$ denotes the set of two-loop diagrams which specify the possible propagator structures of the amplitude. In fig. 1 we show the set $\Delta$ corresponding to planar two-loop four-point massless amplitudes. The set of diagrams $\Delta$ is organized hierarchically with the partial ordering $\Gamma_{1}>\Gamma_{2}$ if $P_{\,\Gamma_{1}}\supset P_{\,\Gamma_{2}}$, i.e., if the diagram $\Gamma_{2}$ (a descendant) is obtained from $\Gamma_{1}$ (an ancestor) by removing one or more propagators. Here $P_{\Gamma}$ denotes the set of propagators associated to diagram $\Gamma$. The loop momenta are denoted by $\ell_{l},\,l=1,2$ and inverse propagators by $\rho_{j}$. The numerators span the full set of independent numerator terms NumPar and are restricted by power counting. For each diagram, the set $M_{\Gamma}$ specifies numerators associated to master integrals and $S_{\Gamma}$ the set of surface terms (i.e., terms that integrate to zero) necessary to fully parameterize the numerators. The $m_{\Gamma,i}(\ell_{l})$ are polynomials in $\ell_{l}$. We work in dimensional regularization such that these, as well as the coefficient functions $c_{\Gamma,i}$, depend on the loop-momentum dimension $D=4-2\epsilon$. Dependence on external kinematics is implicit. For simplicity we omit color information and focus on leading-color contributions.

The physical scattering amplitude is obtained in terms of master integrals $I_{\Gamma,i}$,

[TABLE]

as surface terms drop out from eq. (1) after integration. At this stage, we only need to compute the coefficient functions $c_{\Gamma,i}$. In generalized unitarity, they are determined by solving the linear system of equations (1) for values of the loop momenta $\ell_{l}^{\Gamma}$ for which the internal particles go on-shell, $\rho_{j}\rightarrow 0$ for all $j\in P_{\Gamma}$. In the absence of subleading poles associated with this limit DProp , putting the ansatz (1) and factorization limits together we obtain

[TABLE]

which we call the cut equations. The set $T_{\Gamma}$ labels all tree amplitudes corresponding to the vertices in the diagram $\Gamma$. The state sum runs over all possible internal states, here the $(D_{s}-2)$ gluon helicity states for each internal line of $\Gamma$.

The cut equations apply only for the subset of diagrams $\Gamma\in\Delta^{\prime}$ for which the factorization limit has no subleading poles, i.e., not all propagator structures yield a cut equation. Nevertheless, it is possible to organize the set of cut equations to produce enough linear systems of equations to compute all coefficient functions $c_{\Gamma,i}$ DProp .

This variant of the unitarity approach for computing multi-loop amplitudes is suitable both for analytic and numerical calculations. In what follows, we will focus on its numerical implementation.

Construction of integrand parameterization. For the construction of the integrand (1) we extend the method of ref. IntDec to massless propagators and external momenta. The ansatz (1) relies on a suitable set of integration-by-parts (IBP) relations,

[TABLE]

which we use to construct the numerators $m_{\Gamma,i}(\ell_{l})$. In order to control propagator powers to match the ones of the integrand (1), we construct IBP relations from constrained vector fields IBPGKK which solve,

[TABLE]

where no summation over the index $j$ is implied. Solutions $\{u_{i}^{\nu}\}$ of eq. (5) are referred to as IBP-generating vectors. Alternatively, reduction programs IntDecProg can be used to obtain appropriate IBP relations.

An important tool for solving eq. (5) and analyzing IBP relations are the natural coordinates Baikov ; IntDec ; LZ ,

[TABLE]

which parameterize the strand momenta $\ell_{l}$ ($l=1,2,3$, see fig. 2) in terms of inverse propagators $\rho_{li}$, and auxiliary variables $\alpha^{li}$ and $\mu_{l}^{i}$. The vectors $q_{li}$ are fixed by momentum conservation after imposing $q_{10}=q_{20}=0$. The vectors $n^{i}$ form an orthonormal basis transverse to the scattering plane, i.e $n^{i}\cdot p_{j}=0$. Labels in $B^{\epsilon}$ refer to directions beyond four dimensions and labels in $B^{ct}$ denote transverse directions within four dimensions. Within the scattering plane, we setup distinct basis vectors $v_{l}^{i}$ for each strand of the diagram; the vectors $p_{i}$ with ${i\in B_{l}^{p}}$ which exit the strand $l$ are completed with additional physical momenta $p_{i}$ with ${i\in B_{l}^{t}}$, so as to span the whole scattering plane. The vectors $v_{l}^{i}$ are dual to the $p_{j}$, $v_{l}^{i}=(G_{l})^{ij}p_{j}$, with $(G_{l})^{ij}$ the inverse of the Gram matrix $(G_{l})_{ij}=p_{i}\cdot p_{j}$. The new coordinates $\alpha^{li}$, $\mu_{l}^{i}$ and $\rho_{li}$ are not independent: momentum conservation $\ell_{1}+\ell_{2}+\ell_{3}+p_{b}=0$ removes redundant $\alpha_{li}$ and $\mu_{3}^{i}$ variables. We often label $\rho$ variables and independent $\alpha$ variables by a single subscript referring to the pair of strand index $l$ and label $i$. Furthermore, we have the constraint

[TABLE]

for the extra-dimensional momentum squares.

In these coordinates, IBP-generating vectors are IntDec

[TABLE]

where we denote basis vectors along coordinate lines by the respective partial derivatives. Repeated indices are summed over. We impose rotation symmetry in the extra-dimensional components as manifest in the integrand (1). The vectors must be polynomial when written in canonical coordinates of $\ell_{l}$, requiring first that the coordinate functions $f_{i}$, $f_{l}^{l^{\prime}}$ and $u_{i}$ to be polynomial in $\alpha^{i}$ and $\rho_{i}$, and second that the vectors are compatible with the relations in eq. (8). After redundant variables have been eliminated, the compatibility conditions are quadratic in $\alpha^{i}$ and $\rho_{i}$, and read

[TABLE]

We use $\bar{u}$ to denote the $\mu$-independent directions of $u$ in eq. (9) and $\mu_{12}\equiv(\mu_{33}-\mu_{11}-\mu_{22})/2$. The generating set of solutions for the unknowns $f_{i},f_{l^{\prime}}^{l}$ and $u_{j}$ are obtained using Singular Singular which we use in addition to universal vectors IntDec . We consider the solutions that do not vanish on-shell and keep a generating set of these. An automated implementation of the on-shell variant of the above equations is also available Azurite .

The simplest IBP-relations are obtained from diagonal rotations IntDec . The respective relations are generalizations of the one-loop numerators OPP ; NumUnitarity being traceless completions of the transverse monomials. They can be systematically constructed from traceless completions

[TABLE]

with $j\in B^{ct}$. After multiplying with additional $\alpha$’s either from $B^{t}_{l}$, $B^{t}_{l^{\prime}}$ or $B^{ct}$ the rational constants $b_{1}$ and $b_{2}$ are fixed to yield expressions that are transverse to the Lorentz structures after loop integration. Monomials with odd powers of $\alpha$’s from $B^{ct}$ integrate to zero and can be used for numerators as they are. The completion of the IBP relations is obtained by multiplying the IBP-generating vectors with the independent $\alpha$ variables. The numerators are,

[TABLE]

taking into account the Jacobian factor from eq. (4), see e.g. Zhang:2016kfo . In eq. (12), $n_{\alpha}$ denotes the number of independent $\alpha_{i}$-variables and $\nu_{i}$ the propagator power.

The numerator decomposition into surface terms and master integrands must be non redundant and complete. The spanning set of IBP-relations is obtained modulo inverse propagators, which allows to check the linear independence with on-shell conditions $\rho_{i}=0$ imposed. Finally, the numerators of master integrands are given by the span of irreducible numerators modulo surface terms. In natural coordinates, the irreducible numerators are the monomials in $\alpha^{i}$ consistent with power counting. It is sufficient to compare the span of IBP-relations and irreducible numerators on a given phase-space point with on-shell conditions imposed to determine a decomposition in terms of master integrals as in eq. (2) IntDec ; LZ .

Numerical implementation. We first decompose the amplitude into master integrals and surface terms as in eq. (1) with the method described above. We find two master integrals associated with the first diagram in the first line of fig. 1 (the double-box) and one master integral associated with each of the first two diagrams of the third line, the first two diagrams of the fourth line, and the diagram in the last line. For the double-box diagram we choose the scalar and the (irreducible) numerator insertion, and for all other diagrams we choose the scalar integrals.

The coefficients $c_{\Gamma,i}$ in eq. (3) must now be fixed. They are determined numerically by sampling of the on-shell phase spaces $\ell_{l}^{\Gamma}$, for all $\Gamma\in\Delta^{\prime}$ DProp . This was implemented in a C++ library, for which the required analytic information - color decomposition (see ref. Color2Loop ), hierarchy of cut equations organized to handle subleading poles and IBP relations - is produced in Mathematica. For the four-dimensional spinor algebra we have used tools from the BlackHat library BlackHat .

We built an (over)constraining system of linear equations for $c_{\Gamma,i}$ by computing the subtracted product of trees (through off-shell recursions BGrec ) at randomly sampled phase-space points. Though analytically $D$-dimensional, phase-space samples can be constructed with 6-dimensional momenta due to the rotational invariance beyond the 4-dimensional physical slice. The on-shell phase spaces are generated by nested one-loop parameterizations BlackHat ; NumUnitarity . We single out one of the loops and construct loop momenta on its on-shell phase space. We then input these momenta into the second loop to find two-loop on-shell configurations.

In order to numerically solve the linear systems we constructed for the $c_{\Gamma,i}$, we employ standard linear algebra techniques as implemented in the LAPACK library LAPACK which uses the BLAS routines BLAS . For higher-precision arithmetics, we also use the associated routines from the MPACK libraries MPACK . For every $n\times n$ system of equations $Mc=a$ we employ a $PLU$ factorization of the square matrix $M$, in which $L$ is a lower-triangular matrix, $U$ is an upper-triangular matrix and $P$ is a permutation matrix. For over constrained systems of equations, we employ a $QR$ factorization ($Q$ a rectangular orthogonal matrix, and $R$ a squared upper-triangular one) for minimizing $||Mc-a||$. These factorizations allow for efficient and numerically stable solutions to the systems of equations, which are typically order 100-dimensional.

By solving these systems, we numerically determine the coefficients of master integrals for fixed values of $D$ and $D_{s}$. The $D_{s}$ dependence of master-integral coefficients is at most quadratic and can be reconstructed DNumUnitarity by evaluating the coefficients with state sums (3) for three different values of $D_{s}$. The $D$ dependence of the coefficients is rational and originates from the $D$-dependent IBP relations, see eqs. (11) and (12). A priori, the polynomial degrees of the numerator and denominator of the rational function are unknown. However, as this is phase-space independent, it can be determined once and for all in a dedicated run using rational reconstruction techniques in high precision Peraro:2016wsq ; Abramowitz . Further, because the coefficients of the denominators are rational numbers, they can be exactly reconstructed using continued fractions. In this way, after a warm-up phase, one needs only reconstruct a polynomial dependence of known rank.

Finally, to obtain values for the amplitudes (2) we combine the integral coefficients with our own implementation of the loop integrals Integrals4pt using the function library of ref. HPL . By setting $D_{s}=D=4-2\epsilon$ in the coefficients and expanding around $\epsilon=0$, we recover the value of the bare amplitude as a Laurent series in $\epsilon$ in the HV HV variant of dimensional regularization.

Validation of results. We have computed the two independent helicity configurations with non-vanishing tree amplitudes in the leading-$N_{c}$ color limit, omitting contributions from closed fermion loops. We normalize the results to the respective tree-level amplitudes ${\cal A}_{0}$. We provide numerical values for the coefficients of the Laurent series in $\epsilon$ in the HV scheme. For the point $p_{1,3}=-1/2\,(1,0,0,\pm 1)$, $p_{2,4}=1/2\,(1,0,\pm\sqrt{3}/2,\pm 1/2)$, we find the following values for the bare amplitudes ($g_{s}=1$ and $\mu=1$):

[TABLE]

We have checked that the above results agree with the expected universal infrared-pole structure CataniFormula , and that they match the results obtained from the analytic expressions of ref. 4gluonBFD . These comparisons validate the C++ implementation of our numerical unitarity algorithm.

In fig. 3, we show the stability of our calculation by looking at minus the base 10 logarithm of the relative error for the numerical calculation with respect to the analytic result. This is done over 10,000 phase-space points evenly distributed, with a minimum $p_{T}$ cut on the final-state partons set to $1/100$ the total energy in the center-of-mass frame. The numerical computation has employed single-double precision operations, except when extracting the maximal numerators for which we choose to use double-double precision arithmetics. We have also introduced a rescue system, based on comparisons to the known pole structure of the amplitude CataniFormula , and when failing we recompute the full amplitude in higher floating-point precision. We observe that our calculation is precise to better than a per mill relative error for the bulk of phase space and that the rescue system works properly by shifting the left tail of the distribution towards 16 digits to the right. When computing with single-double precision over phase space and after reconstruction of the dimensional regulators, our results for master integral coefficients are commonly accurate to 8 digits.

Finally, we have produced results with quadruple-double precision QD , which gives enough numerical precision to fully reconstruct the analytic form of the amplitudes from the purely numerical computations. This is achieved with the same techniques as for the regulator reconstruction since, after rescaling, the integral coefficients are a function of a single variable.

Conclusions. We have presented a method based on generalized unitarity for automated analytic and numerical computations of two-loop scattering amplitudes in QCD. We have validated our algorithm by comparing the numerical results to the ones of known analytic amplitudes. The method is numerically stable and refinements can improve this further in the future. Furthermore, our method can be used to extract the analytic expressions for the amplitude from the purely numerical results. In the last decade, similar developments at one-loop level have led to important predictions for high-multiplicity processes at hadron colliders. We are thus optimistic that our approach will be instrumental in yielding new two-loop matrix elements necessary for QCD phenomenology in the near future.

Acknowledgements: We thank Z. Bern, L.D. Dixon and D.A. Kosower for many helpful discussions. We particularly thank Z. Bern for providing analytic expressions from ref. 4gluonBFD . We thank C. Duhr for the use of his Mathematica package PolyLogTools. We thank J. Dormans and E. Pascual for technical discussions. H.I is grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and its partial support during the completion of this work. 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. The work of M.Z. is supported by the U.S. Department of Energy under Award Number DE-SC0009937. 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.