Top-quark effects in diphoton production through gluon fusion at NLO in QCD

At hadron colliders, the leading production mechanism for a pair of photons is from quark-anti-quark annihilation at the tree level. However, due to large gluon-gluon luminosity, the loop-induced process $gg\to \gamma \gamma$ provides a substantial contribution. In particular, the amplitudes mediated by the top quark become important at the $t \bar t$ threshold and above. In this letter we present the first complete computation of the next-to-leading order (NLO) corrections (up to $\alpha_S^3$) to this process, including contributions from the top quark. These entail two-loop diagrams with massive propagators whose analytic expressions are unknown and have been evaluated numerically. We find that the NLO corrections to the top-quark induced terms are very large at low diphoton invariant mass $m(\gamma \gamma)$ and close to the $t \bar t$ threshold. The full result including five massless quarks and top quark contributions at NLO displays a much more pronounced change of slope in the $m(\gamma \gamma)$ distribution at $t \bar t$ threshold than at LO and an enhancement at high invariant mass with respect to the massless calculation.

Diphoton production is one important channel at hadron colliders. Gluon fusion into diphoton via quark loops provide a substantial contribution. When the top quark mediated contribution is negligible below the top pair threshold, it becomes important at the threshold and above. We present the first complete computation of the next-to-leading order(NLO) corrections to this process, including contributions from the top quark. The analytical expression for relevant two-loop diagrams are unknown and we developed a numerical algorithm to compute them. We found that NLO corrections is large over all kinematic regions. The top quark contribution leads to change of slope below and above the top pair threshold, and it is more visible at NLO. We further find that above the top quark threshold, NLO corrections is larger after including the top quark contribution.

Introduction
The production of a pair of photons (diphoton) is an important process at hadron colliders. The final state signature is experimentally clean, and experimental measurements are available at Tevatron [1,2] and LHC [3,4]. It is one golden channel for the discovery of the Higgs boson [5,6], and H → γγ remains one of the cleanest final states to study the properties of the Higgs boson. Diphoton channel is also very important channel on searching new physics beyond the Standard Model [7,8], including new scalar or spin-2 reasonance, multiple resonances from extra-dimension/clockwork models [9,10,11,12], or peak-dip structures due to new particles in loops [13].
At leading order, it is produced through quark-antiquark annihilation. Next-to-Leading order(NLO) corrections to this process has been known since 2000 [14]. Recently, NNLO corrections are also know and available in public code [15,16,17]. In addition to quark-antiquark annihilation, diphoton can also be produced through gluon fusion. Formally it can be counted as part of NNLO corrections to qq → γγ. However, due to large gluon-gluon luminosity, the contribution from gluon fusion is anomalous large. Since the gluon fusion process is separately gauge-invariant and IR finite, it can be treated as a standalone channel, and thus NLO corrections can be defined and calculated, without consider full NNNLO corrections to qq → γγ. However, at NLO, only the contribution from massless quarks are known [18,16]; the contribution from the top quark is still missing. Naively counting electric charge, after including top quark the cross section can be 1.86 times larger. Therefore, including the contribution from top quark is essential for precise prediction. The cross section for gg → γγ is given by

Method
where dσ LO is the LO cross section, dσ R is the corrections with one more external legs(gg → γγg, q(q)g → γγq(q), qq → γγg), dσ V is the corrections with one more loop, and dσ C is the counter term from the parton distribution function. We show representative Feynman diagrams in Fig. 1.
The last three terms are IR divergent, and only their sum is finite. To handle the IR divergence, we adopt the dipole subtraction method [19], which introduces extra terms for the last three terms so that they become finite, while the extra terms sum up to zero.
With the subtraction method available, we also need to calculate the matrix elements. For the real corrections, since they are one-loop process, automated tools are available to compute them. In particular, we adopt MADGRAPH5_AMC@NLO [20] and Recola2 [21,22] to compute those matrix elements.
The main challenge is on the calculation of the two-loop virtual amplitude. In the case that internal lines are massless quarks, the results are known since 2001 [23]. However, for massive quark contribution, the results remain unknown. The main reason is that not all master integrals are known analytically; only part of them are known(mostly planar). Instead of analytical approach, we developed a numerical method to calculate them [24].
We use QGRAF [25] to generate corresponding Feynman diagrams, which yields 138 diagrams. We use FORM [26] to deal further proceed the expressions of the amplitude. We adopt projector methods to decompose the amplitudes into form factors. In particular, we work under ddimension, and we get 10 different tensor bases. We organize the Feynman integrals according to their propagator denominators, obtaining 33 different integral families. After projection, we obtain the form factors written in terms of around 40000 scalar integrals at two-loop. We adopt FIRE5 [27] with LiteRed [28,29] to perform integration-by-parts(IBP) reduction, and we get 1180 master integrals, not counting relations among families. If we consider symmetries among different families, indeed we only have 161 master integrals.
With IBP reductions, it was shown [30,31] that the derivatives of master integrals can be written in terms of linear combination of themselves, with coefficients are rational expressions on the dimension d and kinematics. Such system of differential equations provide an important way on analytical computation of master integrals. On the other hand, numerical methods for differential equations are well-studied in mathematics. Even if we have multiple kinematics, and thus such system of differential equations is a system partial differential equations by nature, we can integrate all kinematic iteratively. Therefore, numerical algorithms for ordinary differential equations can be applied. We adopt numerical methods for initial value problems to solve the differential equations, which requires numerical value of master integrals in a specific point as the initial condition to fully fix the solution.
We adopt the sector decomposition method [32] to obtain the initial condition. The sector decomposition method is a systematic method that resolving IR and UV divergence for Feynman integrals. After resolving divergences, the remaining integrals can be expanded order-by-order in ε, and numerical integration can be carried out on those integrals. When working in the Euclidean region, saying that the −i0 prescription can be discarded, the integrands are well-behaved, and numerical integration is quite efficient with suitable methods. In particular we adopt the quasi-Monte Carlo methods [33]. We implement the sector decomposition method in NIFT [34].
In Table 1, we show numerical results of some relevant Feynman integrals obtain by NIFT and compare it with analytical expression. The corresponding diagrams are shown in Fig. 2. To obtain numerical results to O(10 −7 ) precision, only several seconds are spent.
Since the initial condition is given in the unphysical region, we need to design an integration contour to evolve it to the physical region. The main requirement of the integration contour is that it should not cross branch cuts. It can be also understood that along the integration contour the Feynman integrals can be defined without the −i0 prescription, and the integration contour should connect to the target point according to such prescription. For the integrals appear in the two-loop    amplitude of gg → γγ, we show the integration contour and branch cuts in Fig. 3.
In Table 2, we compare numerical results obtained via our algorithms with the one from analytical expressions. We can see that O(10 −7 ) precision can be obtained within 1 second. In Table  3, we show results for I 2 and I 3 , whose analytical results are unknown. From the discrepancy of results from two differential initial condition, we conclude that our numerical precision is still at O(10 −7 ). For cross check, we also show results obtained from pySecDec, which is much slower and the precision is much lower.
To further speed up the numerical evaluation, we first compute the results at a set of predefined kinematic points in the physical region. After such initial condition table is setup, when evaluating the master integrals at each point, we pick a point in such set and use that result as the initial condition to compute the results. Since the initial condition belongs to the physical region now, it is much closer to the target point, and the integration contour is much simpler and shorter, which   Table 3: Numerical results obtained with our algorithm from two different choices of initial conditions for the Feynman integral I 2 and I 3 are shown. The results from pySecDec [35] are also shown for comparison.
yield a significant speedup. Furthermore, since the results of the initial condition table are only need to be computed once, we demand higher accuracy when computing them.
In addition, we observed that t is unchanged during s evolution. Therefore, we reorganize the expressions into a form which is rational expressions on ε and s, with coefficients are polynomials on t. The coefficients are unchanged when evolving the differential equations, and thus they only need to be computed once in each evolution. Similar reorganized can be done for t evolution. For better numerical stability in the backward region(|u| |t|, |s|), we use u to compute the right-hand side of differential equations. That is to say, we adopt the following forms in the backward region.
In total, 10000 CPU hours are spent to obtain the initial conditions with O(10 −9 ) precision. We expand up to O(ε) for 7-propagator master integrals, and O(ε 2 ) for other master integrals. For each phasespace point, it costs around 1 second to evaluate the two-loop amplitude.

PoS(RADCOR2019)024
Top quark effects in gg → γγ Xiaoran Zhao   (d) Differential cross section for diphoton production through gluon fusion at the 13TeV LHC in the threshold region are shown for the massless quark only case(5F only), and also including top quark contribution(full).

Results
In Fig. 4a, we show the differential cross section for 13 TeV LHC, only including the top quark contribution. Below top pair threshold, the cross section is strongly suppressed, since the amplitude scaling as O(m −4 t ). The differential cross section peaks around top pair threshold, and NLO corrections is large in all regions. In the threshold region, there is a peak in the K-factor, due to the contribution comes from one Coulomb gluon effects exchange.
Indeed, Coulomb gluon exchange yields corrections as α S /v, which may spoil perturbativity. The contribution from one Coulomb gluon exchange is already included in the two-loop amplitude, and multiple Coulomb gluon exchange is corresponding to the formation of toponium T (nS). The contribution of toponium can be included explicitly by adding the contribution of Feynman diagram(see fig. 4b which corresponding to the following: In fig. 4c, we show the differential cross section in the threshold region for the top quark only contribution. NLO corrections are quite large in such region. Furthermore, including multiple Coulomb gluon exchanges leads a peak slightly below top pair threshold, which is corresponding to the lightest toponium T (1S). Away from top pair threshold, the matched results reduce to the fixed order result.
In fig. 4d, we show the differential cross section in the threshold region for the light quark only case and after including top quark contribution. Due to negative interference effects, after including top quark the differential cross section gets reduced. In addition, below and above the top pair threshold we can see that the slope is differential, and such effect is more visible at NLO than at LO.
In the low energy region, the contribution from the top quark can be described as an EFT. Taking only the Abelian part for simplicity, which is corresponding to QED case, the EFT Lagrangian is given by: The Wilson coefficients can be obtained by matching the amplitude in the low-energy limit in the full theory and EFT. However, naively extrapolating the amplitude to the low-energy limit(s = 0) leads to low precision. Alternatively, we calculate the low energy limit based on Cauchy integral formula, which only need the amplitude with non-zero s. In Table 4 we compare the results obtain via the above approach with analytic results [36], and find that our results are agree with the one in literature, with good precision. In fig. 5, we show the differential cross section for m(γγ) ∈ [100, 1000] GeV. We can see that below top pair threshold, the top quark contribution is negligible. In the threshold region including top quark decreases the cross section, and above threshold region the cross section becomes larger after including top quark contribution. Moreover, the K-factor for the full case is larger than light quark only case. As the invariant mass increases, the ratio between full case and light quark only case slowly approaching naive six-flavour limit(1.86).

Conclusion
We consider NLO corrections to gg → γγ. We include both light quark contributions and the top quark contribution. We developed numerical methods for the two-loop massive amplitudes. We find the NLO corrections are large for the top quark contribution. We also consider resume Coulomb gluon effects in the threshold region and match it with fixed order results. We further examine the low-energy behavior and match it with EFT operators. We find that the slope change below and above the top pair threshold is more visible at NLO than at LO. Going above the top pair threshold, the differential cross section get further enhanced after including the top quark contribution, and more closer to the naive six flavour case at NLO than at LO.