Photon pair production in gluon fusion: Top quark effects at NLO with threshold matching

We present a calculation of the NLO QCD corrections to the loop-induced production of a photon pair through gluon fusion, including massive top quarks at two loops, where the two-loop integrals are calculated numerically. Matching the fixed-order NLO results to a threshold expansion, we obtain accurate results around the top quark pair production threshold. We analyse how the top quark threshold corrections affect distributions of the photon pair invariant mass and comment on the possibility of determining the top quark mass from precision measurements of the diphoton invariant mass spectrum.


Introduction
The production of pairs of photons in hadronic collisions has attracted interest from both the experimental and the theory side for several decades. Most prominently, the diphoton final state served as one of the key discovery channels for the Higgs boson [1,2], which can decay into two photons. As a very clean experimental channel, it is also well suited for precision studies of the Standard Model (SM) and in particular the Higgs sector. For example, there is the possibility to constrain the Higgs boson width from interference effects of the continuum gg → γγ spectrum with the signal gg → H → γγ [3][4][5][6][7][8][9][10]. Furthermore, various New Physics models predict the production of photon pairs, where the study of angular correlations between the decay photons can provide information about the spin of the underlying resonances [11,12]. Another interesting aspect of diphoton production is the possibility of measuring the top quark mass via the top quark pair production threshold effects manifest in the diphoton invariant mass spectrum [13,14]. While current LHC measurements [1,2] are not yet able to provide the necessary statistics for such a threshold scan, the feasibility at the High-Luminosity LHC, and even more so at a future 100 TeV collider, is worth investigating.
Direct diphoton production 1 in hadronic collisions occurs via the leading order (LO) α 0 s process qq → γγ. The next-to-leading order (NLO) QCD corrections to this process, including fragmentation contributions at NLO, were implemented in the public program Diphox [15]. The loop induced gg → γγ process enters as a next-to-next-to-leading order (NNLO) QCD (order α 2 s ) correction to the pp → γγ cross section. The process gg → γγ has been calculated at LO including both massless and massive quark loops in Ref. [3] and is included in Diphox at one loop for massless quark loops. Even though the gg → γγ contribution is a higher-order correction to the total pp → γγ cross section, its contribution is similar in size to the LO result at the LHC, due to the large gluon luminosity. A calculation that includes also the effects of transverse-momentum resummation to direct photon production is implemented in the program ResBos [16]. NLO QCD corrections to the gluon-fusion channel with massless quarks, i.e. O(α 3 s ) corrections, have been first calculated in Refs. [17,18] and implemented in the code 2γMC [18] as well as in MCFM [19]. Very recently, the NLO QCD corrections to the gluon-fusion channel including massive top quark loops have become available [20], where the master integrals have been calculated numerically based on the numerical solution of differential equations [21,22]. Analytic results for the planar two-loop box integrals with massive top quarks have been presented in Ref. [23,24]. Regarding the non-planar contributions, 3-point topologies containing elliptic integrals have been calculated in Ref. [25,26]. Other 3-point topologies have been calculated earlier in the context of Higgs production and decay [27,28]. The NNLO QCD corrections to the process pp → γγ were first calculated in Ref. [29], including the gg → γγ contribution at order α 2 s with massless quark loops. For a phenomenological study see also Ref. [30]. The NNLO QCD corrections to pp → γγ have also been calculated and implemented in MCFM in Ref. [31], supplemented by the gg initiated loops proportional to n f at LO and NLO for five massless quark flavours, and at LO for massive top quark loops. Diphoton production at NNLO with massless quarks is also available in Matrix [32]. The aim of this paper is twofold. Firstly, we provide an independent calculation of the QCD corrections to the process gg → γγ including massive top quark loops, confirming the results of Ref. [20] for the central scale choice. Secondly, we combine our results with threshold resummation as advocated in Ref. [14], such that the top quark pair production threshold region in the diphoton invariant mass spectrum can be predicted with high accuracy. The calculation can thus serve as a starting point for investigating the possibility of a top quark mass measurement from the diphoton invariant mass spectrum. This work is structured as follows. In Section 2 we describe our calculation of the NLO corrections including both massless and massive fermion loops. Section 3 contains a description of our treatment of the top quark pair production threshold region. In Section 4 we present our numerical results. Finally, in Section 5 we summarise and present an outlook on the possibility of measuring the top quark mass from the diphoton spectrum.

Building blocks of the fixed order calculation
We consider the following scattering process, with on-shell conditions p 2 j = 0, j = 1, ..., 4. The helicities λ i of the external particles are defined by taking the momenta of the gluons p 1 and p 2 (with colour indices a 1 and a 2 , respectively) as incoming and the momenta of the photons p 3 and p 4 as outgoing. The Mandelstam invariants associated with eq. (2.1) are defined by

Projection operators
We define the tensor amplitude M µ 1 µ 2 µ 3 µ 4 by extracting the polarisation vectors from the amplitude M, where the ε µ i λ i denote the polarisation vectors. The amplitude is computed through projection onto a set of Lorentz structures related to linear polarisation states of the external massless bosons. An appropriate set of D-dimensional projection operators is constructed following the approach proposed in Ref. [33], which has been applied recently in the calculation of Ref. [34], and which we will summarise briefly in the following. A physical polarisation vector ε(p) of a massless vector boson with (on-shell) momentum p fulfils the transversality and (imposed) normalisation conditions, ε(p) · p = 0, ε(p) · ε(p) = −1. (2.4) These conditions fix two components of the polarisation vectors in four space-time dimensions. Now we construct explicitly a basis of the space of polarisation states defined by (2.4) for the external massless vector bosons. First, we introduce a polarisation basis vector ε X , valid for both intial-state gluons, which can be written in terms of the linearly independent momenta of the process where the Lorentz invariant coefficients c X i are determined by the system of equations Note that the conditions above constitute a gauge choice in which the reference momentum of either incoming gluon is set to be the momentum of the other gluon. A polarisation vector ε T for both outgoing photons can be constructed analogously: A third basis vector ε Y , pointing out of the scattering plane, is needed to span the space of all possible polarisation vectors for this process: In four dimensions, such a vector can be constructed using the Levi-Civita tensor: Since we consider only QCD corrections to a QED process, neither γ 5 nor Levi-Civita tensors are introduced by the relevant Feynman rules. Consequently, a completely Ddimensional tensor decomposition of this scattering amplitude can be expressed solely in terms of metric tensors and external momenta. Therefore, a contraction of the tensor amplitude with an odd number of ε Y evaluates to zero. A product of two Levi-Civita tensors, however, can be rewritten in terms of metric tensors using which has a straightforward D-dimensional continuation. For a detailed discussion of the subtleties related to the manipulation of Levi-Civita tensors in the construction of projectors for more general cases we refer to Ref. [33].
Applied to the scattering process (2.1), this construction leads to eight projectors where the square bracket [·, ·] means either entry and where only the combinations containing an even number of ε Y are considered. Let us emphasize again that, in order to avoid possible ambiguities in the application of these projectors, all pairs of Levi-Civita tensors are replaced according to the contraction rule (2.10) before being used for the projection of the amplitude. Then the aforementioned projectors are expressed solely in terms of external momenta and metric tensors whose open Lorentz indices are all set to be D-dimensional.
The usual helicity amplitudes can be constructed as circular polarisation states from the linear ones using the relations Analytic results for the LO amplitudes of (2.1) were obtained quite some time ago in Refs. [35][36][37] for massless quark loop contributions and in Refs. [38,39] with massive quark loop contributions. With the linear polarisation projectors defined in (2.11), we re-computed these LO amplitudes analytically, with both massless and massive quark loops. These expressions were implemented in our computational setup for the NLO QCD corrections to the considered process, which we describe below.

UV renormalisation
The bare scattering amplitudes of the process (2.1), denoted byM, beyond LO contain poles in the dimensional regulator ≡ (4−D)/2 arising from ultraviolet (UV) as well as soft and collinear (IR) regions of the loop momenta. In our computation, we renormalise these UV divergences using the MS scheme, except for the top quark mass which is renormalised on-shell.
The bare virtual amplitudeM is a function of the bare QCD couplingα s and the bare top quark massm t . The UV renormalisation ofM is achieved by the replacement and by renormalising the gluon wave function. Here, S = (4π) e − γ E , with γ E the Euler constant. The strong coupling is given by α s = g 2 s /(4π) andμ is an auxiliary massdimensionful parameter introduced in dimensional regularisation to keep the coupling constants dimensionless. The usual renormalisation scale is denoted µ R , and we will useμ = µ R in the following. Both the bare virtual amplitudes and the UV renormalisation constants are expanded in a s ≡ α s (µ R )/(4π). We may write the renormalisation constants as Under the MS scheme for α s with n f massless quark flavours and top-quark loops renormalised on-shell, the renormalisation constants needed in our computation read We write the scattering amplitude for the process gg → γγ, up to second order in a s , in the following formM Here, M B,ren (m t ) and M V,ren (m t ) are the one-loop and UV renormalised two-loop amplitudes, respectively, with the Born kinematics given in (2.1). The mass counterterm amplitudeM CT (m t ) is obtained by inserting a mass counter-term into the heavy quark propagators where a, b, c are colour indices in the fundamental representation. The mass counterterm can also be obtained by taking the derivative of the one-loop amplitude with respect tom t .

Definition of the IR-subtracted virtual part
The UV renormalised virtual amplitude M V,ren still contains divergences arising from soft and collinear configurations of the loop momenta, which appear as poles in the dimensional regulator. We employ the FKS subtraction approach [40] to deal with the intermediate IR divergences, as implemented in the POWHEG-BOX-V2 framework [41][42][43].
For the process gg → γγ, the corresponding integrated subtraction operator is given by To second order in a s the UV renormalised and IR subtracted virtual amplitude is given by Note that the LO amplitude M B,ren needs to be computed to O( 2 ) as it is multiplied by coefficients containing 1/ 2 poles. In practice, we need to supply only the finite part of the born-virtual interference, under a specific definition [43] in order to combine it with the FKS-subtracted real radiation generated within the GoSam/POWHEG-BOX-V2 framework. Explicitly, we compute The renormalisation scale dependence of V fin can be derived from the above definitions, it is given by 23) where µ 0 stands for an arbitrarily chosen initial renormalisation scale.

Evaluation of the virtual amplitude
For the two-loop QCD diagrams contributing to our scattering process there is a complete separation of quark flavours due to the colour algebra and Furry's theorem. Consequently we have n f +1 sets of two-loop diagrams which can be treated separated from each other. The two-loop amplitude has been obtained with the multi-loop extension of the program GoSam [44] where Reduze 2 [45] is employed for the reduction to master integrals. In particular, each of the linearly polarised amplitudes projected out using (2.11) is eventually expressed as a linear combination of 39 massless integrals and 171 integrals that depend on the top quark mass, distributed into three integral families. All massless two-loop master integrals involved are known analytically [17,37,46], and we have implemented the analytic expressions into our code. Regarding the twoloop massive integrals which are not yet fully known analytically, we first rotate to an integral basis consisting partly of quasi-finite loop integrals [47]. Our integral basis is chosen such that the second Symanzik polynomial, F, appearing in the Feynman representation of each of the integrals is raised to a power, n, where |n| ≤ 1 in the limit → 0. This choice improves the numerical stability of our calculation near to the tt threshold, where the F polynomial can vanish. The integrals are then evaluated numerically using pySecDec [48,49]. Examples of contributing two-loop Feynman diagrams are shown in Figure 1. The phase-space integration of V fin is achieved by reweighting unweighted Born events. The accuracy goal imposed on the numerical evaluation of the virtual two-loop amplitudes in the linear polarisation basis in pySecDec is 1 per-mille on both the relative and the absolute error. We have collected 6898 phase space points out of which 862 points fall into the diphoton invariant mass window m γγ ∈ [330, 360] GeV. We have also calculated a further 2578 phase space points restricted to the threshold region.

Computation of the real radiation contributions
The real radiation matrix elements are calculated using the interface [50] between GoSam [51,52] and the POWHEG-BOX-V2 [41][42][43], modified accordingly to compute the real radiation corrections to loop-induced Born amplitudes. Only real radiation contributions which contain a closed quark loop at the amplitude level are included. We also include the qq initiated diagrams which contain a closed quark loop, even though their contribution is numerically very small. Examples of Feynman diagrams contributing to the real radiation amplitude are shown in Figure 2. The diagrams in which one of the photons is radiated off a closed fermion loop and the other photon is radiated off an external quark line vanish due to Furry's theorem.

Treatment of the threshold region
When the partonic centre-of-mass energy is close to the threshold for the production of a tt pair, the top quarks are produced with a non-relativistic velocity such that Coulomb interactions between the top quarks can play a significant role. In the case of the top-loop induced contribution to diphoton production, the Coulomb singularity appears in the form of a logarithmic dependence on the velocity first at two-loop order, due to the exchange of a soft gluon between the top quarks in the loop. To overcome this issue and correctly describe the threshold, we employ the so-called non-relativistic QCD (NRQCD) [53][54][55][56], which is an effective field theory designed to describe nonrelativistic heavy quark-antiquark systems in the threshold region.

NRQCD amplitude
To the order which we consider here, the amplitude can be expressed as a coherent sum of light quark loop contributions and the top quark loop contributions, where α e = e 2 /(4π) and Q q denotes the electric charge of quark q. In our computation, the NRQCD expansion of the amplitude M t near the tt threshold is performed according to the formalism explained in more detail in Refs. [14,57]. Near the production threshold of an intermediate tt pair, m γγ 2m t , we define 2) and the scattering angle is given by Close to threshold, the amplitude M t can be parametrised as [14,57] where E = E + iΓ t includes the top-quark decay width Γ t 2 . Note that the P-wave contribution B t,P (β) G P ( 0; E) starts at O(β 3 ). In this parametrisation, the amplitude M NR t is split into two parts: B t (β) G( 0; E), which contains the tt bound state effects, and A t (θ), which does not. The term B t (β) G( 0; E) contains the effects from resumming the non-relativistic static potential interactions, where the Green's function G( 0; E) is obtained by solving the non-relativistic Schrödinger equation describing a colour-singlet tt bound state: with the QCD static potential [59,60] The mass m t appearing in (3.5) is the pole mass of the top quark. G( 0; E) is the r → 0 limit of the Green's function G( r; E). The real part of the NLO Green's function at r = 0 is divergent and therefore has to be renormalised. We adopt the MS scheme, thus introducing a scale µ into the renormalised Green's function [61][62][63][64]. The coefficient B t (β) can be obtained from the Wilson coefficients of the ggtt and γγtt operators [14] in the NRQCD effective Lagrangian for the process gg → γγ. The term A t (θ) encompasses the non-resonant corrections, resulting from quark loops with large virtuality which can be systematically computed order by order in α s . Both A t and B t can be expanded perturbatively in α s . For the process gg → γγ, corrections to B t have been calculated up to O(α s ) and O(β 2 ) in Ref. [14], where explicit expressions of B t at the leading order for all relevant helicity configurations can be found. Here we repeat for completeness the expressions for the S-wave tt resonance we are considering. For the S-wave the B t coefficients are independent of the scattering angle. We use the notation G(β) ≡ G( 0; E) and Note that an overall factor of α s already has been extracted from the amplitude (see Eq.
The expansion of the Green's function in α s is given by where [64,68] For A t (θ), we can make use of a partial-wave decomposition in terms of Wigner functions d J hh (θ), (3.14) where h = −λ 1 + λ 2 and h = λ 3 − λ 4 .

Matched amplitude
We would like to retain NRQCD resummation effects and, at the same time, keep the cross section accurate up to NLO in the fixed-order power counting. We define the "NRQCD-matched" amplitude as [14] M match where the first term is the fixed-order amplitude, the second term describes the threshold according to NRQCD and the third term M OC ≡ B t G( 0; E) subtracts double counted contributions included in both the fixed-order amplitude and NRQCD contribution. The M OC term in a fixed-order computation should be expanded to the same order as the fixed-order amplitude. Expanding (3.15) to next-to-leading order, we have (3.17) Inserting into the matched amplitude we obtain, The NLO-matched cross section is obtained by squaring the matched amplitude and adding the corresponding real-radiation. Upon squaring the matched amplitude we obtain,

Matched cross section
We define our NLO-matched cross section as follows 21) where N ij contains the flux factor and the average over spins and colours of the initial state partons of flavour i and j, e.g. N gg = 1 2s 1 64 1 4 . And we have introduced the luminosity factors L ij , defined as is the parton distribution function (PDF) of a parton with momentum fraction x and flavour i (including gluons) and µ F is the factorisation scale. The 2-and 3-particle phase-space integration measures are denoted by dΦ 2 and dΦ 3 . The symbol c ≡ 32π α e Q 2 t T R δ a 1 a 2 collects constants which have been extracted in the definition of M t . The real-radiation contributions with the factors of a s extracted are symbolically denoted by M R, [ij] and the collinear-subtraction counterterm is denoted by σ C . We do not include resummation effects in the real-radiation because it is suppressed by a factor of β. The M OC to respectively O(β 3 ) and O(β 2 ) using the expressions stated in Section 3.1. At the two-loop order, the UV-renormalised and IR-subtracted fixed-order amplitude M t has a Coulomb singularity which is logarithmically divergent in the limit β → 0. This singularity is, however, subtracted by the expanded term M OC , while a resummed description of the Coulomb interactions is added back by the term B t G( 0; E). For this purpose, we evaluate the Schrödinger equation (3.5) numerically [69] to obtain G( 0; E), where we include O(α s ) corrections to the QCD potential [59,60]. Unlike the calculation in [14], we also include O(α s ) corrections to B t as listed above.

Results
Our numerical results are calculated at a hadronic centre-of-mass energy of 13 TeV, using the parton distribution functions PDF4LHC15 nlo 100 [70][71][72][73] interfaced via LHAPDF [74], along with the corresponding value for α s . For the electromagnetic coupling, we use α = 1/137.035999139. The mass of the top quark is fixed to m t = 173 GeV. The top-quark width is set to zero in the fixed order calculation, and to Γ t = 1.498 GeV in the numerical evaluation of the Green's function G( 0; E, µ) in accordance with Ref. [14]. We use the cuts p min T,γ 1 = 40 GeV, p min T,γ 2 = 25 GeV and |η γ | ≤ 2.5. No photon isolation cuts are applied. The factorisation and renormalisation scale uncertainties are estimated by varying the scales µ F and µ R . Unless specified otherwise, the scale variation bands represent the envelopes of a 7-point scale variation with µ R,F = c R,F m γγ /2, where c R , c F ∈ {2, 1, 0.5} and where the extreme variations (c R , c F ) = (2, 0.5) and (c R , c F ) = (0.5, 2) have been omitted. The dependence on the scale µ introduced by renormalisation of the Green's function G( r; E) in our NRQCD matched results is investigated separately.

Fixed-order calculation
We have validated the massless NLO cross section by comparison to MCFM version 9.0 [75] and find agreement within the numerical uncertainties for all scale choices. We also compared against the results shown in [20] and find agreement for the central scale choice, however we find a smaller scale uncertainty band. We remark that the helicity amplitudes can also be computed via first performing the Lorentz tensor decomposition, using the form factor projectors given in Ref. [37], and then evaluating contractions between the corresponding Lorentz structures and external polarisation vectors in 4 dimensions using the spinor-helicity representations. This amounts to obtaining helicity amplitudes defined in the t'Hooft-Veltman scheme [76]. We confirm numerically that the same finite remainders are obtained for all helicity configurations at a few chosen test points (while the unsubtracted helicity amplitudes do differ starting from the subleading power in ). As a further cross check, we evaluate our amplitude with t ↔ u interchanged and confirm that the helicity amplitudes are permuted as expected.

NRQCD amplitude
Numerical values for the coefficients A J t,{λ i } at leading-order in α s up to J = 4 are given in Ref. [14]. We have used them as a check of our numerical calculation of the Born amplitude. We also evaluated the massive two-loop amplitude at 615 phase space points with m t = 173 GeV in the ranges 0 < cos (θ) < 1 and 0.001 ≤ β ≤ 0.2, using the program pySecDec [48,49]. The amplitude can numerically be fitted to a suitable ansatz in β and cos θ. We have compared the coefficients of terms proportional to ln (β) to the known analytical results based on expanding equation (3.8) and find good agreement. Note that the coefficients of terms not proportional to ln (β) receive contributions from the unknown term A (1) (θ) and can therefore not be checked this way.

Invariant mass distribution of the diphoton system
The distribution of the invariant mass of the photon pair is shown in Fig. 3 for invariant masses up to 1 TeV, where we show purely fixed order results at LO, at NLO with five massless flavours and at NLO including massive top quark loops. The ratio plots show the K-factor including the full quark loop content and the ratio between the full and the five-flavour NLO cross-section. We observe that the scale uncertainties are reduced at NLO, and that the top quark loops enhance the differential cross section for m γγ values far beyond the top-quark pair-production threshold, asymptotically approaching the n f = 6 value [31]. In Fig. 4 we zoom into the threshold region, still showing fixed order results only. We can clearly see that after the top quark pair production threshold, the full result shows a dip and then changes slope, which is due to the fact that the two-loop amplitude contains the exchange of a Coulomb gluon (see top left diagram of Fig. 1), as explained in Section 3. In Ref. [14] it was suggested that this characteristic "dip-bump structure" could be used for a determination of the top quark mass which is free from top quark reconstruction uncertainties, at least at the FCC where the statistical uncertainties for this process would be very small, and the systematic uncertainty due to the finite resolution of the photon energies and angles should be at least as good as at the LHC, where it is at the sub-percent level [77,78]. In Fig. 5 we show the m γγ -distribution in the threshold region which results from a combination of the fixed-order NLO (QCD) calculation with the resummation of Coulomb gluon exchanges as described in Section 3.2. The scale band in this figure are produced by varying only µ, the scale associated to the renormalisation of the Green's function. We observe that the dependence on the scale µ is considerably reduced at NLO compared to the leading-order matched cross-section. The scale band at NLO is comparable to the size of our numerical uncertainties. Further, our leading-order matched cross-section shows a milder dependence on µ than the one presented in [14]. This is due to the inclusion of NLO-terms in the coefficient B t (β), which have been omitted in [14]. We do not consider the effects from a colour-octet tt state because the corresponding Green's function is monotonically increasing in the resonance region [67] and therefore not expected to move the position of the dip significantly. Now let us address the prospects to measure the top quark mass from the threshold behaviour of the m γγ distribution. In Ref. [14] it was argued that the characteristic dipbump structure does not change its location in the m γγ spectrum under scale variations, only the overall normalisation is changing. It was also anticipated that the inclusion of the fixed order two-loop amplitude would reduce this uncertainty. Indeed we find that the NLO corrections reduce the scale uncertainties due to 7-point µ R , µ F -variations from about 20% at LO to just below the 10% level at NLO. However, the pronounced dip-bump structure present in the LO resummed calculation is partly washed out in the NLO NRQCD-improved calculation.

Conclusions and Outlook
We have calculated the production of a photon pair in gluon fusion at order α 3 s , including massive top quark loops. This calculation, which is NLO for the gluon initiated channel, is formally part of the N 3 LO corrections to the pp → γγ process. However, the gluon channel is important at the LHC due to the large gluon luminosity. The top quark loops have a considerable impact on the diphoton invariant mass spectrum, at values of m γγ larger than about 800 GeV they enhance the m γγ differential cross section by more than 50%. The region around the top quark pair production threshold in the diphoton invariant mass spectrum is particularly interesting. The fixed order amplitude has a divergence starting at two loops due to Coulomb gluon exchange. We have used NRQCD methods to resum the bound state effects in order to obtain a more reliable description of the threshold region. Matching the resummed calculation to our fixed order NLO calculation we observe a reduction of the renormalisation and factorisation scale uncertainties in the threshold region by more than a factor of two, and an even more drastic reduction of the scale uncertainty related to the renormalised NLO Green's function. These results are promising in view of the possibility of measuring the top quark mass from the characteristic behaviour of the diphoton invariant mass spectrum around the top quark pair production threshold. In Ref. [14], it was found that the LO resummed result shows a characteristic "dip-bump" structure and the conclusion was that this would allow a precise measurement of the top quark mass with the statistics and photon resolution projected for the FCC, once an NLO calculation is available such that the scale uncertainties are reduced. Now we indeed found that at NLO, the scale uncertainties are reduced, on the other hand the characteristic "dip-bump" structure is partly washed out. Nonetheless the change in slope is clearly visible. A detailed assessment of whether this behaviour is still pronounced enough for a top quark mass measurement once all channels contributing to this observable are included deserves further study. It also requires a detailed study of the prospective experimental uncertainties. Furthermore, it would be interesting to investigate other top quark mass schemes, however this is beyond the scope of this paper and therefore we defer it to future work.