At the LHC, a diphoton mass spectrum \(\mathrm{d}\sigma /\mathrm{d}m_{\gamma \gamma }\) has attracted broad attention for observations of the properties of the Higgs boson in the standard model (SM) [1,2,3,4,5] and searches for new phenomena beyond the SM [6,7,8,9,10,11]. At hadron colliders, pairs of high-\(p_T\) photon are produced by \(q\bar{q}\) annihilation and gluon-fusion mechanisms [12,13,14], and processes which involve fragmentation-photon contributions [15]. The \(gg\rightarrow \gamma \gamma \) process, which is the main focus of this letter, is described by loop diagrams with quarks in the SM. The analytic expression of the one-loop amplitude has been known for a long time for both the massless- and the massive-quark loops [16,17,18,19,20,21,22]. The two-loop amplitude has been calculated only for the massless-quark loops [23,24,25].

The threshold structure of the massive-quark-loop amplitude deserves particular interests [21, 26, 27] where the massive quark is regarded as the top quark or a hypothetical particle beyond the SM. Beyond the one-loop level, the amplitude receives large QCD corrections due to the Coulomb-gluon exchanges between the nearly on-shell and low-velocity heavy quarks in s-channel. Thus, the description of the amplitude requires an elaborate treatment based on the non-relativistic QCD formalism. For \(gg\rightarrow \gamma \gamma \) process, such a study cannot be found in the literature. The aim of this letter is to compile present knowledge of the non-relativistic QCD theory for the description of the bound-state effects in the massive-quark-loop amplitude, and to present a dedicated and quantitative study on the diphoton mass spectrum near the \(t\bar{t}\) threshold. Our framework follows the preceding studies on \(h\rightarrow \gamma \gamma \) [28, 29], and some of our numerical results overlap with that in Ref. [26].

We discuss further to utilize the predicted mass spectrum for a precise determination of the top-quark mass, which is one of the fundamental parameters in the SM. Although the top-quark mass has been measured with an error of sub-GeV level [30], its interpretation in terms of well-defined mass parameters is not settled yet in perturbative QCD. It is known that the well-defined mass parameters can be determined by using the threshold scan method at future \(e^+e^-\) colliders [31,32,33]. We show that the diphoton mass spectrum measurement can be a considerable alternative to it at hadron colliders. The application of the formula for physics beyond the SM will be reported elsewhere.

Table 1 Numerical coefficients of \({\mathcal A}^{J}_{t,\lambda _1\lambda _2\lambda _3\lambda _4}\)
Full size table

We start by introducing the scattering amplitude for \(gg\rightarrow \gamma \gamma \) at the one-loop level with the top quark, and provide an easy-to-use expression for its threshold behavior. By using the all-outgoing convention for the momenta (\(p_i\)) and helicities (\(\lambda _i\)), \(g^{a_1}(-p_1,-\lambda _1) + g^{a_2}(-p_2,-\lambda _2) \rightarrow \gamma (p_3,\lambda _3) + \gamma (p_4,\lambda _4)\), the one-loop amplitude is written as

$$\begin{aligned}&{\mathcal M}_{gg\rightarrow \gamma \gamma } (\{p_i\};\{\lambda _i\};a_1,a_2) = 4\alpha \alpha _s \delta ^{a_1a_2} \nonumber \\&\quad \times \left[ \left( \sum _{j=1}^{n_f}Q_j^2\right) M_{q,\{\lambda _i\}}(\{p_i\}) + Q_t^2 M_{t,\{\lambda _i\}}(\{p_i\};m_t) \right] , \end{aligned}$$
(1)

where \(M_q\) is the contribution from the massless-quark loop with five flavors (\(n_f=5\)), and \(M_t\) from the top-quark loop with the top-quark pole-mass, \(m_t\). The amplitude for the top-quark loop near the threshold is expressed as

$$\begin{aligned}&M_{t,\{\lambda _i\}} = {\mathcal A}_{t,\{\lambda _i\}}(\theta ) + {\mathcal B}_{t,\{\lambda _i\}} G^{(0)}(\mathbf {0};E) + {\mathcal O}(v^2), \end{aligned}$$
(2)

where \(E \equiv m_{\gamma \gamma }-2m_t\simeq m_tv^2\) and \(v = \sqrt{1-4m_t^2/m^2_{\gamma \gamma }}\). \(G^{(0)}(\mathbf {0};E) \equiv -m_t^2/(4\pi )\sqrt{-E/m_t-i\epsilon }\) is the \(t\bar{t}\) Green function in S-wave without QCD effects. The first term which is energy-independent, represents the contribution from the hard-momentum integral. The second term which is \({\mathcal O}(v)\), represents the contribution from the soft-momentum loop where the top-quarks can be on-shell. At the one-loop level, all the imaginary part of the amplitude originates from \(G^{(0)}\) above the threshold, \(E\ge 0\). \({\mathcal A}_t\) depends on the scattering angle \(\theta \), while \({\mathcal B}_t\) is independent of \(\theta \) because only the spin-singlet \(t\bar{t}\) state contributes at this order. For \(\{\lambda _i\}=\lambda _1\lambda _2\lambda _3\lambda _4\), we find \({\mathcal B}_{t,++++} = -{\mathcal B}_{t,--++} = -4\pi ^2/m_t^2\), while \({\mathcal B}_{t,-+++} = {\mathcal B}_{t,-+-+} = 0\). For the other combinations of the helicity, \({\mathcal B}_t\) as well as \({\mathcal A}_t\) can be written in terms of them. For a description of \({\mathcal A}_t\), we make use of the partial-wave decomposition with numerical coefficients. The \({\mathcal A}_t\) term is expanded as

$$\begin{aligned} {\mathcal A}_{t,\{\lambda _i\}} (\theta ) = \sum _{J=0}^{\infty } (2J+1){\mathcal A}^{J}_{t,\{\lambda _i\}} d^{J}_{\mu \mu '}(\theta ) , \end{aligned}$$
(3)

where \(\mu =-\lambda _1+\lambda _2\) and \(\mu '=\lambda _3-\lambda _4\). Because \({\mathcal B}_t\) is constant, \({\mathcal B}_t\) has only the \(J=0\) component, \({\mathcal B}_t={\mathcal B}_t^{J=0}\). In Table 1, we list the numerical values of \({\mathcal A}^{J}_{t}\) for J up to 4. The Wigner d-functions \(d^{J}_{\mu \mu '}\) can be found in the literature. We find that the expansion up to \(J=4\) gives a sufficiently good approximation.

We incorporate the \(t\bar{t}\) threshold effects into the Green function by evaluating it with the QCD potential [34, 35]. The amplitude with the threshold effects is expressed as [28]

$$\begin{aligned} M_{t,\{\lambda _i\}}^\mathrm{thr} = {\mathcal A}^{J=0}_{t,\{\lambda _i\}} + {\mathcal B}_{t,\{\lambda _i\}} G(\mathbf {0};{\mathcal E}) + {\mathcal A}^{J>0}_{t,\{\lambda _i\}}(\theta ) , \end{aligned}$$
(4)

where we define \({\mathcal A}^{J>0}_t(\theta ) = {\mathcal A}_t(\theta ) - {\mathcal A}_t^{J=0}\) and \({\mathcal E}=E+i\Gamma _t\) with the top-quark decay width, \(\Gamma _t\). The Green function is defined by the following Schrödinger equation:

$$\begin{aligned} \left[ \left\{ -\frac{\nabla ^2}{m_t} + V(r) \right\} - {\mathcal E} \right] G(\mathbf {r};{\mathcal E}) = \delta ^3(\mathbf {r}), \end{aligned}$$
(5)

where V(r) is the QCD potential. For the \(t\bar{t}\) system, we can utilize the perturbatively calculated potential. The real part of the Green function at \(\mathbf {r}=\mathbf {0}\) is known to be divergent, thus has to be renormalized. We adopt the \(\overline{\mathrm{MS}}\) renormalization scheme in dimensional regularization [36,37,38]. An artificial scale \(\mu \) is introduced to the renormalized Green function. By matching with the one-loop amplitude, the amplitude is finally expressed as

$$\begin{aligned} M_{t,\{\lambda _i\}}^\mathrm{match}&= M_{t,\{\lambda _i\}} + {\mathcal B}_{t,\{\lambda _i\}}[ G(\mathbf {0};{\mathcal E}) - G^{(0)}(\mathbf {0};E) ] . \end{aligned}$$
(6)

Before moving to the numerical evaluation, we discuss the order of the corrections in the non-relativistic QCD formalism. Taking v and \(\alpha _s\) as the expansion parameters, the leading-order contribution is the \({\mathcal A}_t^{J=0}\) term which is constant, and the \({\mathcal B}_tG(\mathbf {0};{\mathcal E})\) term is at the next-to-leading order (NLO). There is another NLO term in the two-loop amplitude, which is an \({\mathcal O}(\alpha _s)\) correction to \({\mathcal A}_t\). However, this has not been calculated yet for the massive-quark contribution. Indeed, this term is required for the consistent calculation of the threshold corrections up to NLO in order that the scale dependence of the real part of the Green function is canceled with the \({\mathcal O}(\alpha _s)\) term of \({\mathcal A}_t\) [29]. In our calculation, we do not include the \({\mathcal O}(\alpha _s)\) \({\mathcal A}_t\)-term, thus the scale dependence remains in the threshold amplitude. We treat it as an uncertainty of our calculation.

Since the leading contribution to the squared amplitude is the absolute square of the sum of \(M_q\) and the \({\mathcal A}_t\) term where both are independent of energy, the uncertainty of \({\mathcal A}_t\) term mainly affects the overall normalization of the diphoton mass spectrum. On the other hand, some of the NNLO corrections improve the description of the \(t\bar{t}\) resonances. Therefore, for the sake of a precise and stable prediction of the resonance structure, it is worthwhile to include the available NNLO corrections even though we cannot reach the full NLO accuracy. The known corrections are (1) the NLO correction to the Green function, (2) the \({\mathcal O}(\alpha _s)\) correction to \({\mathcal B}_t\), and (3) the \({\mathcal O}(\alpha _s)\) correction to \(\Gamma _t\). First, the NLO correction to the Green function is incorporated by solving the Schrödinger equation with the NLO QCD potential [39, 40] given by

$$\begin{aligned} V(r) = - C_F \frac{\bar{\alpha }_s}{r} \left[ 1 + \frac{\bar{\alpha }_s}{4\pi } \left\{ 2\beta _0[ \ln {(\mu _B r)} + \gamma _E ] + a_1 \right\} \right] \!, \end{aligned}$$
(7)

where \(\bar{\alpha }_s=\alpha _s(\mu _B)\), \(\beta _0=11/3\cdot C_A-2/3\cdot n_f\) and \(a_1=31/9\cdot C_A-10/9\cdot n_f\) with \(C_F=4/3\) and \(C_A=3\). We will show later that evaluating the Green function beyond LO is crucial for a reliable prediction. Second, the \({\mathcal O}(\alpha _s)\) correction to \({\mathcal B}_t\) can be derived from the \({\mathcal O}(\alpha _s)\) hard-vertex corrections to the \(gg\rightarrow t\bar{t}\) and \(t\bar{t}\rightarrow \gamma \gamma \) processes. The hard-vertex factor to the \(gg\rightarrow t\bar{t}\) cross section in the color-singlet channel reads \(1+(\alpha _s/\pi )h_1\) with [41,42,43]

$$\begin{aligned} h_1 = C_F \left( - 5 + \frac{\pi ^2}{4} \right) + C_A \left( 1 + \frac{\pi ^2}{12} \right) + \beta _0\ln {\left( \frac{\mu _R}{2m_t}\right) } , \end{aligned}$$
(8)

where \(\mu _R\) is the renormalization scale of \(\alpha _s\). The corresponding factor for \(t\bar{t}\rightarrow \gamma \gamma \) reads only the first term of Eq. (8). By using them, \({\mathcal B}_t\) with the \({\mathcal O}(\alpha _s)\) correction is given as \({\mathcal B}_t = {\mathcal B}^{(0)}_t[ 1 + (\alpha _s/\pi ) b_1 ]\) with

$$\begin{aligned} b_1 = C_F \left( - 5 + \frac{\pi ^2}{4} \right) + \frac{C_A}{2} \left( 1 + \frac{\pi ^2}{12} \right) + \frac{\beta _0}{2} \ln {\left( \frac{\mu _R}{2m_t}\right) } . \end{aligned}$$
(9)

Finally, the \({\mathcal O}(\alpha _s)\) correction to \(\Gamma _t\) has been calculated in Refs. [44,45,46]. However, we treat \(\Gamma _t\) as an input parameter in our study. Identification and derivation of the remaining NNLO corrections are beyond the scope of this letter.

Fig. 1
figure 1

\(gg\rightarrow \gamma \gamma \) contribution to the diphoton invariant-mass distribution near the \(t\bar{t}\) threshold at the LHC 13 TeV. Top panel is for the LO Green function, and bottom panel is for the NLO Green function

Full size image

We present numerical studies for the \(gg\rightarrow \gamma \gamma \) amplitude as well as the cross sections at the LHC. In Fig. 1, we plot the \(gg\rightarrow \gamma \gamma \) contribution to the hadronic differential cross section, \(\mathrm{d}\sigma /\mathrm{d}m_{\gamma \gamma }\), for the LHC 13 TeV with kinematical cuts of \(|\eta _\gamma |<2.5\) and \(p_{T}^\gamma >40\) GeV [8]. Both the massive- and the massless-quark loops are included. We use the CT14NLO gluon distribution function [47], and take the renormalization and factorization scales as \(\mu _R=\mu _F=m_{\gamma \gamma }\). The Green function is evaluated by numerically solving Eq. (5) with the LO or NLO QCD potential following the method described in Ref. [48]. The scale of \(\alpha _s\) in the QCD potential is taken as the same as the renormalization scale \(\mu \) of the Green function, which we vary from 20 to 160 GeV. The result with the one-loop amplitude is also plotted for comparison. In the plots, we observe that the distributions show a characteristic structure near \(m_{\gamma \gamma }\simeq 2m_t=346\) GeV; it shows a dip and then a small bump below the threshold [26]. We find that, if we employ the LO Green function, the shape of the distribution changes by the scale choice. In contrast, by using the NLO Green function the shape of the distribution is quite stable apart from the overall normalization. The positions of the dip and the bump are shifted by the choice of \(\mu \) by around 0.6 GeV. A relatively large uncertainty appears as the overall size of the cross section, which amounts to about 10%. This uncertainty originates mainly from the lack of the \(\mathcal{O}(\alpha _s)\) correction in the \({\mathcal A}_t\) term. We note that there exists another source of the uncertainty for the overall normalization, which is the scale choice of \(\mu _R\) and \(\mu _F\). For the LHC 13 TeV, changing these scales from \(m_{\gamma \gamma }/2\) to \(2m_{\gamma \gamma }\) varies the cross section by about 20%.

Fig. 2
figure 2

\(gg\rightarrow \gamma \gamma \) amplitudes for \(\{\lambda _i\}=++++\) in the \(J=0\) channel for \(m_{\gamma \gamma }=330\)–360 GeV with points in a 5-GeV step. For the illustrative convenience, we set \(m_t=172.5\) GeV in this plot

Full size image

For a better understanding of the behavior of the cross section, we plot in Fig. 2 the \(gg\rightarrow \gamma \gamma \) amplitudes in a complex plane for \(\{\lambda _i\}=++++\) in the \(J=0\) channel, by varying \(m_{\gamma \gamma }\) from 330 to 360 GeV. The massless-quark-loop amplitude gives a constant contribution, \(11/9\cdot M_{q,++++}^{J=0}=11/9\). The total amplitude \(M^{J=0}_{++++} = 11/9\cdot M^{J=0}_{q,++++} + 4/9\cdot M^{J=0}_{t,++++}\) with the one-loop-level \(M^{J=0}_{t,++++}\) is drawn in the black line. Below the threshold, the two amplitudes, \(M_q\) and \(M_t\), are pure real, and their relative sign is negative. Therefore, there is a destructive interference, and the total amplitude goes toward the origin by increasing \(m_{\gamma \gamma }\) until the threshold. Above the threshold, the amplitude gains an imaginary part and the real part tends to increase along with \(m_{\gamma \gamma }\). At the high-energy limit, where the top quark can be assumed to be massless, the imaginary part goes to zero and the total amplitude arrives at \(M^{J=0}_{++++}=15/9\).Footnote 1 The amplitude with the threshold corrections calculated with the NLO (LO) Green functions are plotted in colored solid (dotted) lines for \(\mu =40\), 80 and 160 GeV. The imaginary part of the amplitude is non-zero even below the threshold, which comes from the finiteness of \(\Gamma _t\). The size of the imaginary part increases rapidly above \(m_{\gamma \gamma }=340\) GeV with showing a resonance-like curve just below the threshold. The scale dependence of the Green function originates from the two sources, one in the QCD potential and the other from the real-part renormalization. For the NLO Green function, the former is well suppressed and the latter affects only the real part of the amplitude by a constant for any \(m_{\gamma \gamma }\). For the LO Green function, both effects are large and the amplitude shows a complicated scale dependence. Especially, there remains a scale dependence in the imaginary part of the amplitude. This explains the reason that the shape of the invariant-mass distribution is stable by using the NLO Green function in contrast to the LO Green function. Although the uncertainty in the real part of the amplitude is significant, it leads only the 10% level uncertainty to the cross section, due to the presence of the large imaginary part and the light-quark-loop contribution.

Fig. 3
figure 3

Scale dependence of the dip and bump positions in the diphoton mass spectrum at the LHC evaluated with the NLO Green function, and the energy level of the 1S toponium evaluated in the pole-mass scheme at NLO as well as those in the \(\overline{\mathrm{MS}}\)-mass scheme up to N\(^{3}\)LO. \(m_t^\mathrm{pole}=173\) GeV or \(\overline{m}_t = 163\) GeV is used

Full size image

In Fig. 3, we show the scale dependence of the dip and bump positions, \(M_\mathrm{dip}\) and \(M_\mathrm{bump}\), respectively, in the diphoton mass spectrum at the LHC evaluated with the NLO Green function. In addition, we plot the 1S energy level of the \(t\bar{t}\) bound-state (toponium) at NLO [\(\mathcal{O}(\alpha _s^3m_t)\)], \(M^{(1)}_{1S}\), which is in good approximation the resonance peak position in the NLO Green function. We find the scale variation of the Green function affects the difference of the two mass scales, \(M_\mathrm{dip}\) and \(M_{1S}\) (and also, \(M_\mathrm{bump}\) and \(M_{1S}\)), by only around 20 MeV (40 MeV). This indicates that the connection of the dip (bump) position and the 1S resonance mass is sufficiently solid under uncertainties of the Green function. The toponium energy levels have been calculated up to \(\mathcal{O}(\alpha ^5_sm_t)\) [53, 54] in non-relativistic QCD, and it is well known that the prediction becomes significantly accurate when it is expressed in terms of the short-distance mass to cancel the renormalon ambiguity. By using the \(\mathcal{O}(\alpha ^5_sm_t)\) formula for the spin-singlet case [53, 54] and the \(\overline{\mathrm{MS}}\) mass with \(\overline{m}_t=163\) GeV, we also plot the 1S energy level, \(\overline{M}^{(n)}_{1S}\), at N\(^{n}\)LO up to \(n=3\) in Fig. 3. It can be seen that the convergency is good, and the scale uncertainty is reduced to around 100 MeV or below.Footnote 2 By combining these arguments, the dip and bump positions can be accurately predicted by including higher-order corrections with the short-distance mass. More detailed studies will be presented in a future publication.

We propose to use the diphoton mass spectrum near the \(t\bar{t}\) threshold for a precise determination of the top-quark mass in hadron-collider experiments. Figure 4 shows the diphoton mass spectra via \(gg\rightarrow \gamma \gamma \) with different values of \(m_t\) for the LHC at \(\sqrt{s}=13\) TeV (top panel) and the proposed future circular collider (FCC) at \(\sqrt{s}=100\) TeV [56,57,58] (bottom panel). We utilize the NLO Green function with \(\mu =40\) GeV. \(\Gamma _t=1.498\) GeV is fixed for any \(m_t\). The setup for the gluon distribution function and acceptance cuts is the same as that for Fig. 1. An additional cut \(p_T^\gamma >0.4m_{\gamma \gamma }\) is applied for the FCC case which enhances the selection efficiency of the \(J=0\) partial-wave contribution. One can clearly see in Fig. 4 that the bump position shifts in proportion to \(m_t\). Consequently, we can extract \(m_t\) from the diphoton mass spectrum. Since a photon is a clean object and not directly affected by final-state QCD interactions, this measurement would be quite transparent experimentally and theoretically. Especially, systematic errors of photon momentum reconstruction are much smaller than those of jet momentum which are the major source of the systematic error in the current \(m_t\) measurement. These virtues are shared with leptonic-observable methods proposed in Refs. [59,60,61].

Fig. 4
figure 4

\(gg\rightarrow \gamma \gamma \) contribution to the diphoton invariant-mass distribution for different \(m_t\). Top panel is for the LHC 13 TeV and bottom panel is for the FCC 100 TeV

Full size image

In order to estimate the sensitivity of the method, we perform pseudo-experiments assuming the LHC 13 TeV with 3 ab\(^{-1}\) data, and the FCC 100 TeV with 1 ab\(^{-1}\) and 10 ab\(^{-1}\) data. We prepare event samples of the signal \(gg\rightarrow \gamma \gamma \) events and the background events by other sources for the range \(m_{\gamma \gamma }=[300,400]\) GeV with applying the above acceptance cuts. The signal events are generated based on the predicted distribution assuming \(m_t^\mathrm{true}=173\) GeV. The background events are generated by Diphox [15] at LO with \(q\bar{q}\rightarrow \gamma \gamma \), one-direct–one-fragmentation, and two-fragmentation contributions. The total number of events is fixed by using the observed data to take into account detector efficiency and a K-factor from higher-order corrections. We read off a corresponding correction factor of \(C\simeq 1.2\) from the LHC 13 TeV diphoton analysis by the ATLAS Collaboration [10]. For simplicity we apply the same C for the FCC case. The signal-to-background ratio, which is crucial to the mass sensitivity, is subject to theoretical uncertainties of the cross-section calculations, such as the choice of scales \(\mu _R\), \(\mu _F\) and \(\mu \), non-calculated higher-order corrections, and also a definition of isolated photons [62]. Based on the LO calculations for both the signal and the background processes, the ratio is estimated to be 10% at the LHC 13 TeV and 30% at the FCC 100 TeV. On the other hand, the ratio is estimated to be 5% at the LHC where the QCD NNLO corrections are included in the background calculation [12, 14], while 10% at the FCC where the QCD NLO corrections are included in the background calculation [58]. We note that a recent study in Ref. [14] indicates that the ratios become closer to the LO estimates when the NLO corrections are included additionally to the signal process. Considering these estimations, we take the ratio to be 5–10% at the LHC, while 10–30% at the FCC in this study.

The sample \(m_{\gamma \gamma }\) distributions are fitted with the sum of the signal prediction which depends on \(m_t\) plus an analytic smooth function for the background, taking the signal-to-background ratio as a fitted parameter. The background function is taken as \((1-x^{1/3})^a\) where \(x=m_{\gamma \gamma }/\sqrt{s}\) and a is a parameter to be fitted. Notice that our fitting procedure does not rely on the value of the signal-to-background ratio nor the accurate prediction of the background shape. We perform least-squares fits to the binned \(m_{\gamma \gamma }\) distribution in the interval [300, 400] GeV with the bin width of 1 GeV. By repeating the pseudo-experiment, we obtain the expected statistical error \(\Delta m_{t}\) from the distribution of the fitted \(m_t\). For the LHC 3 ab\(^{-1}\), the obtained \(m_t\) distribution is not Gaussian, while it has a peak at \(m_t=m_t^\mathrm{true}\). We approximate the distribution as Gaussian and obtain \(\Delta m_t\simeq 2\) to 3 GeV for the signal ratio 10–5%. For the FCC, by assuming the signal-to-background ratio to be 30%, the distribution behaves as Gaussian and we obtain \(\Delta m_t=0.2\) GeV (0.06 GeV) for 1 ab\(^{-1}\) (10 ab\(^{-1}\)). When the ratio is assumed to be 10%, we obtain \(\Delta m_t=0.6\) GeV (0.2 GeV) for 1 ab\(^{-1}\) (10 ab\(^{-1}\)). We find that the correlation between two fitted parameters, \(m_t\) and the signal-to-background ratio, is weak.

Before closing, we present several comments. The systematic error of photon energy scale is about 0.5% [63] in the ATLAS detector and about 0.3% [64] in the CMS detector. Thus we naively expect the systematic error of \(\delta m_t^\mathrm{sys.}\lesssim 1\) GeV at the future LHC measurement. For more realistic estimation at the LHC as well as at the FCC, simulation studies with detailed detector performance are required. Beyond the one-loop level, the mass renormalization scheme becomes explicit. With the signal distribution expressed in terms of theoretically well-defined masses, the top-quark \(\overline{\mathrm{MS}}\) mass can be extracted directly from the diphoton mass spectrum. Measuring the short-distance mass from the resonance structure is conceptually equivalent with the threshold scan method in \(e^+e^-\rightarrow t\bar{t}\). In the \(e^+e^-\) case, the threshold production cross section is established up to N\(^3\)LO in non-relativistic QCD [65, 66]. In the diphoton case at hadron colliders, only the one-loop \(gg\rightarrow \gamma \gamma \) amplitude has been known, and thus the NLO calculation has not been completed yet. To complete, one requires the two-loop \(gg\rightarrow \gamma \gamma \), one-loop \(gg\rightarrow \gamma \gamma g\) and \(gq\rightarrow \gamma \gamma q\) amplitudes. In the one-gluon emission processes, corrections via an initial-state gluon emission, color-octet \(t\bar{t}\) effects, and an ultrasoft gluon emission from the on-shell \(t\bar{t}\) state appear. These corrections can be sizable because of the large partonic luminosity of the color-octet gluons. Investigations of these effects are left for future work. However, we expect that these would not severely spoil the characteristic shape of the spectrum in the resonance region, because the initial-state radiation does not affect the bound-state formation and the color-octet \(t\bar{t}\) Green function is known to have a smooth slope in the resonance region. Finally, it might be possible to determine \(\Gamma _t\) simultaneously with \(m_t\) at the FCC.

To conclude, we have studied the \(gg\rightarrow \gamma \gamma \) amplitude with the \(t\bar{t}\) bound-state effects near their mass threshold by collecting the available higher-order corrections in non-relativistic QCD. We have predicted a characteristic structure in the diphoton mass spectrum near the threshold whose shape is stable under the scale uncertainty, while the overall normalization has an uncertainty of 10% level due to the lack of the two-loop amplitude. We have proposed a new method to determine \(m_t\) from the diphoton mass spectrum at the LHC and the FCC. We have shown that the estimated statistical errors are fairly small at the FCC, which deserves further realistic experimental studies and also motivates one to calculate higher-order corrections in theory.