Higgs boson pair production via gluon fusion at N 3 LO in QCD
HHiggs boson pair production via gluon fusion at N LO in QCD
Long-Bin Chen, ∗ Hai Tao Li, † Hua-Sheng Shao, ‡ and Jian Wang § School of Physics and Electronic Engineering, Guangzhou University, Guangzhou 510006, China Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, 87545, USA Laboratoire de Physique Th´eorique et Hautes Energies (LPTHE), UMR 7589,Sorbonne Universit´e et CNRS, 4 place Jussieu, 75252 Paris Cedex 05, France School of Physics, Shandong University, Jinan, Shandong 250100, China
We present next-to-next-to-next-to-leading order (N LO) QCD predictions for the Higgs bosonpair production via gluon fusion at hadron colliders in the infinite top-quark mass limit. Besidesthe inclusive total cross sections at various collision energies, we also provide the invariant massdistribution of the Higgs boson pair. Our results show that the N LO QCD corrections enhancethe next-to-next-to-leading order cross section by 3 .
0% (2 . √ s = 13 (100) TeV, while thescale uncertainty is reduced substantially below 3% (2%). We also find that a judicious scale choicecan significantly improve the perturbative convergence. For the invariant mass distribution, ourcalculation demonstrates that the N LO corrections improve the scale dependence but almost donot change the shape.
INTRODUCTION
The discovery of the Higgs boson [1, 2] marks the com-pletion of the standard model (SM) of particle physicsand the start of a new era for the physics studies at theLHC. The next primary goal of the LHC is to preciselypin down its interactions with other SM particles or itself.In particular, the precision study of the Higgs potentialis ultimately crucial for understanding the electroweaksymmetry breaking mechanism. At the LHC, it has beenfound that its interaction couplings with massive gaugebosons and fermions agree with the SM expectations [3–6], while there are only quite weak constraints on the tri-linear Higgs self-coupling [7, 8]. However, in the future,the experimental probe will be significantly improved asthe increase of the integrated luminosity and the colli-sion energy [9, 10], and/or by employing novel analyzingmethods [11]. Theoretically, there indeed exist beyond-the-SM (BSM) models, in which the trilinear Higgs self-coupling deviates from the SM value by about 100% whilethe Higgs couplings to gauge bosons and fermions are al-most SM-like [12]. Therefore, the precise measurementof the trilinear Higgs self-coupling at the LHC and futurehigh-energy colliders would be of paramount importanceto test the SM and to explore the elusive BSM signals.The direct manner to probe the trilinear Higgs self-coupling at a hadron collider, such as the LHC, is via theHiggs boson pair production. Like the single Higgs case,the gluon-gluon fusion (ggF) channel is dominant, whileother channels like vector-boson fusion (VBF) are at leastone order of magnitude lower in their yields [13, 14]. Sim-ilarly to the cross section of single Higgs boson produc-tion, the ggF di-Higgs cross section is plagued with largetheoretical uncertainties, dominated by the QCD scaleuncertainty [15] and the top-quark mass scheme depen-dence [16, 17]. The computations of the cross sectionhave been carried out both in the infinite top-quark masslimit and with full top-quark mass dependence. In the infinite top-quark mass limit, the next-to-leading order (NLO) QCD correction was known twentyyears ago [15], and the NNLO QCD calculation wasperformed recently [18–21]. In addition, the effect ofsoft gluon resummation has been investigated at next-to-next-to-leading logarithmic accuracy [22–24].On the other hand, there are also many attempts to gobeyond the infinite top-quark mass approximation. Thefull top-quark mass dependence was first included in thereal-emission part at NLO [14, 25]. The NLO virtualcorrections, involving multi-scale two-loop integrals sincethe LO is already a loop-induced process, have been eval-uated by expansion in the heavy top-quark mass limitup to O (1 /m t ) [26–28], in the small top-quark masslimit [29, 30], and in terms of a small Higgs transversemomentum [31] or a small Higgs mass [32]. Recently,the expansion of the three-loop virtual corrections in theheavy top-quark mass limit has been presented [33]. Fi-nally, the full NLO QCD corrections including exact de-pendence on the top-quark mass were computed numer-ically by two groups [17, 34–36], either by using a quasi-Monte Carlo method [37, 38] or via a direct Monte Carlointegration by Vegas . The matching to parton showershas also been carried out [39–41].Although the infinite top-quark mass approximation isusually insufficient for the corresponding phenomenologystudies, a standard way of improving the theoretical pre-diction on the ggF di-Higgs cross section is to use thelower-order result with full top-quark mass dependence,and to augment it with higher-order corrections in theinfinite top-quark mass limit [42].In this Letter, we provide the first next-to-next-to-next-to-leading order (N LO) perturbative QCD predic-tions for the Higgs boson pair production via gluon fusionat hadron colliders in the infinite top-quark mass limit.This result becomes one of a few highest-precision com-putations for scattering processes relevant at the LHC.The existing calculations performed at N LO include the a r X i v : . [ h e p - ph ] M a r inclusive cross sections of the ggF [43, 44], VBF [45] andbottom-quark fusion [46] of single Higgs boson produc-tion, as well as the VBF of di-Higgs production [47].Some differential distributions approximated at the sameorder are also known for the ggF of single Higgs bosonproduction [48–50]. In our calculation, we will provideboth the inclusive cross sections and the invariant-massdistribution of the Higgs boson pair at N LO, where thelatter is the first exact N LO differential distribution forthe ggF channel. After the completion of this work,the N LO corrections to the Drell-Yan processes are pre-sented in [51], where the invariant mass distribution ofthe lepton pair is obtained at N LO.
THEORETICAL FRAMEWORK
In the infinite top-quark mass limit, the effectiveLagrangian describing the coupling of two gluon fieldstrength tensors with one or two Higgs bosons reads L eff = − G aµν G a µν (cid:18) C h hv − C hh h v (cid:19) , (1)where the vacuum expectation value of the Higgs field v can be related to the Fermi constant by v = (cid:0) √ G F (cid:1) − / . C h and C hh are the Wilson coefficients by matching thefull theory to the effective theory, which start from O ( α s )and have been calculated up to O ( α s ) [52–60]. ( a ) ( b ) gg hh ( c ) FIG. 1: Representative Born cut-diagrams for the Higgs bo-son pair production in ggF. The bullets denote the vertices de-scribed by the effective Lagrangian in Eq.(1) and the crossedcircle represents the trilinear Higgs self-coupling.
The ggF di-Higgs cross section can be organized ac-cording to the number of the effective vertex insertions inthe amplitude squared, where three representative Borncut-diagrams are shown in Fig. 1. They are the ones withtwo, three and four effective vertices, and are denoted asclass- a , - b and - c respectively in the following context.Accordingly, the (differential) cross section can be splitinto three parts, dσ hh = dσ ahh + dσ bhh + dσ chh . (2)Their contributions to various α s orders are tabulated inTable I. At N LO in α s , we need to calculate the class- a (class- b and class- c ) contribution at N LO a (NNLO b and NLO c ). Here the subscripts denote the perturbativeorders in the specified class.Because of the similar topologies, the class- a part canbe obtained from the calculation of a single Higgs boson LO NLO NNLO N LOTotal O ( α s ) O ( α s ) O ( α s ) O ( α s )a LO a NLO a NNLO a N LO a O ( α s ) O ( α s ) O ( α s ) O ( α s )b — LO b NLO b NNLO b O ( α s ) O ( α s ) O ( α s )c — — LO c NLO c O ( α s ) O ( α s )TABLE I: The perturbative orders in α s for different classesat the amplitude squared level. The first and second rowsshow the perturbative expansion orders of the cross section ofHiggs pair production, while the individual class possesses itsown expansion orders that are specified by the subscript. Forexample, we call the O ( α s ) contribution in class- b as the LO b in this class since it is the first nonvanishing term, thoughit is an NLO correction to the cross section of Higgs pairproduction. production. They are related by dσ ahh dm hh = f h → hh (cid:18) C hh C h − λv m hh − m h (cid:19) × σ h ( m h → m hh ) , (3)where λ is the Higgs self-coupling and the function f h → hh accounts for the phase space difference between the singleand double Higgs boson production, f h → hh = (cid:112) m hh − m h π v , (4)and σ h ( m h → m hh ) denotes the cross section calculatedusing iHixs2 [61] after replacing the Higgs boson masswith the invariant mass of the Higgs boson pair in thecode. Such a method has also been used in the ear-lier NNLO a calculation of the ggF di-Higgs productionin Ref. [19].The class- b part can be obtained through the q T -subtraction method [62], in which we divide the crosssection into two parts, dσ bhh = dσ bhh (cid:12)(cid:12)(cid:12) p hhT
p veto T , (5)where p hhT represents the transverse momentum of theHiggs pair system. An artificial cutoff parameter p veto T is introduced in order to deal with the infrared diver-gences. In the first part, the transverse momentum ofthe Higgs pair system is required to be less than the cut-off parameter. At higher orders, there can be additionalsoft emissions or collinear emissions along the beam linein the final state with a small transverse momentum.They cause complex infrared divergences, which will can-cel against those in the virtual corrections at the end. In-stead of calculating them directly, we employ the resultin transverse momentum resummation formula. Witha sufficiently small p veto T , safely ignoring all the power-suppressed terms in p veto T , the cross section admits afactorization form that can resum the large logarithmsln n ( m hh /p veto T ) to all orders in α s . Making use of the soft-collinear effective theory [63–67], we write the factorizedcross section as a convolution of the transverse momen-tum dependent (TMD) beam function, soft function andhard function [68]. The rapidity divergences appearingin the calculation of the TMD beam function and the softfunction need an additional regulator besides the dimen-sional regularization [69–72]. However, the final physicalcross section is independent of such a regulator. The two-loop analytical results for these ingredients can be foundin [73–77]. The NNLO hard function can be obtainedby combining the two-loop amplitudes calculated in [78]and one-loop amplitudes we calculate analytically, andhave been expressed in terms of multiple polylogarithms,which can be evaluated by the public Mathematica pack-age
PolyLogTools [79]. Consequently, the first term onthe right hand of Eq.(5) is obtained by expansion of theresummation formula to a fixed order of α s . We haveset up a streamline to combine the various componentstogether in the computations of the NNLO differentialcross sections of W hh [80] and
Zhh [81] associated pro-duction processes. As opposed to the quark anti-quarkinitial states in the previous calculations, we extend ourprogram to the gluon-gluon initial states in this work.In the second part of class- b in Eq.(5), the transversemomentum of the Higgs pair system is imposed to belarger than the cutoff parameter p veto T . In such a case,there must be an additional jet accompanying with theHiggs pair. Therefore, in order to have NNLO b crosssection of class- b , we only need to calculate the NLOcorrections to hh plus a jet, of which the underlyingBorn is represented for example by Fig.1(b) but withan additional gluon emission. In this work, we use the MadGraph5 aMC@NLO [82] framework to perform suchcalculations. The two Wilson coefficients are also ex-panded in a series of α s . Since the contribution of thisclass is from the interference between the amplitudes withonly one effective vertex insertion and with two effectivevertices, one has to organize these coefficients and am-plitudes in an appropriate way. We have applied, forthe first time, the framework in [83] that was originallyproposed to handle mixed-coupling scenarios to our casewith a single coupling, α s , but with higher dimensionaleffective operators, and obtained the results order-by-order in α s . To calculate the one-loop amplitudes au-tomatically, we construct the model files by using Feyn-Rules [84],
FeynArts [85] and an in-house
Mathemat-ica program, which has been validated in [86, 87]. Wehave derived the counter-terms in this model, especiallythe rational R terms, which are needed for automaticcomputation of the virtual corrections. We have checkedpart of them extensively with the results in the litera-ture [88, 89], while the other part for the vertices with twoeffective operators are new. As a consequence, the ten-sor integrals appearing in the one-loop amplitudes, which contain five-point rank-seven integrals, can then be eval-uated by MadLoop [82, 90] equipped with
Collier [91],while the real emission contribution is computed withthe module
MadFKS [92, 93] with the FKS subtractionmethod [94, 95]. We want to stress that the inclusion ofthe contribution from class- b is indispensable in the sensethat it not only contributes to the same order in α s butalso cancels the remaining scale dependence in class- a atN LO; see the appendix.Finally, since the NLO cross sections of class- c can beobtained with full-fledged methods, we refrain ourselvesfrom presenting details about them, but they have beenroutinely included in our final results. σ b NN L O b [f b ] p Tveto [GeV] pp → hh+X √ s=13 TeVPDF4LHC15_nnlo_30 µ R = µ F =m hh /2m h =125 GeV p Thh >p Tveto p Thh
FIG. 2: The p veto T dependence of the NNLO b cross section forclass- b at the 13 TeV LHC. The error bars denote the MonteCarlo integration uncertainties. We have performed many cross checks and validationsin our calculations. All the terms except for O ( α s ) termsof class- a and class- b listed in Table I have been crosschecked by at least two independent calculations at theinclusive total cross section level. Specifically, we havereproduced the cross section of a single Higgs boson pro-duction up to NNLO in iHixs2 by using our program.This agreement can check our implementations of thetwo-loop beam and soft functions, as well as the calcula-tion of one-loop amplitudes with one effective vertex. Inaddition, we have calculated the NLO and NNLO correc-tions to Higgs pair production in the infinite top-quarkmass limit, and found agreement with HPair2 [15, 16]and Ref.[21], respectively. This helps to check Eq.(3) andour calculation of one-loop amplitudes with two effectivevertices. These nontrivial checks already ensure the cor-rectness of many components of our calculations. Forthe O ( α s ) term of class- a , we simply used iHixs2 by em-ploying Eq.(3). Such a program has been validated withthe Higgs pair cross sections from LO to NNLO, whichmakes us convinced that the O ( α s ) piece of class- a is cor-rect. For the remaining O ( α s ) part of class- b , we carefullychecked the various pieces that are used in our calcula-tion. In particular, we have checked the scale depen-dence of the finite part in the two-loop amplitudes withtwo effective vertices [78] by the renormalization groupequation that the hard function should satisfy. The one-loop amplitude can also been extracted from the scale-dependent part of the two-loop amplitudes, and it hasbeen compared against the analytical result we calculatedwith the assistance of fire [96] and to the numerical re-sult from MadLoop . Again, we find perfect agreements.Moreover, we have checked the independence of the finalNNLO b results for class- b on the values of p veto T over therange from 4 GeV to 20 GeV, as shown in Fig.2. NUMERICAL RESULTS
In our numerical calculations, we take v = 246 . m h = 125 GeV. The top-quarkpole mass, which enters only into the Wilson coefficients,is m t = 173 . PDF4LHC15 nnlo 30
PDF [97–100] provided by
LHAPDF6 [101], and the as-sociated strong coupling α s . The default central scale ischosen to be the invariant mass of the Higgs pair dividedby 2, i.e. µ = m hh /
2, and the scale uncertainty is eval-uated through the 9-point variation of the factorizationscale µ F and the renormalization scale µ R in the form of µ R,F = ξ R,F µ with ξ R , ξ F ∈ { . , , } . Order √ s
13 TeV 14 TeV 27 TeV 100 TeVLO 13 . +31% − . +31% − . +26% − +19% − NLO 25 . +18% − . +18% − . +16% − +13% − NNLO 30 . +5 . − . . +5 . − . . +4 . − . +4 . − . N LO 31 . +0 . − . . +0 . − . . +0 . − . +0 . − . TABLE II: The inclusive total cross sections (in unit of fb)of Higgs boson pair production at different center-of-mass en-ergies from LO to N LO. The quoted relative uncertaintiesare from the 9-point scale variations µ R,F = ξ R,F m hh / ξ R , ξ F ∈ { . , , } . The errors due to the numerical MonteCarlo integration are well below 1 (cid:104) . We present the inclusive total cross sections (from LOto N LO) of the Higgs boson pair production at differentcenter-of-mass energies in Table II and Fig. 3. Similarlyto the single Higgs case, the QCD higher-order correc-tions are prominent. The NLO corrections increase theLO cross section by 87% (85%) at √ s = 13 (100) TeV.The NNLO corrections increase the NLO cross sectionfurther by 18% (16%), reducing the scale uncertainty bya factor of 2 to 3 to be below 8%. Finally, the N LOcorrections turn out to be 3 .
0% (2 . LO is less than 3% (2%),with another significant reduction of 2-3 times. For thepurpose of the comparison, the PDF parameterizationuncertainty at 13 TeV amounts to ± . σ [f b ] pp → hh+XPDF4LHC15_nnlo_30m hh /4< µ R , µ F 3, the perturbative corrections to the inclusivecross section is small from NLO to N LO.Besides the inclusive total cross section, we are alsoable to obtain the exact N LO results for a differentialdistribution, i.e., the invariant mass m hh distributionshown in Fig.5. As in the total cross section case, theinclusion of the N LO corrections dramatically stabilizesthe perturbative calculation of the invariant mass differ-ential distribution. It can also be seen that the higher-order QCD corrections do not change the peak position,and the K factor of N LO over NNLO is almost flat overa large region of m hh . The N LO result with small scaleuncertainty is completely enclosed within the NNLO un-certainty band. Such a feature consolidates that the per-turbative expansion of this differential cross section in aseries of α s is converging up to this order. σ [f b ] log ( ξ ) pp → hh+X √ s=13 TeVPDF4LHC15_nnlo_30 µ R = µ F = ξ m hh /2m h =125 GeV LONLONNLON LO FIG. 4: The scale dependence of the total cross section forHiggs boson pair production at the LHC with √ s = 13 TeV.We set µ R = µ F = ξµ in this plot. d σ / d m hh [f b / G e V ] pp → hh+X √ s=13 TeVPDF4LHC15_nnlo_30m hh /4< µ R , µ F CONCLUSIONS We have calculated the N LO QCD corrections tothe Higgs boson pair production via gluon fusion athadron colliders in the infinite top-quark mass limit. Wefind that the total cross section at N LO increases by3 . 0% (2 . √ s = 13 (100) TeV with respect to theNNLO result under the central scale choice µ = m hh / LO compared to the previous result at NNLO, whichis now less than 3% (2%). In contrast, the PDF uncer-tainty is ± . 3% at the 13 TeV LHC. Moreover, we havecomputed the invariant mass distribution at N LO forthe first time, and the shape of the distribution is almostunchanged. The perturbative series of both the total in-clusive cross section and the invariant mass distributionare found to be converging up to this order.In the future, for the phenomenological applications,it is essential to combine our N LO calculation in theinfinite top-quark mass limit with the NLO result in-cluding exact top-quark mass dependence [102]. TheNNLO result with finite top-quark mass is not availableat the moment. A preliminary investigation indicatesthat the uncertainty from taking infinite top-quark masslimit at NNLO can be as large as ± 3% ( ± Acknowledgments We thank D. Y. Shao for collaborations at the earlystage of this work. We also thank C. Lee and Y.-Q. Mafor carefully reading the manuscript. LBC was supportedby the National Natural Science Foundation of China(NSFC) under the grants 11747051 and 11805042. HTLwas supported by the Los Alamos National LaboratoryLDRD program. The work of HSS was supported by theILP Labex (ANR-11-IDEX-0004-02, ANR-10-LABX-63).JW was supported by the BMBF project 05H18WOCA1when he was in Technische Universit¨at M¨unchen and bythe program for Taishan scholars. Appendix: Renormalization scale dependence In this appendix, we present the method to obtain the µ R dependence in the framework of soft-collinear effectivetheory (SCET), where the renormalization scale is usually set to be the same as the factorization scale in the expansionof the resummed formula, while they can be distinguished clearly in fixed-order calculations.The cross section is scale ( µ R = µ F = µ ) invariant, dd ln µ σ hh ( µ, µ ) = (cid:18) dd ln µ R σ hh ( µ R , µ F ) + dd ln µ F σ hh ( µ R , µ F ) (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) µ R = µ F = µ = 0 + O ( α s ) . (6)The individual renormalization and factorization scale dependence can be obtained through σ hh ( µ R , µ F ) = σ hh ( µ F , µ F ) + (cid:90) µ R µ F d ¯ µ (cid:32) ddµ R σ hh ( µ R , µ F ) (cid:12)(cid:12)(cid:12)(cid:12) µ R =¯ µ (cid:33) , (7)where the first part on the right hand can be predicted with q T -subtraction method in the framework of SCET andthe second part is obtained by requiring the renormalization scale independence of the total cross section.The N LO cross section for Higgs pair production is renormalization scale invariant up to O ( α s ) corrections, i.e. dd ln µ R σ hh ( µ R , µ F ) = dd ln µ R σ ahh ( µ R , µ F ) + dd ln µ R σ bhh ( µ R , µ F ) + dd ln µ R σ chh ( µ R , µ F ) = 0 + O ( α s ) . (8)For class- a , the differential equation is dd ln µ R σ ahh ( µ R , µ F ) = (cid:90) dm hh f h → hh [ σ h ( µ R , µ F , m h → m hh )] × dd ln µ R (cid:18) C hh C h − λv m hh − m h (cid:19) . (9)where σ h has the expansion σ h = σ (0) h + σ (1) h + . . . with σ ( i ) h ∝ α is . Up to N LO, we need the NLO QCD correctionsto class- c cross section which is standalone and scale invariant. Therefore, for class- b , the renormalization groupequation is dd ln µ R σ bhh ( µ R , µ F ) = − (cid:90) dm hh f h → hh [ σ h ( µ R , µ F , m h → m hh )] × (cid:18) C hh C h − λv m hh − m h (cid:19) (cid:18) dd ln µ R C hh C h (cid:19) . (10)The ratio of C hh over C h can be written in a series of a s ≡ α s ( µ R ) / π , C hh C h = 1 + δ a s + δ ( µ R ) a s + O ( a s ) . (11)where the coefficient δ is scale independent. Therefore, we obtain dd ln µ R C hh C h = (cid:18) da s d ln µ R ∂∂a s + ∂∂ ln µ R (cid:19) C hh C h = − β δ a s + a s dδ d ln µ R + O ( a s ) ≡ a s η + O ( a s ) (12)with β = (11 C A − n f ) / dd ln µ R σ bhh ( µ R , µ F ) = − a s η (cid:90) dm hh f h → hh [ σ (0) h ( µ R , µ F , m h → m hh )] (cid:18) − λv m hh − m h (cid:19) + O ( a s )= − a s η σ b (1) hh ( µ R , µ F ) + O ( a s ) , (13)where the class- b cross section is written as σ bhh = (cid:80) i =1 a is σ b ( i ) hh with σ b ( i ) hh ∝ a s and σ b (1) hh ( µ R , µ F ) = 83 (cid:90) dm hh f h → hh [ σ (0) h ( µ R , µ F , m h → m hh )] (cid:18) − λv m hh − m h (cid:19) . (14)The scale-invariance violation terms in Eq. (13) start from NNLO corrections to class- b Higgs pair production. ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] [1] G. Aad et al. (ATLAS Collaboration), Phys.Lett. B716 , 1 (2012), 1207.7214.[2] S. Chatrchyan et al. (CMS), Phys. Lett. B716 , 30 (2012), 1207.7235.[3] G. Aad et al. (ATLAS), Eur. Phys. J. C76 , 6 (2016), 1507.04548.[4] G. Aad et al. (ATLAS, CMS), JHEP , 045 (2016), 1606.02266.[5] A. M. Sirunyan et al. (CMS), Eur. Phys. J. C79 , 421 (2019), 1809.10733.[6] G. Aad et al. (ATLAS) (2019), 1909.02845.[7] A. M. Sirunyan et al. (CMS), Phys. Rev. Lett. , 121803 (2019), 1811.09689.[8] G. Aad et al. (ATLAS) (2019), 1906.02025.[9] M. Cepeda et al. (HL/HE WG2 group) (2019), 1902.00134.[10] R. Contino et al., CERN Yellow Rep. pp. 255–440 (2017), 1606.09408.[11] J. H. Kim, K. Kong, K. T. Matchev, and M. Park, Phys. Rev. Lett. , 091801 (2019), 1807.11498.[12] S. Kanemura, S. Kiyoura, Y. Okada, E. Senaha, and C. P. Yuan, Phys. Lett. B558 , 157 (2003), hep-ph/0211308.[13] J. Baglio, A. Djouadi, R. Groeber, M. M. Muehlleitner, J. Quevillon, and M. Spira, JHEP , 151 (2013), 1212.5581.[14] R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, P. Torrielli, E. Vryonidou, and M. Zaro, Phys. Lett. B732 ,142 (2014), 1401.7340.[15] S. Dawson, S. Dittmaier, and M. Spira, Phys. Rev. D58 , 115012 (1998), hep-ph/9805244.[16] T. Plehn, M. Spira, and P. M. Zerwas, Nucl. Phys. B479 , 46 (1996), [Erratum: Nucl. Phys.B531,655(1998)], hep-ph/9603205.[17] J. Baglio, F. Campanario, S. Glaus, M. Muehlleitner, M. Spira, and J. Streicher, Eur. Phys. J. C79 , 459 (2019),1811.05692.[18] D. de Florian and J. Mazzitelli, Phys. Lett. B724 , 306 (2013), 1305.5206.[19] D. de Florian and J. Mazzitelli, Phys. Rev. Lett. , 201801 (2013), 1309.6594.[20] J. Grigo, K. Melnikov, and M. Steinhauser, Nucl. Phys. B888 , 17 (2014), 1408.2422.[21] D. de Florian, M. Grazzini, C. Hanga, S. Kallweit, J. M. Lindert, P. Maierhoefer, J. Mazzitelli, and D. Rathlev, JHEP , 151 (2016), 1606.09519.[22] D. Y. Shao, C. S. Li, H. T. Li, and J. Wang, JHEP , 169 (2013), 1301.1245.[23] D. de Florian and J. Mazzitelli, JHEP , 053 (2015), 1505.07122.[24] D. De Florian and J. Mazzitelli, JHEP , 156 (2018), 1807.03704.[25] F. Maltoni, E. Vryonidou, and M. Zaro, JHEP , 079 (2014), 1408.6542.[26] J. Grigo, J. Hoff, K. Melnikov, and M. Steinhauser, Nucl. Phys. B875 , 1 (2013), 1305.7340.[27] J. Grigo, J. Hoff, and M. Steinhauser, Nucl. Phys. B900 , 412 (2015), 1508.00909.[28] G. Degrassi, P. P. Giardino, and R. Groeber, Eur. Phys. J. C76 , 411 (2016), 1603.00385.[29] J. Davies, G. Mishima, M. Steinhauser, and D. Wellmann, JHEP , 048 (2018), 1801.09696.[30] J. Davies, G. Mishima, M. Steinhauser, and D. Wellmann, JHEP , 176 (2019), 1811.05489.[31] R. Bonciani, G. Degrassi, P. P. Giardino, and R. Groeber, Phys. Rev. Lett. , 162003 (2018), 1806.11564.[32] X. Xu and L. L. Yang, JHEP , 211 (2019), 1810.12002.[33] J. Davies and M. Steinhauser (2019), 1909.01361.[34] S. Borowka, N. Greiner, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk, U. Schubert, and T. Zirke, Phys. Rev. Lett. , 012001 (2016), [Erratum: Phys. Rev. Lett.117,no.7,079901(2016)], 1604.06447.[35] S. Borowka, N. Greiner, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk, and T. Zirke, JHEP , 107 (2016), 1608.04798.[36] J. Davies, G. Heinrich, S. P. Jones, M. Kerner, G. Mishima, M. Steinhauser, and D. Wellmann (2019), 1907.06408.[37] Z. Li, J. Wang, Q.-S. Yan, and X. Zhao, Chin. Phys. C40 , 033103 (2016), 1508.02512.[38] J. Dick, Journal of Complexity , 493 (2004), ISSN 0885-064X.[39] G. Heinrich, S. P. Jones, M. Kerner, G. Luisoni, and E. Vryonidou, JHEP , 088 (2017), 1703.09252.[40] S. Jones and S. Kuttimalai, JHEP , 176 (2018), 1711.03319.[41] G. Heinrich, S. P. Jones, M. Kerner, G. Luisoni, and L. Scyboz, JHEP , 066 (2019), 1903.08137.[42] M. Grazzini, G. Heinrich, S. Jones, S. Kallweit, M. Kerner, J. M. Lindert, and J. Mazzitelli, JHEP , 059 (2018),1803.02463.[43] C. Anastasiou, C. Duhr, F. Dulat, F. Herzog, and B. Mistlberger, Phys. Rev. Lett. , 212001 (2015), 1503.06056.[44] B. Mistlberger, JHEP , 028 (2018), 1802.00833.[45] F. A. Dreyer and A. Karlberg, Phys. Rev. Lett. , 072001 (2016), 1606.00840.[46] C. Duhr, F. Dulat, and B. Mistlberger (2019), 1904.09990.[47] F. A. Dreyer and A. Karlberg, Phys. Rev. D98 , 114016 (2018), 1811.07906.[48] F. Dulat, B. Mistlberger, and A. Pelloni, JHEP , 145 (2018), 1710.03016.[49] L. Cieri, X. Chen, T. Gehrmann, E. W. N. Glover, and A. Huss, JHEP , 096 (2019), 1807.11501.[50] F. Dulat, B. Mistlberger, and A. Pelloni, Phys. Rev. D99 , 034004 (2019), 1810.09462.[51] C. Duhr, F. Dulat, and B. Mistlberger (2020), 2001.07717.[52] T. Inami, T. Kubota, and Y. Okada, Z. Phys. C18 , 69 (1983).[53] K. G. Chetyrkin, B. A. Kniehl, and M. Steinhauser, Phys. Rev. Lett. , 353 (1997), hep-ph/9705240.[54] K. G. Chetyrkin, B. A. Kniehl, and M. Steinhauser, Nucl. Phys. B510 , 61 (1998), hep-ph/9708255.[55] Y. Schroder and M. Steinhauser, JHEP , 051 (2006), hep-ph/0512058.[56] K. G. Chetyrkin, J. H. Kuhn, and C. Sturm, Nucl. Phys. B744 , 121 (2006), hep-ph/0512060.[57] B. A. Kniehl, A. V. Kotikov, A. I. Onishchenko, and O. L. Veretin, Phys. Rev. Lett. , 042001 (2006), hep-ph/0607202.[58] P. A. Baikov, K. G. Chetyrkin, and J. H. K¨uhn, Phys. Rev. Lett. , 082002 (2017), 1606.08659. [59] M. Spira, JHEP , 026 (2016), 1607.05548.[60] M. Gerlach, F. Herren, and M. Steinhauser, JHEP , 141 (2018), 1809.06787.[61] F. Dulat, A. Lazopoulos, and B. Mistlberger, Comput. Phys. Commun. , 243 (2018), 1802.00827.[62] S. Catani and M. Grazzini, Phys. Rev. Lett. , 222002 (2007), hep-ph/0703012.[63] C. W. Bauer, S. Fleming, and M. E. Luke, Phys. Rev. D63 , 014006 (2000), hep-ph/0005275.[64] C. W. Bauer, S. Fleming, D. Pirjol, and I. W. Stewart, Phys.Rev. D63 , 114020 (2001), hep-ph/0011336.[65] C. W. Bauer and I. W. Stewart, Phys. Lett. B516 , 134 (2001), hep-ph/0107001.[66] C. W. Bauer, D. Pirjol, and I. W. Stewart, Phys.Rev. D65 , 054022 (2002), hep-ph/0109045.[67] M. Beneke, A. Chapovsky, M. Diehl, and T. Feldmann, Nucl.Phys. B643 , 431 (2002), hep-ph/0206152.[68] T. Becher, M. Neubert, and D. Wilhelm, JHEP , 110 (2013), 1212.2621.[69] J.-y. Chiu, A. Jain, D. Neill, and I. Z. Rothstein, Phys. Rev. Lett. , 151601 (2012), 1104.0881.[70] J.-Y. Chiu, A. Jain, D. Neill, and I. Z. Rothstein, JHEP , 084 (2012), 1202.0814.[71] T. Becher and G. Bell, Phys. Lett. B713 , 41 (2012), 1112.3907.[72] Y. Li, D. Neill, and H. X. Zhu, Submitted to: Phys. Rev. D (2016), 1604.00392.[73] T. Gehrmann, T. Lubbert, and L. L. Yang, Phys. Rev. Lett. , 242003 (2012), 1209.0682.[74] T. Gehrmann, T. Luebbert, and L. L. Yang, JHEP , 155 (2014), 1403.6451.[75] T. Luebbert, J. Oredsson, and M. Stahlhofen, JHEP , 168 (2016), 1602.01829.[76] M. G. Echevarria, I. Scimemi, and A. Vladimirov, JHEP , 004 (2016), 1604.07869.[77] M.-X. Luo, X. Wang, X. Xu, L. L. Yang, T.-Z. Yang, and H. X. Zhu (2019), 1908.03831.[78] P. Banerjee, S. Borowka, P. K. Dhani, T. Gehrmann, and V. Ravindran, JHEP , 130 (2018), 1809.05388.[79] C. Duhr and F. Dulat (2019), 1904.07279.[80] H. T. Li and J. Wang, Phys. Lett. B765 , 265 (2017), 1607.06382.[81] H. T. Li, C. S. Li, and J. Wang, Phys. Rev. D97 , 074026 (2018), 1710.02464.[82] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro,JHEP , 079 (2014), 1405.0301.[83] R. Frederix, S. Frixione, V. Hirschi, D. Pagani, H. S. Shao, and M. Zaro, JHEP , 185 (2018), 1804.10017.[84] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, Comput. Phys. Commun. , 2250 (2014), 1310.1921.[85] T. Hahn, Comput. Phys. Commun. , 418 (2001), hep-ph/0012260.[86] H.-S. Shao, Y.-J. Zhang, and K.-T. Chao, JHEP , 048 (2011), 1106.5030.[87] H.-S. Shao and Y.-J. Zhang, JHEP , 112 (2012), 1205.1273.[88] P. Draggiotis, M. V. Garzelli, C. G. Papadopoulos, and R. Pittau, JHEP , 072 (2009), 0903.0356.[89] B. Page and R. Pittau, JHEP , 078 (2013), 1307.6142.[90] V. Hirschi, R. Frederix, S. Frixione, M. V. Garzelli, F. Maltoni, and R. Pittau, JHEP , 044 (2011), 1103.0621.[91] A. Denner, S. Dittmaier, and L. Hofer, Comput. Phys. Commun. , 220 (2017), 1604.06792.[92] R. Frederix, S. Frixione, F. Maltoni, and T. Stelzer, JHEP , 003 (2009), 0908.4272.[93] R. Frederix, S. Frixione, A. S. Papanastasiou, S. Prestel, and P. Torrielli, JHEP , 027 (2016), 1603.01178.[94] S. Frixione, Z. Kunszt, and A. Signer, Nucl. Phys. B467 , 399 (1996), hep-ph/9512328.[95] S. Frixione, Nucl. Phys. B507 , 295 (1997), hep-ph/9706545.[96] A. V. Smirnov, Comput. Phys. Commun. , 182 (2015), 1408.2372.[97] J. Butterworth et al., J. Phys. G43 , 023001 (2016), 1510.03865.[98] S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump, and C. P. Yuan, Phys.Rev. D93 , 033006 (2016), 1506.07443.[99] L. A. Harland-Lang, A. D. Martin, P. Motylinski, and R. S. Thorne, Eur. Phys. J.