N 3 LO Higgs and Drell-Yan production at threshold: the one-loop two-emission contribution
Ye Li, Andreas von Manteuffel, Robert M. Schabinger, Hua Xing Zhu
aa r X i v : . [ h e p - ph ] D ec MITP/14-029, SLAC-PUB-15946 N LO Higgs and Drell-Yan production at threshold:the one-loop two-emission contribution
Ye Li, ∗ Andreas von Manteuffel, † Robert M. Schabinger, ‡ and Hua Xing Zhu § SLAC National Accelerator Laboratory,Stanford University, Stanford, CA 94309, USA The PRISMA Cluster of Excellence and Mainz Institute of Theoretical Physics,Johannes Gutenberg Universit¨at, 55099 Mainz, Deutschland
Abstract
In this paper, we study phenomenologically interesting soft radiation distributions in masslessQCD. Specifically, we consider the emission of two soft partons off of a pair of light-like Wilsonlines, in either the fundamental or the adjoint representation, at next-to-leading order. Our resultsare an essential component of the next-to-next-to-next-to-leading order threshold corrections toboth Higgs boson production in the gluon fusion channel and Drell-Yan lepton production. Ourcalculations are consistent with the recently published results for Higgs boson production. As anon-trivial cross-check on our analysis, we rederive a recent prediction for the Drell-Yan thresholdcross section using a completely different strategy. Our results are compact, valid to all orders inthe dimensional regularization parameter, and expressed in terms of pure functions. ∗ [email protected] † manteuff[email protected] ‡ [email protected] § [email protected] K -factors and theoretical uncertainties inherent in the cal-culation of such hadron collider observables [2, 4–11]. Unfortunately, this task is also highlynon-trivial; at next-to-next-to-next-to-leading order (N LO), the current state-of-the-art,one must compute three-loop virtual corrections, single-emission, two-loop real-virtual cor-rections, double-emission, one-loop real-virtual corrections, triple-emission real corrections,and collinear counterterms for the parton distribution functions. All of these ingredientshave severe infrared and/or collinear divergences which complicate their calculation andIRC finiteness is only achieved once all contributions have been appropriately combined.While full fixed-order calculations at N LO remain a challenge, significant progress hasbeen made over the last few years. For gluon fusion Higgs production, the Wilson coefficientsobtained by integrating out top quark loops to construct an effective ggh
Lagrangian areknown to three loops [13–15] and the relevant ultraviolet renormalization constants havebeen worked out in Refs. [16–19]. The collinear counter terms at three loops are also knownand have been obtained by the authors of Refs. [20–22]. The three-loop virtual correctionsare under control, as they can be extracted from the three-loop quark and gluon form factorscomputed in Refs. [23–26]. For finite top quark mass, an estimate of the N LO cross sectionfor Higgs production was made in Ref. [27] by arguing that one can exploit informationencoded in the soft gluon and high energy limits to construct an accurate approximation tothe full result.In recent years, rapid progress has been made towards the analytical calculation of therelevant phase space integrals by expanding them in the threshold limit . Two-loop, single-soft currents for Higgs and Drell-Yan production were computed to O ( ǫ ) in Refs. [29, 30], to O ( ǫ ) in Ref. [31], and, finally, to all orders in ǫ in Ref. [32], where ǫ is the usual parameterof dimensional regularization, D = 4 − ǫ . The contribution from the square of the one-loopamplitude for Higgs boson production in association with a single soft parton is also knownto all orders in ǫ [33–36]. Last year, the phase space integrals for soft triple-emission werecalculated in Ref. [37] to sufficiently high orders in ǫ for an N LO calculation, and veryrecently, the phase space integrals for Higgs boson production in association with two softpartons were computed by the same group, thereby completing the threshold corrections toHiggs boson production through gluon fusion at N LO [28]. However, this final piece of thesoft part of the three-loop bare threshold cross section has not yet appeared in a separatepublication. In this paper, we report on a parallel calculation of the eikonal double-emissionphase space integrals relevant to the N LO threshold calculation of gluon fusion Higgs orDrell-Yan lepton production. Our results are remarkably compact and valid to all ordersin the dimensional regularization parameter. The higher order in ǫ terms will be usefulcomponents of future N LO calculations, for example higher order extractions of quark andgluon collinear anomalous dimensions. It should be noted that, as pointed out in Ref. [28], the evaluation of the leading term in the thresholdexpansion by itself is not completely satisfactory. It may be necessary in future work to compute subleadingterms in the threshold expansion as well to determine the size and impact of power corrections.
2e describe our calculational strategy below. At a hadron collider with incoming hadrons N and N and center-of-mass energy √ s , the inclusive cross section for producing a Higgsboson of mass M H or Drell-Yan lepton pair of invariant mass M DY can be written as σ i (cid:0) s, M i (cid:1) = X a,b Z d x d x f a/N (cid:0) x , µ F (cid:1) f b/N (cid:0) x , µ F (cid:1) Z d z δ (cid:18) z − M i ˆ s (cid:19) ˆ σ iab (ˆ s, z, α s ( µ R )) , (1)where i = DY or H and the summation is over all partonic channels. To save space, we willoften use the abbreviations DY and H when discussing the processes that we consider in thispaper. The parton distribution function of parton n in hadron N , f n/N ( x, µ F ), depends onthe momentum fraction, x , and the factorization scale, µ F . ˆ σ iab (ˆ s, z, α s ( µ R )), the partoniccross section, depends on the partonic center-of-mass energy, √ ˆ s = √ x x s , the parameter z defined by Eq. (1), and the renormalized strong coupling constant at renormalizationscale µ R , α s ( µ R ). In this work, we take µ F = µ R = M i for the sake of simplicity. In fullfixed-order calculations, the dependence of partonic cross sections on ˆ s and α s ( M i ) is simplebut their dependence on z is quite non-trivial. Throughout this paper, we shall be concernedonly with a threshold expansion in the z → σ iab (ˆ s, z, α s ( M i )). In this regime,additional partons radiated into the final state are constrained to be soft by virtue of thephase space constraint ˆ s ≈ M i and the partonic cross sections for both H and DY simplifydramatically: ˆ σ H gg (ˆ s, z → , α s ( M H )) = σ gg ( M H ) G H ( z ) (2)ˆ σ DY q ¯ q (ˆ s, z → , α s ( M DY )) = σ q ¯ q ( M DY ) G DY ( z ) . (3)For simplicity, we keep the running coupling dependence of the functions introduced onthe right-hand side of Eqs. (2) and (3) implicit. In D = 4 − ǫ dimensions, σ gg ( M H ) and σ q ¯ q ( M DY ) are given respectively by σ gg ( M H ) = πλ ( α s ( M H ))8(1 − ǫ )( N c −
1) and σ q ¯ q ( M DY ) = 4 π (1 − ǫ ) α ( M DY ) e q N c M , (4)where λ ( α s ( M i )) is the effective coupling of the Higgs boson to gluons in the limit of infinitetop quark mass [39], L int = − λ ( α s ( M H )) HG µν, a G aµν , (5) N c is the number of colors, α ( M DY ) is the fine structure constant at renormalization scale M DY , and e q is the electric charge of quark q . It is worth pointing out that λ ( α s ( M H )) hasa perturbative expansion in the renormalized strong coupling constant of its own (see Ref.[15] for explicit expressions in a slightly different normalization) and, throughout this work,we employ the conventional dimensional regularization scheme. The coefficient functions G H ( z ) and G DY ( z ) contain singular distributions of the form δ (1 − z ) and (cid:20) ln k (1 − z )1 − z (cid:21) + , (6) We follow Moch and Vogt [38] and write Eq. (4) with off-shell photons in mind. g ( z )] + acts on functions regular in the z → Z d z f ( z ) [ g ( z )] + = Z d z (cid:16) f ( z ) − f (1) (cid:17) g ( z ) . (7)The functional form of the coefficient functions can be obtained by taking the z → G i ( z ) factorize into products of the form [44–46] G i ( z ) = H i ¯ S i ( z ) (8)where the H i are hard functions encoding the virtual corrections at threshold and the ¯ S i ( z )are renormalized soft functions encoding the real radiative corrections at threshold. At first,it might seem peculiar that we have not taken the soft functions in Eq. (8) to be a functionsof ω/M i , where ω is twice the energy of the soft QCD radiation in the final state. However,in threshold kinematics, we have ω = (1 − z ) √ ˆ s ≈ (1 − z ) M i (9)and we see that the form of Eq. (8) follows naturally from Eq. (9).In general, hard functions in soft-collinear effective theory are complex squares of hardWilson coefficients obtained by matching QCD onto soft-collinear effective theory. The per-turbative expansions of the Wilson coefficients relevant to Higgs and Drell-Yan productionare known to the order required for an N LO calculation of the G i ( z ). They can be extractedfrom the three-loop quark and gluon form factors [23–26]. In fact, they are written downexplicitly in Ref. [25], Eqs. (7.4) - (7.9). The actual hard coefficients that we need for ouranalysis are derived by taking the complex square of Eq. (7.3) in that work, setting µ = M i ,and then expanding in α s ( M i ) to third order.The bare soft functions can be written as squares of time-ordered matrix elements ofpairs of semi-infinite Wilson line operators, S i ( z ) = M i d i X X s h | T n Y n Y † ¯ n o δ (cid:16) λ − ˆ P (cid:17) | X s ih X s | T n Y ¯ n Y † n o | i . (10)As mentioned above, the soft functions defined in Eq. (10) depend on the ratio of twice theenergy of the soft QCD radiation to M i ( i.e. on 1 − z ), as well as the bare strong couplingconstant α s . d i is a conventional normalization constant, d H = N c − d DY = N c , (11)depending on whether the soft Wilson lines are in the adjoint or the fundamental represen-tation of su ( N c ). The summation in Eq. (10) is over all possible soft parton final states, | X s i . The operator ˆ P acts on the final state | X s i according toˆ P | X s i = 2 E X s | X s i , (12)where E X s is the energy of the soft radiation in final state | X s i . The Wilson line operators, Y n and Y † ¯ n , are respectively defined as in-coming, path-ordered ( P ) and anti-path-ordered (cid:0) P (cid:1) exponentials [47], Y n ( x ) = P exp (cid:18) ig Z −∞ d s n · A ( ns + x ) (cid:19) † ¯ n ( x ) = P exp (cid:18) − ig Z −∞ d s ¯ n · A (¯ ns + x ) (cid:19) . (13)In the above, A µ = A aµ T a , where the T a are either adjoint (for H) or fundamental (for DY) su ( N c ) matrices. As usual, n and ¯ n are light-like vectors whose space-like components areback-to-back and determine the beam axis. For a generic four-vector, k µ , we have n · k = k + and ¯ n · k = k − (14)given the usual definitions k + = k + k and k − = k − k .In the perturbative regime, we can calculate the bare soft functions order-by-order in α s .They admit α s expansions of the form S i ( z ) = δ (1 − z ) + ∞ X L =1 (cid:18) α s S ǫ π (cid:19) L (cid:18) µ M i (cid:19) ǫL − z ) ǫL S iL ( ǫ ) . (15)In Eq. (15), S ǫ is the usual MS factor, S ǫ = (cid:0) πe − γ E (cid:1) ǫ , (16)where γ E = 0 . . . . is the Euler constant and µ is the ’t Hooft scale introduced fromcontinuing the space-time dimension to D = 4 − ǫ dimension. After renormalization inthe MS scheme, it will be replaced by the renormalization scale, µ R = M i . The factor1 / (1 − z ) ǫL can be expanded in ǫ in terms of the plus distributions introduced above(Eqs. (6) and (7)), 1(1 − z ) ǫL = − δ (1 − z )2 ǫL + ∞ X m =0 ( − ǫL ) m (cid:20) ln m (1 − z )1 − z (cid:21) + . (17)The expansion coefficients of the bare soft functions, S iL ( ǫ ), are functions of ǫ and N c only. In the case of DY, the one- and two-loop results are well-known and were calculated toall orders in ǫ in Ref. [48]. The structure of the eikonal approximation is such that, throughtwo-loop order, the color degrees of freedom completely factorize from the part of the bareexpansion coefficients which encodes the non-trivial soft dynamics; by simply replacing C F everywhere with C A , the results reported in [48] can be carried over in a straightforwardmanner to the case of H as well . At N LO, many of the contributions can be converted fromthe DY language to the H language in a similar way but, overall, the situation at three loopsis significantly more complicated than at one and two loops. Although it is reasonable toexpect some correspondence by virtue of the fact that the cut eikonal diagrams contributingto H and DY differ only in the Lie algebra representation of their gluon-soft Wilson lineinteractions, there are sufficiently many color structures for both H and DY at three-looporder that finding such a correspondence requires non-trivial analysis. In fact, with dedicatedeffort, it is possible to use the known results for the single-emission, two-loop real-virtual Here, C A and C F are the usual quadratic Casimir invariants, expressed in terms of N c as C A = N c and C F = N c − N c . (18) For example, both the single-emission, two-loop real-virtual corrections computed in Refs. [31–36] andthe double-emission, one-loop real-virtual corrections discussed in the present paper fall into this category. LO result for DY can bededuced from results available in the literature and the main result of this article. Thisline of analysis was quite useful and facilitated additional non-trivial cross-checks on ourcalculation which we summarize towards the end of this paper. In order to keep the lengthof this article manageable, we defer a detailed discussion of this part of our analysis to aseparate publication which focuses on the triple-emission real contributions.As mentioned briefly in the introduction and in the last paragraph, there are several typesof real radiative corrections in the threshold limit at three-loop order involving the emissionof one or more soft partons. To obtain complete results for the S i ( ǫ ), one must computesoft single-emission, two-loop [31, 32] and one-loop squared real-virtual corrections [33–36], soft triple-emission real corrections [37], and soft double-emission, one-loop real-virtualcorrections, which are not yet publicly available . In what follows, we compute the softdouble-emission contributions to the S i ( ǫ ), referred to hereafter as S i : 2R − V3 ( ǫ ), and fill thislast remaining gap in the literature. (a) ( b )(c) ( d ) FIG. 1: Panels ( a ) and ( b ) show cut eikonal Feynman diagrams with final state gluons andpanels ( c ) and ( d ) show cut eikonal Feynman diagrams with final state fermions.To generate the full set of cut eikonal Feynman diagrams, we use QGRAF [51]. Repre-sentative cut diagrams are shown in Fig. 1. In order to test the correctness of our squared In a very recent paper [28], these contributions were computed, thereby completing the calculation of theN LO Higgs boson production cross section at threshold. R ξ gaugefor our internal gluons and lightcone gauge for our final state gluons. As we explain below,the fact that our problem supplies two natural light-like reference vectors, n and ¯ n , allowsfor a very stringent consistency check on our construction of the integrand. Once the di-agrams are generated, they are processed by an in-house FORM [52] code which dresses thecut eikonal diagrams with Feynman rules and performs all relevant Dirac and color algebra.After that, we use an in-house
Maple code to express all cut eikonal diagrams as a linearcombination of Feynman integrals belonging to one of 28 integral families introduced for theeventual purpose of performing an integration by parts reduction [53, 54] on the integrandusing Laporta’s algorithm [55]. Mapping Feynman integrals to specific integral families isnot entirely straightforward in this case and requires the application of numerous partialfraction identities to the raw integrand produced by our
FORM code. This complication isdue in part to the constraint on the kinematics imposed by the delta function in the operatordefinition of the soft functions (see Eq. (10)). We perform the integration by parts reductionusing the development version of
Reduze 2 [56] which, among other new features, supportsthe reduction of phase space integrals. The idea behind this was worked out some time agoand is usually called the reverse unitarity method [10, 57, 58]. The key insight is that, forthe purpose of integral reduction, one can use the relation δ (cid:0) k (cid:1) = − πi (cid:18) k + i − k − i (cid:19) (19)to replace delta function constraints with propagator denominators. After integration byparts reduction, the result is a rather simple linear combination of nine master integrals.To present our result in a compact fashion, let us first introduce some additional notation.Besides the Casimir invariants C A and C F , already defined above, and the number of masslessflavors, n f = 5, our result for S i : 2R − V3 ( ǫ ) will depend on an additional Casimir invariant whichitself depends on the index i : C H = C A and C DY = C F . (20)Finally, if we make the definitions Z [ k k q ] = − iπ ǫ − e γ E ǫ Z d D k Z d D k Z d D q δ + (cid:0) k (cid:1) δ + (cid:0) k (cid:1) δ (cid:16) − ( n + ¯ n ) · ( k + k ) (cid:17) D = q + i + D = ( k + k − q ) + i + D = 2 q · n + i + D = 2 k · ¯ n + i + D = 2( k − q ) · ¯ n + i + D = ( k − q ) + i + D = ( k + k ) + i + D = 2 k · n + i + D = 2 k · ¯ n + i + D = 2( k + k − q ) · ¯ n + i + , (21)the nine master integrals mentioned above are given by I ( ǫ ) = − D − D − (3 D − (cid:26)Z [ k k q ] 1 D D (cid:27) = e γ E ǫ Γ (1 − ǫ )Γ (1 + ǫ )Γ(1 − ǫ )Γ (1 − ǫ )Γ(1 + 2 ǫ )Γ(1 − ǫ ) 7 1 − ζ ǫ − ζ ǫ + 417 ζ ǫ + (cid:18) ζ ζ − ζ (cid:19) ǫ + (cid:18) − ζ
128 + 4225 ζ (cid:19) ǫ + O (cid:0) ǫ (cid:1) (22) I ( ǫ ) = 12( D − ( D − D − D − (cid:26)Z [ k k q ] 1 D D D (cid:27) = e γ E ǫ Γ (1 + ǫ )Γ (1 − ǫ )Γ(1 − ǫ )Γ(1 + 2 ǫ )Γ(1 − ǫ )= 1 − ζ ǫ − ζ ǫ − ζ ǫ + (cid:18) ζ ζ − ζ (cid:19) ǫ + (cid:18) − ζ
128 + 4489 ζ (cid:19) ǫ + O (cid:0) ǫ (cid:1) (23) I ( ǫ ) = 2( D − (3 D − (cid:26)Z [ k k q ] 1 D D D D (cid:27) = e γ E ǫ Γ (1 − ǫ )Γ (1 + ǫ )Γ(1 + 2 ǫ )Γ(1 − ǫ )= 1 − ζ ǫ − ζ ǫ + 225 ζ ǫ + (cid:18) ζ ζ − ζ (cid:19) ǫ + (cid:18) − ζ
128 + 4761 ζ (cid:19) ǫ + O (cid:0) ǫ (cid:1) (24) I ( ǫ ) = − D − ( D − D − (cid:26)Z [ k k q ] 1 D D D D (cid:27) = e γ E ǫ Γ (1 − ǫ )Γ (1 + ǫ )Γ(1 − ǫ ) F (1 , − ǫ, ǫ ; 1 − ǫ, − ǫ ; 1)Γ (1 − ǫ )Γ(1 + 2 ǫ )Γ(1 − ǫ )= 1 − ζ ǫ − ζ ǫ + 713 ζ ǫ + (cid:18) ζ ζ − ζ (cid:19) ǫ + (cid:18) − ζ
128 + 5005 ζ (cid:19) ǫ + O (cid:0) ǫ (cid:1) (25) I ( ǫ ) = ( D − ( D − D − (cid:26)Z [ k k q ] 1 D D D D (cid:27) = e γ E ǫ ǫ Γ (1 − ǫ )Γ (1 + ǫ )Γ(1 − ǫ ) F (1 , ǫ, ǫ ; 2 − ǫ, ǫ ; 1)Γ (1 − ǫ )Γ(1 + 2 ǫ )Γ(1 − ǫ )Γ(2 − ǫ )Γ(2 + ǫ ) (26)= ζ ǫ + 3 ζ ǫ − ζ ǫ + (cid:18) − ζ ζ ζ (cid:19) ǫ + (cid:18) − ζ − ζ (cid:19) ǫ + O (cid:0) ǫ (cid:1) I ( ǫ ) = 2( D − (2 D − D − (cid:26)Z [ k k q ] 1 D D D D (cid:27) = e γ E ǫ Γ (1 + ǫ )Γ (1 − ǫ )Γ (1 − ǫ ) F ( − ǫ, ǫ, ǫ ; 1 − ǫ, − ǫ ; 1)Γ(1 + 2 ǫ )Γ(1 − ǫ )Γ(1 − ǫ )= 1 − ζ ǫ − ζ ǫ − ζ ǫ + (cid:18) ζ ζ − ζ (cid:19) ǫ (cid:18) − ζ
128 + 5483 ζ (cid:19) ǫ + O (cid:0) ǫ (cid:1) (27) I ( ǫ ) = 3( D −
10 Re (cid:26)Z [ k k q ] 1 D D D D D (cid:27) = 2 e γ E ǫ Γ (1 − ǫ )Γ (1 + ǫ )Γ(1 − ǫ ) F (1 , , − ǫ ; 1 − ǫ, − ǫ ; 1)5Γ (1 − ǫ )Γ(1 + 2 ǫ )Γ(1 − ǫ )= 1 − ζ ǫ − ζ ǫ + 5757 ζ ǫ + (cid:18) ζ ζ − ζ (cid:19) ǫ + (cid:18) − ζ
640 + 38261 ζ (cid:19) ǫ + O (cid:0) ǫ (cid:1) (28) I ( ǫ ) = 3( D − (3 D − D −
13) Re (cid:26)Z [ k k q ] 1 D D D D D D (cid:27) = 9 e γ E ǫ Γ (1 + ǫ )Γ (1 − ǫ )Γ(1 − ǫ ) F (1 , − ǫ, − ǫ, ǫ ; 1 − ǫ, − ǫ, − ǫ ; 1)11Γ(1 + 2 ǫ )Γ(1 − ǫ )= 1 − ζ ǫ − ζ ǫ − ζ ǫ + (cid:18) ζ ζ − ζ (cid:19) ǫ + (cid:18) − ζ ζ (cid:19) ǫ + O (cid:0) ǫ (cid:1) (29) I ( ǫ ) = 9( D − (3 D − D −
13) Re (cid:26)Z [ k k q ] 1 D D D D D D D (cid:27) + 13 D − (cid:16) − I ( ǫ ) + 3 I ( ǫ ) − I ( ǫ ) (cid:17) = e γ E ǫ Γ (1 − ǫ )Γ (1 + ǫ )Γ (1 − ǫ )2Γ( − ǫ )Γ(1 + 2 ǫ )Γ(1 − ǫ )Γ(1 − ǫ )Γ(2 + ǫ ) ×× F (cid:18) — , , − ǫ − ǫ, − ǫ, − ǫ − ǫ ǫ − ǫ ; 1 , (cid:19) − ǫ ) (cid:16) − I ( ǫ ) + 3 I ( ǫ ) − I ( ǫ ) (cid:17) = 1 − ζ ǫ − ζ ǫ − ζ ǫ + (cid:18) − ζ ζ − ζ (cid:19) ǫ + (cid:18) − ζ − ζ (cid:19) ǫ + O (cid:0) ǫ (cid:1) . (30)While integrals I ( ǫ ) through I ( ǫ ) can be simply expressed in terms of gamma functions andgeneralized hypergeometric functions, I ( ǫ ) is more complicated. In particular, it involvesthe Kamp´e de F´eriet function (see Ref. [59], Eq. (1.3.2.1) for our Kamp´e de F´eriet functionconventions) F (cid:18) — , , − ǫ − ǫ, − ǫ, − ǫ − ǫ ǫ − ǫ ; 1 , (cid:19) = ∞ X m =0 ∞ X n =0 m !( − ǫ ) m [( − ǫ ) n ] ( − ǫ ) n n !(1 − ǫ ) n (2 + ǫ ) m ( − ǫ ) n + m , where ( r ) s = Γ( r + s )Γ( r ) . (31)9 few words about our derivation of the master integrals are in order. Although, forthe most part, we wish to defer the discussion of the novel method used to do the integralsto a more technical future publication, let us briefly describe some consistency checks onEqs. (22)-(30) that we found useful. For all integrals, a good first step is to integrate out thevirtual momentum q since, in all cases, one obtains simple expressions built out of gammafunctions and Gauss hypergeometric functions. The numerical sector decomposition [60, 61]code FIESTA 3 [62] was successfully used to check our virtual integrations for sign errors. Unfortunately, more stringent checks on the virtual integrations with
FIESTA 3 could notbe performed with confidence due to apparent stability issues with the numerical integra-tion routines. For all of the integrals except I ( ǫ ) and I ( ǫ ), we found that the phase spaceintegrals over k and k that remain once q has been integrated out could be numericallychecked in Mathematica for appropriate real, negative values of ǫ using a parametrizationemployed in earlier work (see Refs. [63–65]). In all cases where our parametrization ap-plies , we found that our analytical results could be checked numerically at the level of thephase space integrals to at least seven significant digits for several appropriately chosen butessentially random values of ǫ . Fortunately, the two integrals whose phase space integralsdo not admit a straightforward parametrization for technical reasons, I ( ǫ ) and I ( ǫ ), areactually the simplest to set up in our non-standard analytical approach. Our results for I ( ǫ ) and I ( ǫ ) were verified a posteriori by using them to carry out the stringent globalconsistency checks described at the end of this work. To derive the Taylor expansions ofintegrals I ( ǫ ) through I ( ǫ ), we made extensive use of the Mathematica package
HypExp [66, 67]. As a sanity check, we performed a completely independent expansion of the masterintegral I ( ǫ ) using general-purpose integration scripts written by one of us (see Ref. [65])and found complete agreement.The master integral I ( ǫ ) is somewhat more subtle than I ( ǫ )- I ( ǫ ). In an attempt toexpand the second line of Eq. (30) in ǫ , we encountered the single-parameter integral Z d u u − − ǫ (1 − u ) − − ǫ (cid:16) F ( − ǫ, − ǫ ; 1 − ǫ ; 1 − u ) (cid:17) . (32)Historically, such integrals were considered pathological due to the fact that, na¨ıvely, they donot appear to be sector-decomposable. The problem is that the integrand of (32) has small u asymptotics of the too-singular form u − − ǫ . Integrals such as (32) have appeared in otherQCD computations. For example, in Ref. [68], a dedicated unitarity-based sewing methodwas used to avoid classes of integrals with similar power-law singularities. In fact, we foundthat we were able to treat integrals like (32) in a very direct way. First, observe that theabove integral converges for all complex ǫ such that Re( ǫ ) < − /
3. In spite of the fact thatthe integral does not converge in a neighborhood of ǫ = 0, one can expand (32) in a Laurentseries about ǫ = 0 anyway by carefully applying the principle of analytical continuation. Tounderstand this, begin by considering a fictitious, alternative representation of our result,Eq. (33), written in terms of integrals free of power-law divergences which converge in aneighborhood of ǫ = 0. This hypothetical representation should always exist if all of thedivergences in the calculation are regulated by ǫ . Using the fact that all Feynman integralsin dimensional regularization are analytic functions of ǫ , one can perform an analytical It is worth pointing out that, although we have absorbed the operation of taking the real part into thedefinitions of our master integrals, they are complex functions of ǫ at the level of FIESTA and it is thereforeessential to use version 3 since earlier versions of the code support only Euclidean kinematics. Roughly speaking, we expect our methods to apply in a straightforward manner if the dependence of thephase space integrals on the dot product k · k is of the form ( k · k ) p ( ǫ ) for some linear polynomial, p ( ǫ ). Of course, values of ǫ which lead to a convergent numerical integration must be chosen. ǫ which happens to simultaneously lie in the region of convergenceof (32) and I ( ǫ ) through I ( ǫ ). Now, by performing an integration by parts reduction onthe fictitious representation, one can recover the precise form of Eq. (33), written in termsof our preferred basis of master integrals. Due to the fact that, by assumption, the point inthe complex ǫ plane to which the original analytical continuation was performed lies withinthe domain of convergence of all of our master integrals, one can derive the closed formulasfor I ( ǫ )- I ( ǫ ) given above (the second lines of Eqs. (22)-(30)) in this region of the complex ǫ plane. Finally, one can analytically continue back to the point ǫ = 0 and derive Laurentexpansions for all of the masters integrals using standard techniques. Carrying out thisprogram explicitly was necessary to Laurent expand I ( ǫ ) about ǫ = 0 and our analysis willbe presented in a future publication which features an in-depth discussion of the masterintegrals.Our choice of integral basis also requires some explanation. As is clear from the Taylorseries expansions of our master integrals, the definitions we have made are such that all ofthe I j ( ǫ ) are pure functions in the sense of reference [69]. As we shall see, the all-orders-in- ǫ result for S i : 2R − V3 ( ǫ ) assumes a particularly simple form when expressed in terms of theintegral basis defined above. After integration by parts reduction, it can be written as S i : 2R − V3 ( ǫ ) = C i C A (cid:20)(cid:18) − D −
4) + 63689( D − − D −
1) + 1427227( D − + 180( D − − D − − D − + 14089( D − + 25627( D − − D − (cid:19) I ( ǫ )+ (cid:18) − D −
4) + 20483( D − − D − D −
1) + 5286481( D − − D − + 56329( D − − D − − D − (cid:19) I ( ǫ ) + 768( D − I ( ǫ )+ (cid:18) D − − D −
3) + 6445( D − − D − − D − + 3215( D − (cid:19) I ( ǫ ) + (cid:18) − D −
4) + 1283( D − − D −
1) + 1283( D − + 5123( D − − D − (cid:19) I ( ǫ ) + (cid:18) D − − D −
3) + 8961215( D − − D − + 857627( D − − D − − D − + 6415( D − (cid:19) I ( ǫ ) − D − I ( ǫ ) − D − I ( ǫ ) − D − I ( ǫ ) (cid:21) − C i C A D − I ( ǫ )+ C i C A n f (cid:20)(cid:18) − D − − D −
3) + 16 D − − D − − D − + 150481( D − − D − − D − + 2569( D − + 169( D − (cid:19) I ( ǫ )+ (cid:18) D − − D −
3) + 1123( D −
2) + 240641215( D − − D −
11 64( D − + 256027( D − − D − + 1615( D − (cid:19) I ( ǫ ) + (cid:18) − D − D − − D − D −
1) + 643( D − − D − (cid:19) I ( ǫ )+ (cid:18) D − − D −
3) + 32 D − − D − − D − + 3215( D − (cid:19) I ( ǫ ) + (cid:18) − D −
4) + 2563( D − − D − − D − D − − D − + 5129( D − + 17921215( D − (cid:19) I ( ǫ ) (cid:21) + C i C F n f (cid:18) − D −
3) + 160243( D − − D − + 43904243( D − − D − + 243227( D − − D − (cid:19) I ( ǫ )+ C i n f (cid:18) − D − − D − − D − + 643( D − − D − + 25627( D − (cid:19) I ( ǫ ) . (33)The compact form of Eq. (33) is due in large part to the absence of spurious poles in thecoefficients of our master integrals; the coefficients of the basis integrals in Eq. (33) havepoles at positive integer values of D only. By rewriting the above expression in terms of amore conventional basis such as the one supplied by Reduze 2 out of the box, it becomesclear that this intriguing feature of our result is a direct consequence of the fact that the I j ( ǫ ) are pure functions. For a problem with similar kinematics, the computation of thethree-loop form factors in massless QCD, it was observed in Ref. [70] (in slightly differentlanguage) that an integral basis of pure functions could be written down for all masterintegrals. However, to the best of our knowledge, our paper is the first to show that,even for single-scale problems, it is reasonable to expect the pole structure of the reducedintegrand to dramatically simplify when it is expressed in terms of a basis of pure functions.Substituting the Taylor expansions of the master integrals, Eqs. (22)-(30), into Eq. (33) andexpanding in ǫ , we find S i : 2R − V3 ( ǫ ) = C i C A (cid:20) ǫ + 883 ǫ + (cid:18) − ζ (cid:19) ǫ + (cid:18) − ζ − ζ (cid:19) ǫ + (cid:18) − ζ − ζ − ζ (cid:19) ǫ + 870688729 − ζ − ζ − ζ ζ ζ − ζ (cid:18) − ζ − ζ − ζ ζ ζ − ζ − ζ
144 + 125284 ζ (cid:19) ǫ (cid:21) + C i C A (cid:20) ǫ − ζ ǫ − ζ ǫ + 675 ζ ǫ + 54648 ζ ζ ζ (cid:18) − ζ ζ (cid:19) ǫ (cid:21) + C i C A n f (cid:20) − ǫ − ǫ + (cid:18) ζ (cid:19) ǫ + (cid:18) − − ζ
27+ 880 ζ (cid:19) ǫ − − ζ − ζ
27 + 469 ζ (cid:18) − − ζ − ζ
81 + 5941 ζ − ζ ζ + 11712 ζ (cid:19) ǫ (cid:21) + C i C F n f (cid:20) − ǫ − ǫ + (cid:18) − ζ (cid:19) ǫ + (cid:18) − ζ
9+ 2080 ζ (cid:19) ǫ − ζ
27 + 19760 ζ − ζ (cid:18) − ζ
81 + 179920 ζ − ζ − ζ ζ + 80672 ζ (cid:19) ǫ (cid:21) + C i n f (cid:20) − ǫ − ǫ + (cid:18) −
323 + 176 ζ (cid:19) ǫ − ζ
27 + 2080 ζ (cid:18) − ζ + 20800 ζ − ζ (cid:19) ǫ (cid:21) + O (cid:0) ǫ (cid:1) . (34)The Laurent series expansion of S i : 2R − V3 ( ǫ ) given above is the main result of this paper.In order to ensure the correctness of our results, we found it useful to perform severalconsistency checks at various stages of the calculation. First, recall that, while constructingthe raw integrand, we made use of a physical lightcone gauge for our final state gluons. Itturns out that running our code twice, once with n as the lightcone reference vector for allgluons and once with ¯ n as the lightcone reference vector for all gluons, allows for a verystringent check at the stage of integral reduction. The reason for this is as follows. Atthe level of the unreduced integrand, using a physical gauge with some choice of referencevector breaks the n − ¯ n symmetry of the expression expected by virtue of the fact thatour soft Wilson line operators are back-to-back. This symmetry, however, is restored afterintegral reduction and, most importantly, before reducing all Feynman integrals to masters,we observed that many integrals enter one or the other of the two expressions we derived butnot both. Thus, reducing our two seemingly different expressions for the raw integrand downto Eq. (33) furnished a non-trivial check on the calculation. This check was complementaryto a more conventional check on the gauge invariance of Eq. (33); we employed a general R ξ gauge for our internal gluons and confirmed that the gauge parameter, ξ , drops out ofour final results.The final and most important check on our calculation was an explicit comparison to therelevant literature. After extracting the one-, two-, and three-loop expansion coefficients ofthe hard functions from Refs. [23–26, 71, 72] and deducing the one-, two-, and three-loopexpansion coefficients of the renormalized soft functions from our calculation and variouscomplete and partial results existing in the literature [31–37, 48], we were able to use Eq.(8) to reproduce the full tower of plus distributions of the form (cid:20) ln m (1 − z )(1 − z ) (cid:21) + , m = 0 . . . δ (1 − z ) terms in the full N LOresult obtained recently by the authors of Ref. [28] and, for the case of Drell-Yan leptonproduction, we were able to reproduce the prediction for the δ (1 − z ) terms made recentlyin Ref. [76] using completely different methods. Of course, this is really a check on the δ (1 − z ) term coming from the double-emission, one-loop real-virtual corrections treated inthis paper, not the δ (1 − z ) term in the full N LO result because, so far, no independentcalculation of the triple-emission real contributions has appeared confirming that the anal-ysis of [37] is correct. Nevertheless, the fact that we were able to completely reproduce theresults presented in Ref. [28] is highly non-trivial and very encouraging.We have taken a decisive step in this work towards a fully independent calculation of thetotal cross section at threshold for gluon fusion Higgs production and Drell-Yan lepton pairproduction at N LO. In this paper, we focused on the calculation of the double-emission,one-loop real-virtual contributions to the bare three-loop threshold soft functions for gluonfusion Higgs production and Drell-Yan lepton production. As explained above, we foundresults completely consistent with all of the available literature on the subject. During thecourse of our calculation, we made significant technical advances in the evaluation of cuteikonal phase space integrals which allowed us to compute all of our master integrals to allorders in ǫ . It would be very interesting to try and carry out a similar program for themaster integrals that one finds in the computation of the triple-emission real contributionsto the bare three-loop soft functions. This would be useful because complete all-orders-in- ǫ results for the three-loop bare soft functions would save future researchers the trouble ofrevisiting the calculations performed in Ref. [37] if it turns out to be necessary to calculateto one order higher (to N LO) in perturbative QCD. In fact, it may well be that our masterintegrals find useful application as input integrals to the systems of differential equationsthat will presumably be derived in the future when an attempt is made to calculate the full z -dependent partonic cross sections for the processes we have considered at threshold.We also noted in this work that the use of an integral basis composed of pure functionsleads to remarkable simplifications at the level of the reduced integrand, Eq. (33). In starkcontrast to more conventional bases of master integrals, our basis of pure functions is suchthat the coefficients of the master integrals have poles only at D = 1, 2, 3, 4, and 6. Weargued that this property is a direct consequence of the niceness of our integral basis whichwould certainly not hold in general. In fact, the most complicated pure function in our basiswas identified by demanding that spurious singularities at D = 14 / ǫ expressions remains to be seen and should help us to more accurately assess the scope ofapplicability of the methods that we developed to carry out the calculations discussed inthis paper. Our hope is that the N LO threshold production cross section calculations con-sidered in this work and elsewhere [28, 31–37] will, once N LO parton distribution functionfits become available, allow for more accurate inclusive predictions than ever before andsignificantly advance the physics program of the Large Hadron Collider.14
CKNOWLEDGMENTS
We are grateful to Sven-Olaf Moch for useful correspondence. The Feynman diagramsare drawn with JaxoDraw [77], based on AxoDraw [78]. The research of YL and HXZ is sup-ported by the US Department of Energy under contract DEAC0276SF00515. The researchof AvM is supported in part by the Research Center
Elementary Forces and MathematicalFoundations (EMG) of the Johannes Gutenberg University of Mainz and by the German Re-search Foundation (DFG). The research of RMS is supported by the ERC Advanced GrantEFT4LHC of the European Research Council, the Cluster of Excellence Precision Physics,Fundamental Interactions and Structure of Matter (PRISMA-EXC 1098). [1] G. Altarelli, R. K. Ellis and G. Martinelli, “Large Perturbative Corrections to the Drell-YanProcess in QCD,” Nucl. Phys. B , 461 (1979).[2] S. Dawson, “Radiative corrections to Higgs boson production,” Nucl. Phys. B , 283 (1991).[3] R. Hamberg, W. L. van Neerven and T. Matsuura, “A Complete calculation of the order α s correction to the Drell-Yan K factor,” Nucl. Phys. B , 343 (1991) [Erratum-ibid. B ,403 (2002)].[4] D. Graudenz, M. Spira and P. M. Zerwas, “QCD corrections to Higgs boson production atproton proton colliders,” Phys. Rev. Lett. , 1372 (1993).[5] M. Spira, A. Djouadi, D. Graudenz and P. M. Zerwas, “Higgs boson production at the LHC,”Nucl. Phys. B , 17 (1995) [hep-ph/9504378].[6] A. Djouadi, M. Spira and P. M. Zerwas, “QCD corrections to hadronic Higgs decays,” Z.Phys. C , 427 (1996) [hep-ph/9511344].[7] M. Spira, “QCD effects in Higgs physics,” Fortsch. Phys. , 203 (1998) [hep-ph/9705337].[8] S. Catani, D. de Florian and M. Grazzini, “Higgs production in hadron collisions: Soft andvirtual QCD corrections at NNLO,” JHEP , 025 (2001) [hep-ph/0102227].[9] R. V. Harlander and W. B. Kilgore, “Next-to-next-to-leading order Higgs production at hadroncolliders,” Phys. Rev. Lett. , 201801 (2002) [hep-ph/0201206].[10] C. Anastasiou and K. Melnikov, “Higgs boson production at hadron colliders in NNLO QCD,”Nucl. Phys. B , 220 (2002) [hep-ph/0207004].[11] V. Ravindran, J. Smith and W. L. van Neerven, “NNLO corrections to the total cross-sectionfor Higgs boson production in hadron hadron collisions,” Nucl. Phys. B , 325 (2003)[hep-ph/0302135].[12] R. Harlander and P. Kant, “Higgs production and decay: Analytic results at next-to-leadingorder QCD,” JHEP , 015 (2005) [hep-ph/0509189].[13] K. G. Chetyrkin, B. A. Kniehl and M. Steinhauser, “Decoupling relations to O (cid:0) α s (cid:1) and theirconnection to low-energy theorems,” Nucl. Phys. B , 61 (1998) [hep-ph/9708255].[14] Y. Schr¨oder and M. Steinhauser, “Four-loop decoupling relations for the strong coupling,”JHEP , 051 (2006) [hep-ph/0512058].[15] K. G. Chetyrkin, J. H. K¨uhn and C. Sturm, “QCD decoupling at four loops,” Nucl. Phys. B , 121 (2006) [hep-ph/0512060].[16] O. V. Tarasov, A. A. Vladimirov and A. Y. Zharkov, “The Gell-Mann-Low Function of QCDin the Three Loop Approximation,” Phys. Lett. B , 429 (1980).[17] S. A. Larin and J. A. M. Vermaseren, “The Three loop QCD Beta function and anomalous imensions,” Phys. Lett. B , 334 (1993) [hep-ph/9302208].[18] T. van Ritbergen, J. A. M. Vermaseren and S. A. Larin, “The Four loop beta function inquantum chromodynamics,” Phys. Lett. B , 379 (1997) [hep-ph/9701390].[19] M. Czakon, “The Four-loop QCD beta-function and anomalous dimensions,” Nucl. Phys. B , 485 (2005) [hep-ph/0411261].[20] C. Anastasiou, S. Buehler, C. Duhr and F. Herzog, “NNLO phase space master integrals fortwo-to-one inclusive cross sections in dimensional regularization,” JHEP , 062 (2012)[arXiv:1208.3130 [hep-ph]].[21] M. H¨oschele, J. Hoff, A. Pak, M. Steinhauser and T. Ueda, “Higgs boson production at theLHC: NNLO partonic cross sections through order ǫ and convolutions with splitting functionsto N LO,” Phys. Lett. B , 244 (2013) [arXiv:1211.6559 [hep-ph]].[22] S. Buehler and A. Lazopoulos, “Scale dependence and collinear subtraction terms for Higgsproduction in gluon fusion at N LO,” JHEP , 096 (2013) [arXiv:1306.2223 [hep-ph]].[23] P. A. Baikov, K. G. Chetyrkin, A. V. Smirnov, V. A. Smirnov and M. Steinhauser, “Quarkand gluon form factors to three loops,” Phys. Rev. Lett. , 212002 (2009) [arXiv:0902.3519[hep-ph]].[24] R. N. Lee, A. V. Smirnov and V. A. Smirnov, “Analytic Results for Massless Three-LoopForm Factors,” JHEP , 020 (2010) [arXiv:1001.2887 [hep-ph]].[25] T. Gehrmann, E. W. N. Glover, T. Huber, N. Ikizlerli and C. Studerus, “Calculation of thequark and gluon form factors to three loops in QCD,” JHEP , 094 (2010) [arXiv:1004.3653[hep-ph]].[26] T. Gehrmann, E. W. N. Glover, T. Huber, N. Ikizlerli and C. Studerus, “The quark andgluon form factors to three loops in QCD through to O (cid:0) ǫ (cid:1) ,” JHEP , 102 (2010)[arXiv:1010.4478 [hep-ph]].[27] R. D. Ball, M. Bonvini, S. Forte, S. Marzani and G. Ridolfi, “Higgs production in gluon fusionbeyond NNLO,” Nucl. Phys. B , 746 (2013) [arXiv:1303.3590 [hep-ph]].[28] C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog and B. Mistlberger,“Higgs boson gluon-fusion production at threshold in N LO QCD,” arXiv:1403.4616 [hep-ph].[29] T. Gehrmann, M. Jaquier, E. W. N. Glover and A. Koukoutsakis, “Two-Loop QCD Correc-tions to the Helicity Amplitudes for H → , 056 (2012) [arXiv:1112.3554[hep-ph]].[30] S. D. Badger and E. W. N. Glover, “Two loop splitting functions in QCD,” JHEP , 040(2004) [hep-ph/0405236].[31] Y. Li and H. X. Zhu, “Single soft gluon emission at two loops,” JHEP , 080 (2013)[arXiv:1309.4391 [hep-ph]].[32] C. Duhr and T. Gehrmann, “The two-loop soft current in dimensional regularization,” Phys.Lett. B , 452 (2013) [arXiv:1309.4393 [hep-ph]].[33] C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger, “Real-virtual contributionsto the inclusive Higgs cross-section at N LO,” JHEP , 088 (2013) [arXiv:1311.1425 [hep-ph]].[34] Z. Bern, V. Del Duca, W. B. Kilgore and C. R. Schmidt, “The infrared behavior of oneloop QCD amplitudes at next-to-next-to leading order,” Phys. Rev. D , 116001 (1999)[hep-ph/9903516].[35] S. Catani and M. Grazzini, “The soft gluon current at one loop order,” Nucl. Phys. B ,435 (2000) [hep-ph/0007142].[36] W. B. Kilgore, Phys. Rev. D , 073008 (2014) [arXiv:1312.1296 [hep-ph]].
37] C. Anastasiou, C. Duhr, F. Dulat and B. Mistlberger, “Soft triple-real radiation for Higgsproduction at N LO,” JHEP , 003 (2013) [arXiv:1302.4379 [hep-ph]].[38] S. Moch and A. Vogt, “Higher-order soft corrections to lepton pair and Higgs boson produc-tion,” Phys. Lett. B , 48 (2005) [hep-ph/0508265].[39] M. A. Shifman, A. I. Vainshtein, M. B. Voloshin and V. I. Zakharov, “Low-Energy Theoremsfor Higgs Boson Couplings to Photons,” Sov. J. Nucl. Phys. , 711 (1979) [Yad. Fiz. ,1368 (1979)].[40] C. W. Bauer, S. Fleming and M. E. Luke, “Summing Sudakov logarithms in B → X s + γ ineffective field theory,” Phys. Rev. D , 014006 (2000) [hep-ph/0005275].[41] C. W. Bauer, S. Fleming, D. Pirjol and I. W. Stewart, “An Effective field theory for collinearand soft gluons: Heavy to light decays,” Phys. Rev. D , 114020 (2001) [hep-ph/0011336].[42] C. W. Bauer, D. Pirjol and I. W. Stewart, “Soft collinear factorization in effective field theory,”Phys. Rev. D , 054022 (2002) [hep-ph/0109045].[43] M. Beneke, A. P. Chapovsky, M. Diehl and T. Feldmann, “Soft collinear effective the-ory and heavy to light currents beyond leading power,” Nucl. Phys. B , 431 (2002)[hep-ph/0206152].[44] A. Idilbi, X. Ji and F. Yuan, “Resummation of threshold logarithms in effective field theoryfor DIS, Drell-Yan and Higgs production,” Nucl. Phys. B , 42 (2006) [hep-ph/0605068].[45] T. Becher, M. Neubert and G. Xu, “Dynamical Threshold Enhancement and Resummationin Drell-Yan Production,” JHEP , 030 (2008) [arXiv:0710.0680 [hep-ph]].[46] V. Ahrens, T. Becher, M. Neubert and L. L. Yang, “Origin of the Large PerturbativeCorrections to Higgs Production at Hadron Colliders,” Phys. Rev. D , 033013 (2009)[arXiv:0808.3008 [hep-ph]].[47] J. Chay, C. Kim, Y. G. Kim and J. -P. Lee, “Soft Wilson lines in soft-collinear effectivetheory,” Phys. Rev. D , 056001 (2005) [hep-ph/0412110].[48] A. V. Belitsky, “Two loop renormalization of Wilson loop for Drell-Yan production,” Phys.Lett. B , 307 (1998) [hep-ph/9808389].[49] J. G. M. Gatheral, “Exponentiation of Eikonal Cross-sections in Nonabelian Gauge Theories,”Phys. Lett. B , 90 (1983).[50] J. Frenkel and J. C. Taylor, “Nonabelian Eikonal Exponentiation,” Nucl. Phys. B , 231(1984).[51] P. Nogueira, “Automatic Feynman graph generation,” J. Comput. Phys. , 279 (1993).[52] J. A. M. Vermaseren, “New features of FORM,” math-ph/0010025.[53] F. V. Tkachov, “A Theorem on Analytical Calculability of Four Loop Renormalization GroupFunctions,” Phys. Lett. B , 65 (1981).[54] K. G. Chetyrkin and F. V. Tkachov, “Integration by Parts: The Algorithm to Calculate betaFunctions in 4 Loops,” Nucl. Phys. B , 159 (1981).[55] S. Laporta, “High precision calculation of multiloop Feynman integrals by difference equa-tions,” Int. J. Mod. Phys. A , 5087 (2000) [hep-ph/0102033].[56] A. von Manteuffel and C. Studerus, “Reduze 2 - Distributed Feynman Integral Reduction,”arXiv:1201.4330 [hep-ph].[57] C. Anastasiou, L. J. Dixon, K. Melnikov and F. Petriello, “Dilepton rapidity distribution in theDrell-Yan process at NNLO in QCD,” Phys. Rev. Lett. , 182002 (2003) [hep-ph/0306192].[58] C. Anastasiou, L. J. Dixon, K. Melnikov and F. Petriello, “High precision QCD at hadroncolliders: Electroweak gauge boson rapidity distributions at NNLO,” Phys. Rev. D , 094008(2004) [hep-ph/0312266].
59] H. Exton, “Handbook of Hypergeometric Integrals,” Ellis Hornwood Limited, Chichester,(1978).[60] T. Binoth and G. Heinrich, “An automatized algorithm to compute infrared divergent multi-loop integrals,” Nucl. Phys. B , 741 (2000) [hep-ph/0004013].[61] C. Bogner and S. Weinzierl, “Resolution of singularities for multi-loop integrals,” Comput.Phys. Commun. , 596 (2008) [arXiv:0709.4092 [hep-ph]].[62] A. V. Smirnov, Comput. Phys. Commun. , 2090 (2014) [arXiv:1312.3186 [hep-ph]].[63] R. Kelley, M. D. Schwartz, R. M. Schabinger and H. X. Zhu, “The two-loop hemisphere softfunction,” Phys. Rev. D , 045022 (2011) [arXiv:1105.3676 [hep-ph]].[64] R. Kelley, M. D. Schwartz, R. M. Schabinger and H. X. Zhu, “Jet Mass with a Jet Veto atTwo Loops and the Universality of Non-Global Structure,” Phys. Rev. D , 054017 (2012)[arXiv:1112.3343 [hep-ph]].[65] A. von Manteuffel, R. M. Schabinger and H. X. Zhu, “The Complete Two-Loop Inte-grated Jet Thrust Distribution In Soft-Collinear Effective Theory,” JHEP , 139 (2014)[arXiv:1309.3560 [hep-ph]].[66] T. Huber and D. Maitre, “HypExp: A Mathematica package for expanding hypergeomet-ric functions around integer-valued parameters,” Comput. Phys. Commun. , 122 (2006)[hep-ph/0507094].[67] T. Huber and D. Maitre, “HypExp 2, Expanding Hypergeometric Functions about Half-IntegerParameters,” Comput. Phys. Commun. , 755 (2008) [arXiv:0708.2443 [hep-ph]].[68] Z. Bern, L. J. Dixon and D. A. Kosower, “Two-loop g → gg splitting amplitudes in QCD,”JHEP , 012 (2004) [hep-ph/0404293].[69] J. M. Henn, “Multiloop integrals in dimensional regularization made simple,” Phys. Rev. Lett. , no. 25, 251601 (2013) [arXiv:1304.1806 [hep-th]].[70] R. N. Lee and V. A. Smirnov, “Analytic Epsilon Expansions of Master Integrals Correspondingto Massless Three-Loop Form Factors and Three-Loop g − , 102 (2011) [arXiv:1010.1334 [hep-ph]].[71] T. Matsuura, S. C. van der Marck and W. L. van Neerven, “The Calculation of the SecondOrder Soft and Virtual Contributions to the Drell-Yan Cross-Section,” Nucl. Phys. B ,570 (1989).[72] R. V. Harlander, “Virtual corrections to gg → H to two loops in the heavy top limit,” Phys.Lett. B , 74 (2000) [hep-ph/0007289].[73] E. Laenen and L. Magnea, “Threshold resummation for electroweak annihilation from DISdata,” Phys. Lett. B , 270 (2006) [hep-ph/0508284].[74] V. Ravindran, “Higher-order threshold effects to inclusive processes in QCD,” Nucl. Phys. B , 173 (2006) [hep-ph/0603041].[75] V. Ahrens, T. Becher, M. Neubert and L. L. Yang, “Renormalization-Group ImprovedPrediction for Higgs Production at Hadron Colliders,” Eur. Phys. J. C , 333 (2009)[arXiv:0809.4283 [hep-ph]].[76] T. Ahmed, M. Mahakhud, N. Rana and V. Ravindran, “Drell-Yan production at threshold inN LO QCD,” arXiv:1404.0366 [hep-ph].[77] D. Binosi and L. Theussl, “JaxoDraw: A Graphical user interface for drawing Feynman dia-grams,” Comput. Phys. Commun. , 76 (2004) [hep-ph/0309015].[78] J. A. M. Vermaseren, “Axodraw,” Comput. Phys. Commun. , 45 (1994)., 45 (1994).