Two-loop leading colour QCD corrections to q q ¯ →γγg and qg→γγq
Bakul Agarwal, Federico Buccioni, Andreas von Manteuffel, Lorenzo Tancredi
MMSUHEP-21-005OUTP-21-04P
Two-loop leading colour QCD corrections to q ¯ q → γγg and qg → γγq Bakul Agarwal, a Federico Buccioni, b Andreas von Manteuffel, a Lorenzo Tancredi b a Department of Physics and Astronomy, Michigan State University,East Lansing, Michigan 48824, USA b Rudolf Peierls Centre for Theoretical Physics, University of Oxford,Clarendon Laboratory, Parks Road, Oxford OX1 3PU
E-mail: [email protected] , [email protected] , [email protected] , [email protected] Abstract:
We present the leading colour and light fermionic planar two-loop correc-tions for the production of two photons and a jet in the quark-antiquark and quark-gluonchannels. In particular, we compute the interference of the two-loop amplitudes with thecorresponding tree level ones, summed over colours and polarisations. Our calculation usesthe latest advancements in the algorithms for integration-by-parts reduction and multivari-ate partial fraction decomposition to produce compact and easy-to-use results. We haveimplemented our results in an efficient C++ numerical code. We also provide their analyticexpressions in Mathematica format.
Keywords:
QCD, Collider Physics, NNLO calculations a r X i v : . [ h e p - ph ] F e b ontents The production of a pair of isolated highly-energetic photons in proton-proton collisions, pp → γγ + X , represents an important class of processes for the physics programme at theLarge Hadron Collider (LHC) and at future hadron colliders. Among many relevant as-pects, pairs of prompt photons (diphotons) constitute an irreducible background to variousStandard Model (SM) and Beyond Standard Model (BSM) processes, most prominentlyHiggs boson production in its H → γγ decay channel. Indeed, with the increase of collisionenergy, the diphoton invariant mass distribution can provide a powerful tool to search forheavy resonances decaying to pairs of photons [1, 2], while its transverse momentum distri-bution offers a unique probe to investigate their properties. It is therefore crucial to providean accurate theoretical description of the production of a pair of photons recoiling againsthard Quantum Chromodynamics (QCD) radiation, across a vast spectrum of energies.At the technical level, the theoretical description of the production of a pair of photonswith large transverse momentum is non trivial for several reasons. An obvious one is thatit requires computation of 2 → n scattering amplitudes with n ≥ q ¯ q → gγγ , qg → qγγ , and¯ qg → ¯ qγγ , where the second and third channels can be obtained as crossings of the first.The loop-induced gluon-fusion process gg → gγγ , instead, contributes formally only atnext-to-next-to-leading-order (NNLO).The mathematical complexity of loop corrections to the scattering amplitudes aboveis the main reason why pp → γγ + jet is currently known only up to next-to-leading-order– 1 –NLO) QCD [3, 4]. Calculation of cross section and differential distributions for pp → γγ +jet through NNLO QCD requires the computation of two-loop amplitudes for q ¯ q → gγγ and its crossings, together with an efficient subtraction scheme to organise and cancel theinfrared (IR) divergences between real and virtual contributions. In recent years, a lot ofprogress has been achieved on both fronts. On the IR subtraction side, several schemeshave been developed which are able to handle NNLO QCD corrections to the production ofa colour singlet plus one QCD jet [5–11]. On the amplitude side, equally impressive resultshave been achieved. The first two-loop 2 → → all relevant master integrals(MIs) in terms of a well understood class of functions [43–47]. This has been achieved bythe use of the method of differential equations [48–53] augmented by the choice of a canon-ical basis [54, 55]. Alternative approaches to the reduction of two-loop amplitudes basedon ideas from the one-loop generalised-unitarity program as well as finite field arithmetichave also been very successful [12, 17, 18, 31, 56].In this paper, we consider the calculation of two-loop QCD corrections to the produc-tion of a diphoton pair and a jet for the partonic channels q ¯ q → gγγ and qg → qγγ , basedon a Feynman diagrammatic approach. We demonstrate how the use of new algorithms forthe reduction of loop integrals and multivariate partial fraction decomposition allows us tocompute these two-loop corrections in a rather straightforward way and to produce verycompact and efficient numerical implementations for them. While it is common practiceto compute helicity amplitudes for multi-loop processes, the case of diphoton productionplus jet does not require us to keep track of the polarisations of the external photons. Wetherefore expect it to be sufficient for near-term phenomenological applications to considerthe interference of the two-loop amplitudes with the corresponding tree-level amplitude,summed over colours and polarisations. We will show that, also in this case, very compactresults can be obtained, similar to what can be achieved for comparable helicity ampli-tudes. Together with the recently computed three-loop QCD corrections to diphoton pro-duction [57], these amplitudes also provide an essential ingredient towards the calculationof pp → γγ in N LO QCD.The rest of the paper is organised as follows. In Section 2 we describe the processesand their kinematics. In Section 3 we illustrate the general structure of the scatteringamplitudes, and in particular of the interference terms contributing to the squared ampli-tude. Technical aspects of the diagrammatic calculation, the integral reductions and themultivariate partial decompositions are presented in Section 4. In Section 5 we describeour renormalisation and infrared subtraction procedures which define the final form of theresults. In Section 6 we discuss how we optimise our analytic results by using a minimal set– 2 –f rational functions, and in Section 7 we present their implementation in a
C++ numericalcode. We finally draw our conclusions in Section 8.
We consider the production of a pair of photons in association with a gluon in quark-antiquark annihilation q ( p ) + ¯ q ( p ) → g ( p ) + γ ( p ) + γ ( p ) (2.1)up to two-loop order in massless QCD, i.e. corrections up to O ( α ) relative to the tree-level. We focus here on the scattering process in the physical region, i.e. a 2 → qg → qγγ , which can be obtained by crossing symmetry according to q ( p ) + g ( p ) → q ( p ) + γ ( p ) + γ ( p ) . (2.2)Analytic results for both channels, with the appropriate identification of the external mo-menta as in Eqs. (2.1) and (2.2), are provided in the ancillary files of the arXiv submissionof this paper. Kinematics are fixed by imposing that all external particles fulfil the on-shellcondition p i = 0, such that one is left with five independent kinematic invariants, whichwe choose to be the five adjacent ones, s = ( p + p ) , s = ( p − p ) , s = ( p + p ) ,s = ( p + p ) , s = ( p − p ) . (2.3)Note the relative sign between initial- and final-state momenta reflecting our kinematicconfiguration which implies for the physical region s , s , s > s , s <
0. All theother invariants follow by momentum conservation and can be derived from the independentones (2.3) via s = s − s − s , s = s − s − s , s = s − s − s ,s = s − s − s , s = s − s − s . (2.4)It is also useful to introduce the parity-odd invariant (cid:15) = 4i (cid:15) µνρσ p µ p ν p ρ p σ , (2.5)where (cid:15) µνρσ is the totally anti-symmetric Levi-Civita symbol. The invariant (cid:15) is relatedto the determinant of the Gram matrix G ij through( (cid:15) ) = ∆ ≡ det G ij = det (2 p i · p j ) , i, j ∈ { , . . . , } , (2.6)which reads explicitly∆ = ( s s + s s − s s + s s − s s ) − s s s ( s − s − s ) . (2.7)– 3 –s shown in Ref. [58], ∆ < s ij → s ij + i0 + gives Im( √ ∆) > . (2.8)In the physical scattering region of (2.1)–(2.2) the invariants are constrained to fulfil [44] s ≥ s , s − s ≥ s , ≥ s ≥ s − s , s − ≤ s ≤ s +15 , (2.9)with s ± = 1( s − s ) (cid:104) s s + s s ( s − s ) − s ( s s + s s + s s ) ± (cid:112) s s s s ( s + s − s )( s + s − s ) (cid:105) . (2.10) We start by considering the process in (2.1). Its UV-renormalised scattering amplitude canbe written as A aij = T aij A µνρ (cid:15) ∗ µ ( p ) (cid:15) ∗ ν ( p ) (cid:15) ∗ ρ ( p ) = T aij A , (3.1)where we find it convenient to factor out the SU(3) colour generator T aij , with i and j being the colour indices of the quark and anti-quark and a the colour index of the gluonin the adjoint representation. (cid:15) ∗ ( p ) is the polarisation vector of the outgoing gluon and (cid:15) ∗ ( p ), (cid:15) ∗ ( p ) are the ones of the photons. Clearly, for the qg channel in (2.2) the sameapplies upon suitably renaming the external momenta p ↔ p and ommitting the complexconjugate for the polarisation vector of the incoming gluon. The amplitude A stripped ofthe colour generator T a is then perturbatively expanded in the strong coupling constant α s as A = (4 πα ) Q q √ πα s (cid:18) A + (cid:16) α s π (cid:17) A + (cid:16) α s π (cid:17) A + O ( α ) (cid:19) , (3.2)with α the fine-structure constant and Q q the electric charge of the quarks in units ofelectron charge. The decomposition of (3.2) fixes the normalisation of the expansion terms A i . The amplitude squared and summed over colours and polarisations can be expressedas (cid:88) col , pol | A aij | = C (cid:18) W + (cid:16) α s π (cid:17) W ] + (cid:16) α s π (cid:17) ( W + 2Re [ W ]) + O ( α ) (cid:19) , (3.3)where we introduced the interference terms W ij = (cid:88) pol A ∗ i A j , (3.4)and a global factor C , that accounts for colour sums, which is defined as C = (4 π ) α α s ( N − Q q . (3.5)– 4 – igure 1 . Examples of two-loop diagrams which contribute to the various colour factors in (3.9). Note that all colour and helicity degrees of freedom are summed over in (3.3) rather thanperforming an average for the initial states. The perturbative corrections to the expansioncoefficients for the amplitudes A i , and to their interferences W ij , can be decomposed interms of their colour and electroweak factors which can be expressed as polynomials in C A = N, C F = N − N , n γf = 1 Q q n f (cid:88) i Q i , n γγf = 1 Q q n f (cid:88) i Q i , (3.6)where C A and C F are the quadratic Casimirs for the colour group SU( N ), ( N = 3 for QCD),and n f is the number of (massless) quarks running in closed fermionic loops. The factors n γf and n γγf correspond to the number of quarks, weighted with their electric charge, runningin closed fermion loops, with one and two external photons attached to, respectively. Thecoefficients of various powers of the factors in (3.6) constitute gauge-invariant subsets ofthe final result, thus it is natural to decompose the amplitude in terms of them.Within the prescription of (3.4)–(3.5), the LO result W is a rational function ofinvariants only, while the tree with one-loop interference admits the decomposition W = C A W (1)01 + C F W (2)01 + n γγf W (3)01 + n f W (4)01 . (3.7)The last contribution is introduced by the renormalisation procedure, see Section 5 fordetails. The tree-two loop interference has a richer structure, which we decide to organiseas W = (cid:88) i =1 c i W ( i )02 , (3.8)with the individual colour factors given by c = C A , c = C A C F , c = C F , c = C A n f ,c = C F n f , c = n γγf n f , c = C A n γγf , c = C F n γγf ,c = (8 C A − C F ) n γf , c = n f , (3.9)where c arises from a contraction of type d abc d abc and c follows again from ultraviolet(UV) renormalisation. A selection of representative diagrams contributing to each colour– 5 –tructure is shown in Fig. 1, where diagram ( i ) contributes to the factor c i and the colourfactor c follows from diagrams of type (9) as well. Eventually, one can then expand thepolynomial in N for c , , and write C A W (1)02 + C A C F W (2)02 + W (3)02 = N (cid:102) W (1)02 + (cid:102) W (2)02 + 1 N (cid:102) W (3)02 , (3.10)where (cid:102) W (1)02 is the leading colour (LC), (cid:102) W (2)02 is the next-to-leading colour (NLC) and (cid:102) W (3)02 is the next-to-next-to-leading colour (NNLC) contribution. The expression in terms ofCasimir operators can always be recovered through the set of identities N = C A , C A − C A C F , N = C A − C A C F + 4 C F . (3.11)In this paper we will focus on the calculation of the LC contribution (cid:102) W (1)02 and the threefermionic contributions W ( i )02 , i = 4 , ,
6, which can be obtained by considering only theplanar two-loop diagrams.
We perform our calculation in conventional dimensional regularisation (CDR) [59, 60],working in d = 4 − (cid:15) dimensions. UV and IR singularities will then manifest as poles inthe regulator (cid:15) .In order to compute the coefficients W , W and W , we start by producing allrelevant Feynman diagrams with QGRAF [61]. We find 6 diagrams at tree level, 80 di-agrams at one loop and 1716 diagrams at two loops. We use FORM [62] to manipulatethe diagrams, interfere them with their tree-level counterparts and perform the relevantLorentz and Dirac algebra. Starting at two loops, it is particularly convenient to groupthe diagrams depending on the graph they can be mapped to, before performing all sym-bolic manipulations. This allows us to avoid repeating expensive operations multiple times.After this step is complete, we can express our two-loop interfered amplitude in terms ofscalar Feynman integrals drawn from two different integral families. We define the integralsas I fam n ,...,n = e (cid:15)γ E (cid:90) (cid:89) i =1 (cid:18) d d k i iπ d/ (cid:19) D n ...D n , (4.1)where γ E ∼ . any of the colour factors (3.9). We find that after applying thesymmetrisations above, and modulo crossings of the external legs, the interference of thetwo-loop amplitude with the tree-level amplitude requires the reduction of 1811 integralsin family A and 508 integrals in family B , which include all the integrals required to– 6 –rop. den. Family A Family B D k k D ( k + p ) ( k − p ) D ( k + p + p ) ( k − p − p ) D ( k + p + p + p ) ( k − p − p − p ) D k k D ( k + p + p + p ) ( k − p − p − p − p ) D ( k + p + p + p + p ) ( k − k ) D ( k − k ) ( k − k + p ) D ( k + p + p + p + p ) ( k − p ) D ( k + p ) ( k − p − p ) D ( k + p + p ) ( k − p − p − p ) Table 1 . Definition of the two integral families used in the calculation. We note that also boththe non-planar hexagon-box and the double-pentagon topologies can be described by family B andcrossings thereof. compute both the hexagon-box and the double-pentagon topologies. In both families weneed at most rank-5 integrals (up to 5 irreducible scalar products in the numerator). Fromthis point on, we focus only on the planar colour factors and therefore only on integralsbelonging to family A . We performed their reduction using both Kira [35, 65, 66] and anin-house implementation, Finred, which employs finite field sampling, syzygy techniques,and denominator guessing, see [37, 67] for more details. It is worth stressing that the planarreductions up to rank 5 did not constitute an issue for either program here, and took e.g.40 hours on 36 cores with Kira and a similar runtime with Finred. The reduction lists forfamily A produced in this way are not extremely complicated, with a size of around 390MB in total. Note that we performed the reduction directly in terms of the pre-canonicalset of master integrals defined in [44]. We stress here that these lists do not include thecrossings necessary to reduce all diagrams.To arrive at a complete set of reduction identities and to render their inclusion assimple as possible, we proceed as follows. The integration-by-parts solvers deliver eachintegral coefficient as a rational function in a common-denominator representation. Wefind it useful to convert the rational functions to a partial fraction decomposed form.Due to the choice of master integrals, we encounter only irreducible denominator factorswhich depend on either d or the kinematic variables. For the d dependence of the rationalfunctions, a simple univariate partial fraction decomposition is sufficient. In contrast, thedecomposition involving the kinematic denominators is computationally non-trivial. We– 7 –ncounter the irreducible denominator factors { d , . . . , d } = { s , s , s , s , s , s + s − s , s − s − s , s + s − s ,s − s + s , s − s − s , s + s , s − s , s + s , s − s ,s − s , s + s − s − s , s + s , s − s , s − s ,s + s − s − s , s + s − s − s , s + s ,s − s − s + s , s − s − s − s , s + s } . (4.2)For them, we employ the MultivariateApart package [67] to perform a multivariate par-tial fraction decomposition; see also [68–70] for related decomposition techniques. Thedecomposition of [67] is based on replacing all irreducible denominator factors d k ( { s ij } ) ina given expression according to 1 d k ( { s ij } ) = q k , (4.3)and reducing the resulting polynomial with respect to the ideal generated by { d ( { s ij } ) q − , . . . , d ( { s ij } ) } (4.4)using a suitable monomial ordering. Due to the specific structure of our denominator list,the monomial block ordering of [67] coincides with a lexicographic ordering of the q i be-fore any s ij are considered, which produces a Le˘ınartas decomposition [71, 72]. Note thatwe included only the local denominators of each coefficient in the initial partial fractiondecomposition. This run took a couple of days and reduced the size of the reduction listby almost a factor 5, down to around 80 MB. In particular, the reduction in complex-ity is particularly pronounced for the most complicated rank-5 identities, for which up toa factor of 40 reduction in size is seen. While it has been known for a long time thatintegration-by-parts identities can become substantially simpler even using naive variantsof multivariate partial fraction decompositions, we would like to point out a systematicstudy of the impact of these new algorithms on the size of the reduction identities whichhas recently appeared in [70]. Starting from these substantially simpler identities, we per-form the relevant crossings and then a second partial fraction decomposition. The secondpartial fraction decomposition employs a global Gr¨obner basis for the ideal (4.4), takinginto account the denominator factors of all coefficients at once. We emphasize that thismethod allows to decompose terms of a sum locally, but still ensures a globally uniquerepresentation of the rational functions in the results. With this multi-step procedure,we obtained all identities and their crossings in a form suitable for insertion in the dia-grammatic calculation, while keeping the dimension of the expressions under control ateach step. After insertion of the reduction identities into the interference terms, a finalquick last partial fraction decomposition handles the remaining few additional denominatorfactors. In practice, we find the partial fraction decomposition of the original, uncrossedidentities to be the most time-consuming step, which nevertheless could be handled in acouple of days. The remaining steps to arrive at the reduced and fully partial fractionedinterference terms took a few hours. – 8 –n order to optimize our results for numerical stability, we found it useful to tune ourmonomial ordering for the partial fraction decomposition in such a way that spurious polesand unnecessarily high powers of the denominators are avoided. This choice of orderingwas performed for each colour factor and loop order separately. In this way, we obtaineda rather compact expression for the two-loop interference terms in terms of canonicalmaster integrals, with full dependence on the space-time dimensions d . As we will see inSection 6, by substituting the explicit results for the master integrals, expanding about d = 4, subtracting the poles, and properly simplifying the remaining rational functions, weare able to obtain extremely compact, fully analytic results for the finite remainders of theinterference terms. We work in fully massless QCD and keep the number of light fermion flavours n f generic;we do not consider any loop corrections due to massive quarks. We perform UV renormal-isation in the standard MS scheme, which allows one to express renormalised amplitudesin terms of unrenormalised ones by simply replacing the bare coupling constant α bs withthe running coupling α s ( µ ), evaluated at the scale µ , α bs µ (cid:15) S (cid:15) = α s µ (cid:15) (cid:20) − β (cid:15) (cid:16) α s π (cid:17) + (cid:18) β (cid:15) − β (cid:15) (cid:19) (cid:16) α s π (cid:17) + O ( α ) (cid:21) , (5.1)where S (cid:15) = (4 π ) (cid:15) e − (cid:15)γ E , and β and β denote the first two perturbative orders of the QCDbeta function, β = 116 C A − T r n f , β = 176 C A − C A T r n f − C F T r n f , (5.2)with T r = 1 /
2. The renormalised interference terms in (3.3) are explicitly obtained as W = S − (cid:15) W b − β (cid:15) W , W = S − (cid:15) W b − β (cid:15) S − (cid:15) W b − (cid:18) β (cid:15) − β (cid:15) (cid:19) W , (5.3)where with a superscript b we denote bare quantities. After UV renormalisation the resultsin (5.3) contain only IR poles. We subtract them according to Catani’s factorisationformula [73], which makes it possible to reorganise the interference terms as W = I ( (cid:15), µ ) W + W fin01 , W = I ( (cid:15), µ ) W + I ( (cid:15), µ ) W + W fin02 , (5.4)where the action of the operators I and I produce the complete IR structure of therenormalised amplitudes. The finite remainders W fin01 , W fin02 will be the main result of ourpaper. Before discussing these in detail, let us show the explicit formulae for the Catani– 9 –perators for the process under consideration. By direct calculation it is easy to see thatfor the q ¯ q channel I ( (cid:15), µ ) = e (cid:15)γ E − (cid:15) ) (cid:20) ( C A − C F ) (cid:18) (cid:15) + 32 (cid:15) (cid:19) (cid:18) − µ s (cid:19) (cid:15) − (cid:18) C A (cid:18) (cid:15) + 34 (cid:15) (cid:19) + γ g (cid:15) (cid:19) (cid:18)(cid:18) − µ s (cid:19) (cid:15) + (cid:18) − µ s (cid:19) (cid:15) (cid:19) (cid:21) , (5.5)where γ g = β , and I ( (cid:15), µ ) = − I ( (cid:15), µ ) (cid:18) I ( (cid:15), µ ) + 2 β (cid:15) (cid:19) +e − (cid:15)γ E Γ(1 − (cid:15) )Γ(1 − (cid:15) ) (cid:18) β (cid:15) + K (cid:19) I (2 (cid:15), µ ) + e − (cid:15)γ E (cid:15) Γ(1 − (cid:15) ) H , (5.6)where K is universal, K = (cid:18) − π (cid:19) C A − T r n f , (5.7)and H is a renormalisation and process-dependent factor, in our case given by H = 2 H ,q + H ,g , (5.8)with H ,q and H ,g which read explicitly [74] H ,g = (cid:18) ζ π (cid:19) C A + 2027 T r n f − (cid:18) π
36 + 5827 (cid:19) C A T r n f + C F T r n f , (5.9) H ,q = (cid:18) −
38 + π − ζ (cid:19) C F + (cid:18) − π
48 + 13 ζ (cid:19) C A C F + (cid:18) π − (cid:19) C F T r n f . (5.10)In order to obtain the corresponding expressions for the qg channel, it is enough to swapthe indices 2 ↔ W retaining exact d dependence and we expand W up to O ( (cid:15) ). This is necessary for thecorrect subtraction of the IR poles according to (5.4).After UV renormalisation and IR subtraction the one-loop finite remainder takes theform: W fin01 = C A W (1) , fin01 + C F W (2) , fin01 + T r n γγf W (3) , fin01 + T r n f W (4) , fin01 , (5.11)where, as anticipated in Section 3, the last term on the right-hand side originates from therenormalisation and subtraction procedures and it is proportional to the LO result. As forthe two-loop finite remainder, we find it convenient to cast it in the form W fin02 = (cid:88) i =1 ˜ c i (cid:102) W ( i ) , fin02 + (cid:88) i =4 c i W ( i ) , fin02 , (5.12) Note that the form of I and I is dictated entirely by the presence of a q ¯ q pair and one gluon ascoloured external states. – 10 –here the colour factors c ,..., are the same as defined in (3.9) and we introduced˜ c = N , ˜ c = 1 , ˜ c = N − . (5.13)The factors ˜ c , , clearly organise the colour tower, while c is introduced by UV renormal-isation and IR subtraction. As we have already argued at the end of Sec. 3, the coefficients (cid:102) W (1) , fin02 , W (4) , fin02 , W (5) , fin02 and W (6) , fin02 entail only planar two-loop corrections. These fourfinite contributions constitute the main object of this paper.A first test that we performed on our results at the analytical level, is that we reproduceall the (cid:15) poles of the renormalised interference terms, at one and two loops, as predictedby Catani’s factorisation formula in (5.4). The LO contribution defined in (3.3) for process (2.1) has a very compact expression,which, computed in CDR, reads W = 8 s s s s s (cid:0) s + s + (cid:15) (2 − (cid:15) )( s + s + s + s )+2 (cid:15) (1 − (cid:15) ) s s + (cid:15) (3 + (cid:15) )( s − s )( s − s ) (cid:1) + (cid:0) → , → , → (cid:1) + (cid:0) → , → , → (cid:1) , (6.1)where we kept full (cid:15) -dependence and in the last row we denote cyclic permutations of theindices labelling the kinematic invariants. An equivalent expression holds for process (2.2)upon replacing 2 ↔ O ( (cid:15) ), contain only simple logarithms and dilogarithms andthey are free of the parity odd-invariant (cid:15) . At the two-loop order, as well as for the higher (cid:15) powers of the one loop result, the class of transcendental functions needs to be extendedbeyond classical polylogarithms and (cid:15) enters explicitly. We employ the representation interms of so-called Pentagon Functions presented in [47] . Thus we write a generic finiteremainder as W ( m ) , fin0 n = (cid:88) k r k ( { s ij , (cid:15) } ) f k ( { s ij , (cid:15) } ) , (6.2)where m denotes some given colour factor, n = 1 , r k is a rational function of the kinematicinvariants and f k indicates a pentagon function. Equation (6.2) holds formally at LO, n = 0, as well, upon putting f k = 1. We remind the reader that at this step we representeach rational function r k in its partial fractions decomposed form. We further stress that For further details about the definition and implementation of the pentagon functions see [47]. – 11 –he algorithm presented in [67] avoids the presence of spurious denominators, therefore ourfinal expressions contain all and only the physical denominators of the amplitudes.Despite that the partial fraction decomposition allows us to obtain rather compactresults already at this level, we find that the representation in (6.2) is not the most compactrepresentation yet. This is due to the fact that the rational functions r k are not linearlyindependent and, in fact, the set of independent rational functions is substantially smaller.The minimal set of rational functions can be obtained as follows. First of all, in order toguarantee a globally unique representation of all rational functions, we employ a partialfraction decomposition with respect to a common Gr¨obner basis. As a result, each rationalfunction r k is represented as a polynomial of maximum degree p in the 5 kinematic variables s ij (2.3) and the 25 inverse denominators q i (4.3), r k = (cid:88) m + ... + m n ≤ p a k,m ...m M m ...m , where M m ...m ≡ q m · · · q m s m · · · s m (6.3)and the coefficients a k,m ...m are rational numbers. We emphasize that only certain M k,m ...m occur as irreducible monomials in our partial fractioned results; in practiceit is convenient to enumerate only those. We are interested in determining linear relationsbetween different rational functions, 0 = (cid:88) k b k r k . (6.4)and employing them to re-express all r k in terms of a linearly independent subset of them.By inserting the partial fractioned form (6.3) and observing that the monomials M m ...m are independent, we obtain a system of linear relations for the b k ,0 = (cid:88) k a k,m ...m b k (6.5)for each of the irreducible monomials M m ...m . In our case, there are many more monomi-als than rational functions, and the system is over-constrained. Nevertheless, many of theequations turn out to be linearly dependent, allowing us to find a solution to the systemin terms of a basis of independent rational functions. Note that this can be achieved by arow reduction of a matrix of rational numbers, which can easily be done e.g. with a finitefield solver like Finred . Similarly to what happens for the reduction to master integrals,there is not a unique solution for the basis of rational functions and a more natural choicecan lead to more compact results for the scattering amplitude. We find that very compactanalytic expressions can be found by simplifying the rational functions for each loop orderand colour factor separately. In practice, we order the rational functions according to thenumber of monomials they contain, and at each step we remove the one with the largestnumber.This approach makes it possible to drastically reduce the size of the final results. As anexample, for the LC coefficient (cid:102) W (1) , fin02 we start with an expression of around 64 MB which,after moving to a basis of independent rational functions, can be reduced to around 3 MB.This makes it possible not only to have a more compact result, but, most importantly, a– 12 –ore efficient and potentially more stable numerical implementation. Clearly, we cannotdemonstrate that a more compact representation does not exist for a different choice ofindependent rational functions. On the contrary, we expect that more compact expressionscould be obtained if one does not insist on using a set of independent kinematics invariants,see for example Eq. (6.1) and Ref. [22]. Nevertheless, our current representation is morethan suitable for practical use.As last a remark, we stress that the procedure above does not produce the minimalnumber of rational functions. Indeed, starting from the results simplified according to theprocedure above, we can attempt to further relate rational functions across different colourstructures and different loop orders. In this way, more relations can be found and a minimalset of rational functions can be identified. We find this second representation particularlysuitable for implementation of the interference terms in a numerical code, as described inthe following section. We provide analytical expressions for various colour factors at thedifferent loop orders with the arXiv submission of this manuscript. We implemented our results in the
C++ numerical code aajamp , which we distribute throughthe git repository [75]; it can be easily downloaded with git clone [email protected]:vmante/aajamp.git
Details on the installation procedure and usage of the package are given in the git repository.The main purpose of the code is to evaluate the finite remainders W , all W ( i ) , fin01 of (5.11)and W ( i ) , fin02 of (5.12) for i ∈ { , , , , } . The evaluation of all transcendental functionsin our code, both at one and two-loop level, is handled by the PentagonFunctions-cpp library of [47, 76, 77]. In order to implement our formulae in an efficient fashion we adoptedthe optimised code generation provided by
FORM [78].We then checked the implementation of our LO and NLO results against
OpenLoops 2 [79] and find agreement.In Table 2 we provide benchmark results for a kinematic point in the physical regiondefined by s = 157 , s = − , s = 83 , s = 61 , s = − , µ = 100 . (7.1)To assess the performance of our code, we measured the evaluation time of the NLOand NNLO results in double precision on a single Intel i7-9750H CPU @ 2.60GHz coreusing gcc . × − ms and 1.2 s per phasespace point for the NLO and NNLO contributions, respectively. We note that most of theevaluation time for a given point is spent on the computation of the pentagon functions.The aajamp package allows for numerical evaluations in double and quadruple floating-point precision, i.e. 16 and 32 decimal digits representation. The quadruple precisionevaluation is adapted to follow the strategy of the PentagonFunctions-cpp library, which– 13 – ( i )0 n q ¯ q → gγγ qg → qγγ W . . W (1) , fin01 . − . W (2) , fin01 − . − . W (3) , fin01 . − . W (4) , fin01 − . − . W (1) , fin02 . − . W (4) , fin02 − . . W (5) , fin02 . . W (6) , fin02 − . . W (10) , fin02 . − . Table 2 . Benchmark results for the tree-level, one-loop and two-loop interference terms for thepartonic channels q ¯ q and qg for the kinematic point in (7.1). Only the real part of W ( i )0 n is presented. utilizes the qd library [80] for higher-precision arithmetic. Therefore, aajamp relies on qd as well. We note that the PentagonFunctions-cpp and qd libraries offer octuple precision,but we reckon that quadruple precision is sufficient for phenomenological applications.Although the user can choose at the interface level between a purely double or quadru-ple arithmetic evaluation, we have implemented a primitive precision control system. Thisis motivated by the fact that individual terms appearing in the squared matrix elements candevelop spurious singularities inside the physical phase space, leading to more or less severenumerical instabilities. These are typically associated with small Gram determinants, orin our case, with denominators in the rational functions which become small comparedto the typical energy scale of the process. We observe that in our final expressions, theparity odd invariant (cid:15) does not appear in any denominator and, correspondingly, smallGram determinants are not of prior concern. Instead, we focus on small denominators inthe rational functions, which are not associated with a physical singularity, but lead to amajor source of numerical instabilities in the evaluations. Let us stress that, in principle,also other types of numerical instabilities could arise, but our denominator based analysisindeed works well in practice, as we will show. Examples of these spurious singularitiesarise in events with all particle having a large transverse momentum and entailing collinearphoton pairs or collinear gluon-photon pairs. We therefore find it natural to activate aquadruple precision evaluation if any of the 25 denominators becomes smaller than a giventhreshold which the user can tune. As can be seen in Fig. 2, this precision-control systemsignificantly improves the reliability of the numerical evaluations. Here, we assess the levelof numerical stability via a rescaling test, i.e. we consider the quantity D defined as D = log (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ξ n W ( i )02 ( ξs ij , ξµ ) W ( i )02 ( s ij , µ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (7.2) Strictly speaking the 32 decimal digits representation of qd is based on a double-double arithmetic. – 14 – − − − − − − − accuracy D min − − − − − f r a c t i o n o f e v e n t s q ¯ q → gγγ at O ( α ) f W (1) , fin02 dp W (4) , fin02 dp f W (1) , fin02 rescue W (4) , fin02 rescue Figure 2 . Probability of finding an event with an instability level D defined in (7.2) larger than D min for the two colour factors W (1) , fin02 and W (4) , fin02 . The labels “dp” and “rescue” are describedin the text. Events were generated according to (7.3) where we perform two evaluations of the colour factor W ( i ) , fin02 with input kinematics ( s ij , µ )and with the same kinematics rescaled by a factor ξ . The extra factor ξ n accounts for themass dimensionality of W ( i ) , fin lm , and in our case we just have n = 1. The definition of D in (7.2) intuitively provides an indicator for the number of digits one can trust for a givencomputation, thus we identify it as our level of instability. In Fig. 2 we show a cumulativehistogram of phase-space points, which provide an evaluation of W ( i ) , fin02 with an instability D larger than D min . We plot two selected colour factors, both evaluated in pure doubleprecision, labelled as “dp”, and with the precision control system activated, labelled as“rescue”. We generated 10 uniformly distributed events with Rambo [81] subject to theconstraints E com = 1 TeV , p T ,g >
30 GeV , p T ,γ >
30 GeV , p T ,γ >
30 GeV , (7.3)where E com is the energy in rest frame of the colliding partons and p T ,i is the transversemomentum of particle i . One can see that for this ensemble of events, the precision controlsystem is able to capture and cure the worst instabilities, effectively guaranteeing a verygood level of accuracy. Once again, we stress that a more sophisticated stability systemcould be devised in principle, especially in regions of soft or collinear emissions. In this paper, we presented the calculation of the leading colour and light fermionic two-loop corrections for the production of two photons and a jet in quark-antiquark annihilationand quark-gluon scattering. This calculation has been made possible by the combinationof state-of-the-art techniques for the reduction of multiloop Feynman integrals, new ideas– 15 –bout their representations in terms of multivariate partial fractions, and recent results forthe relevant master integrals. In particular, we have shown how the spin summed inter-ference between the two-loop and the tree amplitude can be computed from its Feynmandiagrammatic representation, resulting in very compact analytic expressions. We havechecked that our two-loop corrections display the correct pole structure, as first predictedby Catani for QCD amplitudes at this perturbative order.In order to demonstrate the flexibility and usability of our result, we have also imple-mented the tree-level, one-loop and two-loop finite remainders in a
C++ library, providingall (cid:15) expansions through to transcendental weight four. Our library links against the
PentagonFunctions-cpp library and allows the user to evaluate the loop corrections indouble and quadruple precision. We have also introduced a simple precision control sys-tem that allows the code to identify phase space points prone to loss of precision, such thatquadruple precision evaluations are restricted to a minimum. We envisage that the algo-rithms developed for this calculation can be extended to solve future cutting-edge problemsin the computation of multiloop multileg scattering amplitudes relevant for collider physicsphenomenology.
Acknowledgments
We thank Fabrizio Caola and Tiziano Peraro for various discussions on different aspectsof the calculation. We further thank Fabrizio Caola for comments on the manuscript. FBis thankful to Jean-Nicolas Lang for useful suggestions on several aspects of the numericalcode implementation. The research of FB is supported by the ERC Starting Grant 804394HipQCD. BA and AvM are supported in part by the National Science Foundation throughGrant 2013859. LT is supported by the Royal Society through Grant URF/R1/191125.
References [1]
ATLAS collaboration,
Search for resonances in diphoton events at √ s =13 TeV with theATLAS detector , JHEP (2016) 001 [ ].[2] CMS collaboration,
Search for physics beyond the standard model in high-mass diphotonevents from proton-proton collisions at √ s =
13 TeV , Phys. Rev. D (2018) 092001[ ].[3] V. Del Duca, F. Maltoni, Z. Nagy and Z. Trocsanyi, QCD radiative corrections to promptdiphoton production in association with a jet at hadron colliders , JHEP (2003) 059[ hep-ph/0303012 ].[4] T. Gehrmann, N. Greiner and G. Heinrich, Photon isolation effects at NLO in gammagamma + jet final states in hadronic collisions , JHEP (2013) 058 [ ].[5] A. Gehrmann-De Ridder, T. Gehrmann and E. Glover,
Antenna subtraction at NNLO , JHEP (2005) 056 [ hep-ph/0505111 ].[6] A. Daleo, T. Gehrmann and D. Maitre, Antenna subtraction with hadronic initial states , JHEP (2007) 016 [ hep-ph/0612257 ]. – 16 –
7] M. Czakon,
A novel subtraction scheme for double-real radiation at NNLO , Phys. Lett. B (2010) 259 [ ].[8] M. Czakon and D. Heymes,
Four-dimensional formulation of the sector-improved residuesubtraction scheme , Nucl. Phys. B (2014) 152 [ ].[9] J. Gaunt, M. Stahlhofen, F. J. Tackmann and J. R. Walsh,
N-jettiness Subtractions forNNLO QCD Calculations , JHEP (2015) 058 [ ].[10] R. Boughezal, C. Focke, W. Giele, X. Liu and F. Petriello, Higgs boson production inassociation with a jet at NNLO using jettiness subtraction , Phys. Lett. B (2015) 5[ ].[11] F. Caola, K. Melnikov and R. R¨ontsch,
Nested soft-collinear subtractions in NNLO QCDcomputations , Eur. Phys. J. C (2017) 248 [ ].[12] S. Badger, C. Broennum-Hansen, H. B. Hartanto and T. Peraro, First look at two-loopfive-gluon scattering in QCD , Phys. Rev. Lett. (2018) 092001 [ ].[13] S. Abreu, F. Febres Cordero, H. Ita, B. Page and M. Zeng,
Planar Two-Loop Five-GluonAmplitudes from Numerical Unitarity , Phys. Rev.
D97 (2018) 116014 [ ].[14] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita and B. Page,
Analytic Form of PlanarTwo-Loop Five-Gluon Scattering Amplitudes in QCD , Phys. Rev. Lett. (2019) 082002[ ].[15] S. Abreu, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov,
Planar Two-Loop Five-PartonAmplitudes from Numerical Unitarity , JHEP (2018) 116 [ ].[16] D. Chicherin, T. Gehrmann, J. M. Henn, P. Wasser, Y. Zhang and S. Zoia, Analytic resultfor a two-loop five-particle amplitude , Phys. Rev. Lett. (2019) 121602 [ ].[17] S. Badger, C. Brønnum-Hansen, H. B. Hartanto and T. Peraro,
Analytic helicity amplitudesfor two-loop five-gluon scattering: the single-minus case , JHEP (2019) 186 [ ].[18] S. Badger, D. Chicherin, T. Gehrmann, G. Heinrich, J. Henn, T. Peraro, P. Wasser et al., Analytic form of the full two-loop five-gluon all-plus helicity amplitude , Phys. Rev. Lett. (2019) 071601 [ ].[19] D. Chicherin, T. Gehrmann, J. M. Henn, P. Wasser, Y. Zhang and S. Zoia,
The two-loopfive-particle amplitude in N = 8 supergravity , JHEP (2019) 115 [ ].[20] H. B. Hartanto, S. Badger, C. Brønnum-Hansen and T. Peraro, A numerical evaluation ofplanar two-loop helicity amplitudes for a W-boson plus four partons , JHEP (2019) 119[ ].[21] S. Abreu, B. Page, E. Pascual and V. Sotnikov, Leading-Color Two-Loop QCD Correctionsfor Three-Photon Production at Hadron Colliders , .[22] G. De Laurentis and D. Maˆıtre, Two-Loop Five-Parton Leading-Colour Finite Remainders inthe Spinor-Helicity Formalism , .[23] H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, Two-loop leading-color helicityamplitudes for three-photon production at the LHC , .[24] H. A. Chawdhry, M. L. Czakon, A. Mitov and R. Poncelet, NNLO QCD corrections tothree-photon production at the LHC , JHEP (2020) 057 [ ].[25] S. Kallweit, V. Sotnikov and M. Wiesemann, Triphoton production at hadron colliders inNNLO QCD , . – 17 –
26] F. Tkachov,
A Theorem on Analytical Calculability of Four Loop Renormalization GroupFunctions , Phys.Lett.
B100 (1981) 65.[27] K. Chetyrkin and F. Tkachov,
Integration by Parts: The Algorithm to Calculate betaFunctions in 4 Loops , Nucl.Phys.
B192 (1981) 159.[28] S. Laporta,
High precision calculation of multiloop Feynman integrals by difference equations , Int.J.Mod.Phys.
A15 (2000) 5087 [ hep-ph/0102033 ].[29] J. Gluza, K. Kajda and D. A. Kosower,
Towards a Basis for Planar Two-Loop Integrals , Phys. Rev. D (2011) 045012 [ ].[30] R. M. Schabinger, A New Algorithm For The Generation Of Unitarity-CompatibleIntegration By Parts Relations , JHEP (2012) 077 [ ].[31] H. Ita, Two-loop Integrand Decomposition into Master Integrals and Surface Terms , Phys.Rev.
D94 (2016) 116015 [ ].[32] K. J. Larsen and Y. Zhang,
Integration-by-parts reductions from unitarity cuts and algebraicgeometry , Phys. Rev.
D93 (2016) 041701 [ ].[33] J. B¨ohm, A. Georgoudis, K. J. Larsen, M. Schulze and Y. Zhang,
Complete sets oflogarithmic vector fields for integration-by-parts identities of Feynman integrals , Phys. Rev.D (2018) 025023 [ ].[34] J. B¨ohm, A. Georgoudis, K. J. Larsen, H. Sch¨onemann and Y. Zhang, Completeintegration-by-parts reductions of the non-planar hexagon-box via module intersections , JHEP (2018) 024 [ ].[35] P. Maierh¨ofer, J. Usovitsch and P. Uwer, Kira—A Feynman integral reduction program , Comput. Phys. Commun. (2018) 99 [ ].[36] J. Klappert and F. Lange,
Reconstructing rational functions with FireFly , Comput. Phys.Commun. (2020) 106951 [ ].[37] B. Agarwal, S. P. Jones and A. von Manteuffel,
Two-loop helicity amplitudes for gg → ZZ with full top-quark mass effects , .[38] A. von Manteuffel and R. M. Schabinger, A novel approach to integration by parts reduction , Phys. Lett.
B744 (2015) 101 [ ].[39] A. von Manteuffel and R. M. Schabinger,
Quark and gluon form factors to four-loop order inQCD: the N f contributions , Phys. Rev.
D95 (2017) 034030 [ ].[40] T. Peraro,
Scattering amplitudes over finite fields and multivariate functional reconstruction , JHEP (2016) 030 [ ].[41] T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflowgraphs , JHEP (2019) 031 [ ].[42] A. V. Smirnov and F. S. Chuharev, FIRE6: Feynman Integral REduction with ModularArithmetic , Comput. Phys. Commun. (2020) 106877 [ ].[43] C. G. Papadopoulos, D. Tommasini and C. Wever,
The Pentabox Master Integrals with theSimplified Differential Equations approach , JHEP (2016) 078 [ ].[44] T. Gehrmann, J. Henn and N. Lo Presti, Pentagon functions for massless planar scatteringamplitudes , JHEP (2018) 103 [ ]. – 18 –
45] D. Chicherin, T. Gehrmann, J. Henn, N. Lo Presti, V. Mitev and P. Wasser,
Analytic resultfor the nonplanar hexa-box integrals , JHEP (2019) 042 [ ].[46] D. Chicherin, T. Gehrmann, J. Henn, P. Wasser, Y. Zhang and S. Zoia, All Master Integralsfor Three-Jet Production at Next-to-Next-to-Leading Order , Phys. Rev. Lett. (2019)041603 [ ].[47] D. Chicherin and V. Sotnikov,
Pentagon Functions for Scattering of Five Massless Particles , JHEP (2020) 167 [ ].[48] Z. Bern, L. J. Dixon and D. A. Kosower, Dimensionally regulated pentagon integrals , Nucl.Phys.
B412 (1994) 751 [ hep-ph/9306240 ].[49] A. Kotikov,
Differential equations method: New technique for massive Feynman diagramscalculation , Phys.Lett.
B254 (1991) 158.[50] E. Remiddi,
Differential equations for Feynman graph amplitudes , Nuovo Cim.
A110 (1997)1435 [ hep-th/9711188 ].[51] T. Gehrmann and E. Remiddi,
Differential equations for two loop four point functions , Nucl.Phys.
B580 (2000) 485 [ hep-ph/9912329 ].[52] C. G. Papadopoulos,
Simplified differential equations approach for Master Integrals , JHEP (2014) 088 [ ].[53] A. Primo and L. Tancredi, On the maximal cut of Feynman integrals and the solution oftheir differential equations , Nucl. Phys.
B916 (2017) 94 [ ].[54] A. Kotikov,
The Property of maximal transcendentality in the N=4 SupersymmetricYang-Mills , .[55] J. M. Henn, Multiloop integrals in dimensional regularization made simple , Phys.Rev.Lett. (2013) 251601 [ ].[56] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, M. Kraus, B. Page, E. Pascual et al.,
Caravel: A C++ Framework for the Computation of Multi-Loop Amplitudes with NumericalUnitarity , .[57] F. Caola, A. von Manteuffel and L. Tancredi, Di-photon amplitudes in three-loop QuantumChromodynamics , .[58] N. Byers and C. Yang, Physical Regions in Invariant Variables for n Particles and thePhase-Space Volume Element , Rev. Mod. Phys. (1964) 595.[59] G. ’t Hooft and M. J. G. Veltman, Regularization and Renormalization of Gauge Fields , Nucl. Phys.
B44 (1972) 189.[60] C. Bollini and J. Giambiagi,
Dimensional Renormalization: The Number of Dimensions as aRegularizing Parameter , Nuovo Cim.
B12 (1972) 20.[61] P. Nogueira,
Automatic Feynman graph generation , J.Comput.Phys. (1993) 279.[62] J. Vermaseren,
New features of FORM , math-ph/0010025 .[63] A. von Manteuffel and C. Studerus, Reduze 2 - Distributed Feynman Integral Reduction , .[64] C. Studerus, Reduze-Feynman Integral Reduction in C++ , Comput.Phys.Commun. (2010) 1293 [ ].[65] P. Maierh¨ofer and J. Usovitsch,
Kira 1.2 Release Notes , . – 19 –
66] J. Klappert, F. Lange, P. Maierh¨ofer and J. Usovitsch,
Integral Reduction with Kira 2.0 andFinite Field Methods , .[67] M. Heller and A. von Manteuffel, MultivariateApart: Generalized Partial Fractions , .[68] A. Pak, The Toolbox of modern multi-loop calculations: novel analytic and semi-analytictechniques , J. Phys. Conf. Ser. (2012) 012049 [ ].[69] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov,
Analytic Form ofthe Planar Two-Loop Five-Parton Scattering Amplitudes in QCD , JHEP (2019) 084[ ].[70] J. Boehm, M. Wittmann, Z. Wu, Y. Xu and Y. Zhang, IBP reduction coefficients madesimple , JHEP (2020) 054 [ ].[71] E. K. Le˘ınartas, Factorization of rational functions of several variables into partial fractions , Soviet Math. (Iz. VUZ) (1978) 35.[72] A. Raichev, Le˘ınartas’ partial fraction decomposition , .[73] S. Catani, The Singular behavior of QCD amplitudes at two loop order , Phys. Lett. B (1998) 161 [ hep-ph/9802439 ].[74] L. Garland, T. Gehrmann, E. Glover, A. Koukoutsakis and E. Remiddi,
Two loop QCDhelicity amplitudes for e+ e- — > three jets , Nucl. Phys. B (2002) 227[ hep-ph/0206067 ].[75] https://gitlab.msu.edu/vmante/aajamp .[76] https://gitlab.com/pentagon-functions/PentagonFunctions-cpp .[77] https://gitlab.com/VasilySotnikov/Li2pp .[78] J. Kuipers, T. Ueda and J. A. M. Vermaseren,
Code Optimization in FORM , Comput. Phys.Commun. (2015) 1 [ ].[79] F. Buccioni, J.-N. Lang, J. M. Lindert, P. Maierh¨ofer, S. Pozzorini, H. Zhang and M. F.Zoller,
OpenLoops 2 , Eur. Phys. J. C (2019) 866 [ ].[80] Y. Hida, X. S. Li and D. H. Bailey, “Quad-double arithmetic: Algorithms, implementation,and application.” , 2000.[81] R. Kleiss, W. Stirling and S. Ellis, A new monte carlo treatment of multiparticle phase spaceat high energies , Computer Physics Communications (1986) 359 .(1986) 359 .