1 Introduction↩︎

The associated production of a photon with a top-antitop quark pair (\(t\bar{t}\gamma\)) is one of the important processes at hadron colliders and next-generation experimental facilities. Evidence for \(t\bar{t}\gamma\) production was first observed in \(p\bar{p}\) collisions at the Tevatron by the CDF Collaboration [1]. Subsequently, in \(pp\) collisions at the Large Hadron Collider (LHC), both the ATLAS and CMS Collaborations observed and studied \(t\bar{t}\gamma\) production [2]–[9]. The \(t\bar{t}\gamma\) process probes the top-photon interaction and provides a direct way to measure the top-quark electric charge [10], [11]. Any possible modification of the top-photon vertex due to new physics would directly affect the production of \(t\bar{t}\gamma\) (see, e.g., Refs. [12]–[15]), making precise measurements of this interaction an additional test of the Standard Model (SM). Furthermore, \(t\bar{t}\gamma\) process also plays an important role in investigating the top-quark charge asymmetry [16], [17].

Theoretical predictions for \(t\bar{t}\gamma\) production at next-to-leading order (NLO) in quantum chromodynamics (QCD) with stable top quarks were first presented in Refs. [18]–[20], and later studied including top quark decays in Ref. [21]. Matching to parton showers was addressed in Ref. [22]. Off-shell effects of the top quarks in this process were investigated in Refs. [23]–[26]. NLO electroweak (EW) corrections were reported in Ref. [27]. Subsequently, the so-called complete NLO corrections with stable top quarks were presented in an automated calculation in Ref. [28], and later extended to include top quark decay in Ref. [29]. More recently, soft gluon resummation was performed and used to derive the approximate next-to-next-to-leading order (NNLO) predictions including both NLO QCD and NLO EW corrections [30]. Beyond these developments, to match the unprecedented precision expected at current and future LHC runs, the full NNLO QCD corrections are essential and represent the current frontier. In these NNLO corrections, besides the two-loop amplitudes, higher order contributions from one-loop amplitudes in the dimensional regulator \(\epsilon\) are also required. Recently, Ref. [31] obtained the one-loop contributions up to \(\mathcal{O}(\epsilon^2)\) for both the quark-antiquark annihilation and gluon fusion channels. In addition, a class of two-loop planar master integrals relevant for \(t\bar{t}\gamma\) production was computed in Refs. [32]–[34]¹. Purely massless two-loop amplitudes have been considered in Ref. [35]. However, an exact calculation of the full two-loop amplitudes is still beyond the reach of current methods.

In this paper, we present approximate two-loop amplitudes for \(t\bar{t}\gamma\) production. We first compute the infrared (IR) poles of the two-loop amplitudes using the method presented in Refs. [36], [37], which incorporates the full mass dependence of the top quark. For the finite terms, instead of performing an exact calculation, we work in the high-energy boosted limit where all kinematic invariants are much larger than the top-quark mass, i.e., \(|s_{ij}|\gg m_t^2\) for all pairs of external legs \(i\neq j\). In this regime, a general mass-factorization formula valid at leading power (LP) has been established in Refs. [38]–[41]. To gain insight into power-suppressed effects—analogous to the study of \(t\bar{t}H\) production in Ref. [42]—we keep part of the top-quark mass dependence in the amplitudes, leading to what we call “semi-massive scheme”. Concretely, we start from the fully massive amplitudes but treat all scalar Feynman integrals as massless. In parallel, we construct amplitudes directly from LP mass factorization, which we refer to as the “massless scheme”. These two approaches thus provide complementary approximate two-loop amplitudes for \(t\bar{t}\gamma\) production. While significantly simpler than evaluating the fully massive amplitudes, both strategies are enabled by recent progress in integration-by-parts (IBP) reduction techniques [43]–[49] and in the computation of two-loop five-point master integrals [50]–[57].

The paper is organized as follows. In Sec. 2, we introduce the notation and define the quantities used in our calculation. In Sec. 3, we derive the IR poles of the two-loop amplitudes with full mass dependence and present approximate formulas for the amplitudes in the high-energy boosted limit. Section 4 discusses the treatment of power-suppressed contributions and the mapping between massive and massless phase-space points. Numerical results are shown in Sec. 5, and conclusions are given in Sec. 6. Lengthy expressions are collected in the appendices.

2 Notation and definitions↩︎

For photon production in association with a top-antitop quark pair, we consider the following two partonic processes: \[\begin{align} q_\beta(p_1)+\bar{q}_\alpha(p_2) &\to t_k(p_3)+\bar{t}_l(p_4)+\gamma(p_5) \,, \\ g_a(p_1)+g_b(p_2) &\to t_k(p_3)+\bar{t}_l(p_4)+\gamma(p_5) \,, \end{align}\] where \(\alpha, \beta, a, b, k, l\) are color indices, and \(p_i\) denote the momenta of the external partons with \(p_1^2=p_2^2=p_5^2=0\) and \(p_3^2=p_4^2=m_t^2\). The kinematic variables are defined as \(s_{ij}=(p_i + \sigma_{ij} p_j)^2\), where \(\sigma_{ij}=+1\) if \(p_i\) and \(p_j\) are both incoming or both outgoing, and \(\sigma_{ij}=-1\) otherwise.

To facilitate the calculation of the amplitudes, we use the color-space formalism [58], [59], in which the amplitudes are treated as vectors \(\ket{{\cal M}_{q,g}}\). The subscripts \(q\) or \(g\) denote the quark-antiquark annihilation (\(q\bar{q}\)) channel or the gluon fusion (\(gg\)) channel, respectively. For the \(q\bar{q}\) channel, we choose the independent color structures as \[\label{eq:qqbasis} \ket{c^q_1} = \delta_{\alpha\beta}\,\delta_{kl} \,, \qquad \ket{c^q_2} = (t^a)_{\alpha\beta}\,(t^a)_{kl} \,.\tag{1}\] For the \(gg\) channel, we adopt the color basis \[\label{eq:ggbasis} \ket{c^g_1} = \delta^{ab}\,\delta_{kl} \,, \qquad \ket{c^g_2} = if^{abc}\,(t^c)_{kl} \,, \qquad \ket{c^g_3} = d^{abc}\,(t^c)_{kl} \,.\tag{2}\] Note that \(\ket{c^q_I}\) and \(\ket{c^g_I}\) are orthogonal but not normalized. The amplitudes \(\ket{\mathcal{M}_{q,g}}\) can therefore be decomposed in color space as \[\begin{align} \ket{\mathcal{M}_{q,g}} = \sum_{I}d_I^{q,g} \ket{c^{q,g}_I}\,, \end{align}\] where \(d_I^{q,g}\) are the decomposition coefficients of the \(I\)-th color structure in the \(q\bar{q}\) or \(gg\) channel. We define the matrix elements of the squared amplitudes in this color space as \[\begin{align} \label{eq:matirxE0} \left(\boldsymbol{H}_{q,g}\right)_{IJ}&= \frac{\braket{c_I|{\cal M}_{q,g}}}{\braket{c_I|c_I}}\frac{\braket{{\cal M}_{q,g}|c_J}}{\braket{c_J|c_J}} \nonumber \\ &= d_{I,q,g}\left(d_{J,q,g}\right)^{\dagger} \,, \end{align}\tag{3}\] where the spin sum over external partons has been performed with the polarization sums for the two initial gluons and the final photon chosen as \[\begin{align} \sum_s\varepsilon^\mu(p_1,s)\varepsilon^\nu(p_1,s)&=\sum_s\varepsilon^\mu(p_2,s)\varepsilon^\nu(p_2,s)=-g^{\mu\nu}+\frac{p_1^\mu p_2^\nu+p_1^\nu p_2^\mu}{p_1\cdot p_2} \,, \nonumber \\ \sum_s\varepsilon^\mu(p_5,s)\varepsilon^\nu(p_5,s)&=-g^{\mu\nu} \,. \end{align}\] The squared amplitudes in terms of \(\left(\boldsymbol{H}_{q,g}\right)_{IJ}\) are then given by \[\begin{align} \label{eq:sqamp0} \braket{{\cal M}_{q,g}|{\cal M}_{q,g}}=\sum_{IJ}\left(\boldsymbol{H}_{q,g}\right)_{IJ}\braket{c_J|c_I} \,. \end{align}\tag{4}\]

The ultraviolet (UV) divergences in the bare amplitudes \(\ket{\mathcal{M}_{q,g}}\) are removed by renormalization, \[\begin{align} \Ket{{\cal M}^R_{q,g}(\alpha_s, e, m_t, \mu, \epsilon)} = \left(\frac{\mu^2 e^{\gamma_E}}{4\pi}\right)^{-\epsilon} Z_{q,g} Z_Q \Ket{{\cal M}^{\text{bare}}_{q,g}(\alpha_s^0, e^0, m_t^0,\epsilon)} , \label{eq:UVren} \end{align}\tag{5}\] where \(Z_g\), \(Z_q\), and \(Z_Q\) are the on-shell wave-function renormalization constants for gluons, light quarks, and the heavy (top) quark, respectively, and dependence on other kinematic variables is left implicit. The top-quark mass is renormalized in the on-shell scheme, \(m_t^0 = Z_m m_t\), while the strong coupling constant \(\alpha_s\) is renormalized in the \(\overline{\text{MS}}\) scheme with \(n_f=n_l+n_h\) active flavors: \[\label{eq:renorcoupling} \alpha_s^0 = \left( \frac{\mu^2 e^{\gamma_E}}{4\pi} \right)^{\epsilon} Z_{\alpha_s} \alpha_s \, .\tag{6}\] Since we restrict ourselves to QCD corrections, the quantum electrodynamics (QED) coupling does not require renormalization, and we simply set \(e^0=e\). Explicit expressions for the renormalization constants up to NNLO are collected in Appendix 7.

After UV renormalization, the remaining IR divergences can be subtracted by a multiplicative factor \(\boldsymbol{Z}^{-1}(\epsilon, m_t, \mu)\), where the boldface symbol denotes an operator in color space. Explicitly, we have \[\begin{align} \label{eq:key} \boldsymbol{Z}_{q,g}^{-1}(\alpha_s, m_t, \mu, \epsilon) \Ket{{\cal M}^R_{q,g}(\alpha_s, e, m_t, \mu, \epsilon)} = \text{finite} \,. \end{align}\tag{7}\] Both the UV-renormalized amplitudes and the IR subtraction factors can be expanded in powers of \(\alpha_s\): \[\begin{align} \Ket{{\cal M}_{q,g}^R} &= 4\pi\alpha_s Q_te \left[ \Ket{{\cal M}_{q,g}^{(0)}} + \frac{\alpha_s}{4\pi} \Ket{{\cal M}_{q,g}^{(1)}} + \left( \frac{\alpha_s}{4\pi} \right)^2 \Ket{{\cal M}_{q,g}^{(2)}} + \cdots \right] , \nonumber \\ \boldsymbol{Z}_{q,g} &= \boldsymbol{1} + \frac{\alpha_s}{4\pi} \boldsymbol{Z}_{q,g}^{(1)} + \left( \frac{\alpha_s}{4\pi} \right)^2 \boldsymbol{Z}_{q,g}^{(2)} + \cdots \, , \end{align}\] where \(Q_t = 2/3\) is the electric charge of the top quark. The IR subtraction factors up to NNLO can be found, for example, in Ref. [37], and read \[\begin{align} \boldsymbol{Z}^{(1)}_{q,g} &= \frac{\Gamma_0'^{q,g}}{4\epsilon^2} + \frac{\boldsymbol{\Gamma}^{q,g}_0}{2\epsilon} \,,\tag{8}\\ \boldsymbol{Z}^{(2)}_{q,g} &= \frac{(\Gamma_0'^{q,g})^2}{32\epsilon^4} + \frac{\Gamma_0'^{q,g}}{8\epsilon^3} \left( \boldsymbol{\Gamma}^{q,g}_0 - \frac{3}{2}\,\beta^{n_l}_0 \right) + \frac{\boldsymbol{\Gamma}^{q,g}_0}{8\epsilon^2} \left( \boldsymbol{\Gamma}^{q,g}_0 -2\beta^{n_l}_0 \right) + \frac{\Gamma_1'^{q,g}}{16\epsilon^2} + \frac{\boldsymbol{\Gamma}^{q,g}_1}{4\epsilon} \nonumber \\ &\qquad - \frac{2T_F}{3}n_h\, \Bigg[ \Gamma_0'^{q,g} \left( \frac{1}{2\epsilon^2}\,\ln\frac{\mu^2}{m_t^2} + \frac{1}{4\epsilon} \left[ \ln^2\!\frac{\mu^2}{m_t^2} + \frac{\pi^2}{6} \right] \right) + \frac{\boldsymbol{\Gamma}_0^{q,g}}{\epsilon}\,\ln\frac{\mu^2}{m_t^2} \Bigg] \,, \tag{9} \end{align}\] where the coefficients \(\boldsymbol{\Gamma}_n\) and \(\Gamma_n'\) are defined via the expansion \[\begin{align} \boldsymbol{\Gamma} &= \sum_{n\ge 0}\,\left( \frac{\alpha_s}{4\pi} \right)^{n+1} \, \boldsymbol{\Gamma}_n \,, \qquad \Gamma_n' = \frac{\partial\boldsymbol{\Gamma}_n}{\partial\ln\mu}\,. \end{align}\] For \(t\bar{t}\gamma\) production, the anomalous-dimension matrices \(\boldsymbol{\Gamma}^{q,g}\) up to \(\alpha_s^2\) in the color bases of Eqs. 1 and 2 coincide with those of \(t\bar{t}H\) process [42], [60]. We emphasize that in Eq. 9 , \(\beta^{n_l}_0 = 11\,C_A/3-4\,T_F n_l/3\) is the \(\beta\)-function coefficient, with \(n_l\) denoting the number of active light flavors. The dependence on \(n_h\) is entirely contained in the last line of Eq. 9 .

With these ingredients, the IR singularities of the amplitudes can be extracted order by order in \(\alpha_s\): \[\begin{align} \label{eq:IRamps} \Ket{{\cal M}_{q,g}^{(1)\,,R,\,\text{sing}}} &= \boldsymbol{Z}_{q,g}^{(1)} \Ket{{\cal M}_{q,g}^{(0)R}} , \nonumber \\ \Ket{{\cal M}_{q,g}^{(2)\,,R,\,\text{sing}}} &= \left[ \boldsymbol{Z}_{q,g}^{(2)} - \left( \boldsymbol{Z}_{q,g}^{(1)R}\right)^2 \right] \Ket{{\cal M}_{q,g}^{(0)R}} + \left( \boldsymbol{Z}_{q,g}^{(1)} \Ket{{\cal M}_{q,g}^{(1)R}} \right)_{\text{poles}} \,. \end{align}\tag{10}\] In particular, predicting the IR poles at two-loop order requires computing the UV-renormalized one-loop amplitudes up to \({\cal O}(\epsilon^1)\).

Having established the structure of the IR poles of the amplitudes, we now turn to the main objects of study in this work: the two-loop renormalized amplitudes \(\ket{\mathcal{M}^{(2)R}_{q,g}}\). Their singular parts follow directly from Eq. 10 , but the finite remainders require evaluating genuinely two-loop five-point Feynman integrals involving six physical scales—a task that is currently beyond exact analytic methods. To make progress, we adopt the high-energy boosted limit and employ the mass-factorization approach [38] to obtain approximate expressions for the finite parts. There is, however, a special class of two-loop diagrams with two fermion loops that is comparatively simple: the corresponding integrals can be factorized into products of two one-loop integrals.² These contributions correspond to the \(n_l^2\), \(n_h^2\), and \(n_ln_h\) terms in the color decomposition of Eq. 23 , which are collectively denoted as \(n_f^2\) terms in the following. We compute these contributions exactly, thereby providing both a check on the IR structure of \(\ket{\mathcal{M}^{(2)R}_{q,g}}\) inferred from Eq. 10 and a useful benchmark for our approximate results in the high-energy boosted limit.

3 Calculation↩︎

3.1 Massive amplitudes at one loop↩︎

As discussed in the previous section, predicting the IR poles at two loops requires the full massive one-loop amplitudes expanded up to \(\mathcal{O}\left(\epsilon^1\right)\). We generate the amplitudes with full top-quark mass dependence using FeynArts [61], and manipulate the resulting expressions with FeynCalc [62], [63]. The amplitudes then need to be expressed in terms of scalar integrals. Since our main interest lies in the IR poles of the interference between \(\mathcal{M}^{(2)}\) and \(\mathcal{M}^{(0)}\) via Eq. 10 , we can directly multiply the one-loop amplitudes by the tree-level ones. This allows the Lorentz contractions and Dirac traces to be carried out straightforwardly, while still keeping track of the color information necessary for predicting the two-loop amplitudes in the high-energy boosted limit according to the mass-factorization formula Eq. 15 . Alternatively, one may apply a complete set of projectors to the one-loop amplitudes and extract the coefficients as linear combinations of scalar integrals. Although this method is more involved for our purposes, it is useful for obtaining the one-loop squared amplitude. We have performed the calculation using both approaches, and found consistent results.

The one-loop scalar integrals can be categorized into \(17\) families (topologies) neglecting permutations of the external momenta among \(p_1\), \(p_2\), and \(p_5\), and/or between \(p_3\) and \(p_4\). Each family is specified by five propagators. The \(j\)-th scalar integral in the \(i\)-th family is defined as \[\begin{align} \label{eq:def95integral951} F^{(1)}_i\left[\vec{a}_j\right] = \int \frac{d^dl_1}{(2\pi)^d} \, f^{(1)}_i\left[\vec{a}_j\right]\,, \qquad f^{(1)}_i\left[\vec{a}_j\right] = \frac{1}{D_{i,1}^{a_{j,1}} D_{i,2}^{a_{j,2}} D_{i,3}^{a_{j,3}} D_{i,4}^{a_{j,4}} D_{i,5}^{a_{j,5}}}\,, \end{align}\tag{11}\] where \(f^{(1)}_i\left[\vec{a}_j\right]\) denotes the integrand of the \(j\)-th one-loop Feynman integral in the \(i\)-th family, and \(D_{i,k}\) (\(k=1,\ldots,5\)) are the corresponding propagator denominators. The choice of \(D_{i,k}\) for each family is not unique. For instance, for the first family we define \[\begin{align} &T^{(1)}_1=\left\{l_1^2,\; (l_1-p_1)^2-m_t^2,\;(l_1-p_1-p_2)^2,\;(l_1-p_1-p_2+p_3)^2, \;(l_1-p_5)^2\right\} \,. \end{align}\] The definitions of all \(17\) families are provided in ancillary electronic files attached with this paper. According to Eq. 3 , the NLO matrix elements of the bare squared amplitudes can be written as \[\begin{align} \label{eq:sqamp1} \left(\boldsymbol{H}^{(1)\text{bare}}_{q,g}\right)_{IJ}&= \left(d^{(0)\text{bare}}_{J,q,g}\right)^{\dagger}d^{(1)\text{bare}}_{I,q,g} \nonumber \\ &=\int \frac{d^dl_1}{(2\pi)^d}\sum_{ij}c^{(1)\text{bare}}_{ij,IJ,q,g}f^{(1)}_i\left[\vec{a}_j\right]\,. \end{align}\tag{12}\]

After applying IBP reduction with the package Kira [64], [65], we identify \(107\) master integrals (MIs), including mappings among the MIs of all families.³ Up to linear transformations arising from different definitions of integral families and choices of MIs, all of these integrals have already been evaluated numerically in the literature using the method of differential equations [31], [66]. In this work, we calculate all \(107\) MIs numerically with the package AMFlow [67], which allows us to obtain the amplitudes up to \(\mathcal{O}(\epsilon^2)\), as required for the NNLO calculation. We have also compared our NLO squared amplitudes up to \(\mathcal{O}(\epsilon^0)\) with those produced by the package GoSam [68], and we find complete agreement, modulo the choice of dimensional regularization scheme.⁴

3.2 Massive amplitudes at two loops↩︎

To obtain approximate expressions for the finite terms of the two-loop amplitudes using the mass-factorization approach [38] in the high-energy boosted limit, we require the amplitudes up to two-loop order with all parton masses set to zero. For generality, we first generate the two-loop amplitudes with full mass dependence using FeynArts [61] and process the expressions with FeynCalc [62], [63]. The amplitudes are then expressed in terms of scalar integrals by multiplying the two-loop amplitudes by the tree-level ones.

Similar to the one-loop case, the two-loop scalar integrals can be classified into \(467\) families (before accounting for permutations of the external momenta). Each family is defined by \(11\) propagator denominators. The \(j\)-th scalar integral in the \(i\)-th family is defined as \[\begin{align} \label{eq:def95integral952} F^{(2)}_i\left[\vec{a}_j\right] = \int \frac{d^dl_1}{(2\pi)^d}\frac{d^dl_2}{(2\pi)^d} f^{(2)}_i\left[\vec{a}_j\right] \,, \qquad f^{(2)}_i\left[\vec{a}_j\right]=\prod_{r=1}^{11}\frac{1}{D_{i,r}^{a_{j,r}} }\, , \end{align}\tag{13}\] where \(f^{(2)}_i\left[\vec{a}_j\right]\) denotes the integrand of the \(j\)-th two-loop scalar integral in the \(i\)-th family. For the first family, the propagator denominators are chosen as \[\begin{align} &T^{(2)}_1=\big\{l_1^2,\;\left(l_1-p_1\right)^2,\;\left(l_1+p_2-p_3-p_4\right)^2,\;\left(l_1-p_3-p_4\right)^2,\;l_2^2,\;\left(l_2-p_4\right)^2-m_t^2, \nonumber \\ &\qquad\quad\;\left(l_1-l_2\right)^2,\; \left(l_1-p_4\right)^2-m_t^2,\;\left(l_2-p_1\right)^2,\;\left(l_2+p_2-p_3-p_4\right)^2,\;\left(l_2-p_3-p_4\right)^2 \big\} \,. \end{align}\] The definitions of all \(467\) families are provided in ancillary electronic files accompanying this paper. As a result, the matrix elements of the interference between \(\mathcal{M}^{(2)}\) and \(\mathcal{M}^{(0)}\) can be written as \[\begin{align} \label{eq:sqamp2} \left(\boldsymbol{H}^{(2)\text{bare}}_{q,g}\right)_{IJ}&= \left(d^{(0)\text{bare}}_{J,q,g}\right)^{\dagger}d^{(2)\text{bare}}_{I,q,g} \nonumber \\ &= \int\frac{d^dl_1}{(2\pi)^d}\frac{d^dl_2}{(2\pi)^d}\sum_{ij}c^{(2)\text{bare}}_{ij,IJ,q,g}f^{(2)}_i\left[\vec{a}_j\right] \,. \end{align}\tag{14}\] Without further simplification, we find that roughly \(1.7 \times 10^6\) two-loop scalar integrals \(F^{(2)}_i\left[\vec{a}_j\right]\) are needed to express \(\boldsymbol{H}^{(2)\text{bare}}_{q,g}\) for both the \(q\bar{q}\) and \(gg\) channels. These two-loop Feynman integrals involve six scales and are extremely challenging, rendering an exact calculation of the two-loop amplitudes for \(t\bar{t}\gamma\) production beyond the reach of current methods. Apart from the IR singularities of the two-loop amplitudes, which can be predicted from Eq. 10 , in what follows we employ the mass-factorization approach to approximate their finite terms.

3.3 The massive amplitudes from mass factorization↩︎

In this section, we present approximate amplitudes for \(t\bar{t}\gamma\) production in the high-energy boosted limit. To this end, we adopt the mass-factorization formula [38]–[41]: \[\begin{gather} \label{eq:QCDfac1} \Ket{\bar{\mathcal{M}}^R_{q,g}\left(\epsilon,\{p\},m_t,\mu\right)} = \mathcal{Z}_{[q,g]}^{(m|0)}(\epsilon,m_t,\mu) \, \mathcal{Z}_{[t]}^{(m|0)}(\epsilon,m_t,\mu) \\ \times \boldsymbol{\mathcal{S}}(\epsilon,\{\tilde{p}\},m_t,\mu)\Ket{\tilde{\mathcal{M}}^R_{q,g}\left(\epsilon,\{\tilde{p}\},\mu\right)} \,, \end{gather}\tag{15}\] where \(\tilde{\mathcal{M}}^R\) is the renormalized massless amplitude with all parton masses set to zero, and \(\bar{\mathcal{M}}^R\) approximates the full amplitude \(\mathcal{M}^R\) up to power-suppressed corrections in \(m_t^2\).⁵ The \(\mathcal{Z}\)-factors \(\mathcal{Z}_{[j]}^{(m|0)}\) and the soft function \(\boldsymbol{\mathcal{S}}\) contain all logarithmic mass dependence, \(\ln(\mu^2/m_t^2)\). General results for \(\mathcal{Z}_{[j]}^{(m|0)}\) and \(\boldsymbol{\mathcal{S}}\) up to NNLO can be found in Refs. [38]–[41], [73], [74], while the specific expressions relevant for \(t\bar t \gamma\) coincide with those for \(t\bar t H\) production and are given in Appendix C of Ref. [42]. We use \(\{p\}\) and \(\{\tilde{p}\}\) to denote the sets of external momenta in the massive and massless cases, respectively. The momenta in \(\tilde{\mathcal{M}}^R\) satisfy \(\tilde{p}_i^2 = 0\), and we define \(\tilde{s}_{ij} \equiv 2\sigma_{ij} \, \tilde{p}_i \cdot \tilde{p}_j\). In the high-energy boosted limit, \(s_{ij}\approx\tilde{s}_{ij}\). There remains some freedom in how the massless momenta \(\{\tilde{p}\}\) are related to the massive ones \(\{p\}\). Before specifying our choice, we turn to the calculation of the two-loop massless amplitude \(\tilde{\mathcal{M}}^R\).

Analogous to the massive case, we renormalize the massless amplitudes as \[\begin{align} \label{eq:renormalmassless} \Ket{{\tilde{\mathcal{M}}}^R_{q,g}(\alpha_s, e, \{\tilde{p}\}, \mu, \epsilon)} = \left(\frac{\mu^2 e^{\gamma_E}}{4\pi}\right)^{-\epsilon} \Ket{\tilde{\mathcal{M}}^{\text{bare}}_{q,g}(\alpha_s^0, e^0, \{\tilde{p}\}, \epsilon)} \,. \end{align}\tag{16}\] In the massless case, the on-shell wave-function renormalization constants for gluons and quarks are unity, and the renormalization of \(\alpha_s\) is the same as in Eq. 6 . The bare massless amplitudes \(\tilde{\mathcal{M}}^{\text{bare}}_{q,g}(\{\tilde{p}\})\) can be obtained by taking the massless limit of the massive ones \(\mathcal{M}^{\text{bare}}_{q,g}(\{p\}, m_t^0)\) at the integrand level, i.e. before loop integration. Consequently, the N\(^k\)LO (\(k=1,2\)) matrix elements of the bare massless squared amplitudes can be written analogously to Eq. 14 :⁶ \[\begin{align} \label{eq:sqamp3} \left(\tilde{\boldsymbol{H}}^{(k)\text{bare}}_{q,g}\right)_{IJ}= \int\left(\prod_{r=1}^{k}\frac{d^dl_r}{(2\pi)^d}\right)\sum_{ij}\tilde{c}^{(k)\text{bare}}_{ij,IJ,q,g}\tilde{f}^{(k)}_i\left[\vec{a}_j\right]\,, \end{align}\tag{17}\] with \[\begin{align} \label{eq:transmass} \tilde{c}^{(k)\text{bare}}_{ij,IJ,q,g}=\left.c^{(k)\text{bare}}_{ij,IJ,q,g}\right|_{m_t\to 0, \,p_i\to \tilde{p}_i} \,, \quad \tilde{f}^{(k)}_i\left[\vec{a}_j\right]=\left.f^{(k)}_i\left[\vec{a}_j\right]\right|_{m_t\to 0, \,p_i\to \tilde{p}_i}\,. \end{align}\tag{18}\] Note that \(\tilde{c}^{(k)\text{bare}}_{ij,IJ,q,g}\) and \(\tilde{f}^{(k)}_i\left[\vec{a}_j\right]\) are exactly the LP contributions of \(c^{(k)\text{bare}}_{ij,IJ,q,g}\) and \(f^{(k)}_i\left[\vec{a}_j\right]\) in the high-energy boosted limit, respectively.

The NNLO massless squared amplitudes can therefore be expressed as linear combinations of two-loop massless scalar Feynman integrals. In this case, the number of two-loop scalar integrals is reduced to about \(7.8 \times 10^5\), compared to the massive case. These scalar integrals are categorized into \(187\) integral families, which reduce to four independent topologies after accounting for permutations of external momenta. We carry out the IBP reduction using the package Kira [64], [65] with the help of FireFly [75], [76], expressing all scalar integrals in terms of master integrals⁷ studied in Refs. [50]–[57], which can be evaluated numerically using PentagonMI [77].⁸ Using the same program, we also compute the one-loop UT bases up to weight \(4\) required for the one-loop amplitudes up to \(\mathcal{O}\left(\epsilon^2\right)\). Finally, we validate the massless squared amplitudes by checking their IR poles against those predicted by the massless IR subtraction factors [78] and the corresponding lower-order squared amplitudes. We find excellent agreement within numerical accuracy.

4 Approximate massive squared amplitudes and mapping between massive and massless phase-space points↩︎

As noted in Eq. 15 , the factorization formula is valid at LP and receives power-suppressed corrections. Consequently, ambiguities of \(\mathcal{O}(m_t^2)\) arise, depending on the specific convention adopted to implement the factorization. In particular, when applying the formula to compute physical cross sections, two types of ambiguities occur, which we discuss in this section.

4.1 Approximate massive squared amplitudes↩︎

The first ambiguity concerns how the factorization formula is applied to compute squared amplitudes. Since the tree-level and one-loop massive amplitudes can be calculated exactly, the only ambiguity in the NNLO squared amplitudes lies in the interference between the two-loop and tree-level amplitudes. To validate our approximate approach, we also present results at NLO in the high-energy boosted limit. We consider two representative schemes and compare their outcomes, while keeping in mind that other choices are possible.

In the first scheme, which we refer to as the “massless scheme”, we directly square the approximate amplitudes \(\ket{\bar{\mathcal{M}}^R_{q,g}}\) given in Eq. 15 and extract the NLO and NNLO terms. Here we only consider the interferences of the amplitudes with the tree-level ones, namely \(\braket{\bar{\mathcal{M}}^{(0)R}_{q,g}|\bar{\mathcal{M}}^{R}_{q,g}}\). Note that \(\ket{\bar{\mathcal{M}}^{(0)R}_{q,g}}\) are identical to the massless tree-level amplitudes \(\ket{\tilde{\mathcal{M}}^{(0)R}_{q,g}}\). All mass dependence is contained in the factors \(\mathcal{Z}_{[j]}^{(m|0)}\) and the soft function \(\boldsymbol{\mathcal{S}}\). Combining Eqs. 4 , 15 , 16 , and 17 , we obtain \[\begin{align} \label{eq:massless95scheme} \braket{\bar{\mathcal{M}}^{(0)R}_{q,g}|\bar{\mathcal{M}}^{R}_{q,g}}&=\mathcal{Z}_{[q,g]}^{(m|0)}(\epsilon,m_t,\mu) \, \mathcal{Z}_{[t]}^{(m|0)}(\epsilon,m_t,\mu) \nonumber \\ &\times \sum_{IJK}\Big[\boldsymbol{\mathcal{S}}(\epsilon,\{\tilde{p}\},m_t,\mu)\Big]^{q,g}_{KI}\Big[\tilde{\boldsymbol{H}}^R_{q,g}(\epsilon,\{\tilde{p}\},\mu)\Big]_{IJ}\braket{c_J|c_K} \,, \end{align}\tag{19}\] where \[\begin{align} \label{eq:softfunctionME} \boldsymbol{\mathcal{S}}^{q,g}_{KI} = \frac{\Bra{c^{q,g}_K}\boldsymbol{\mathcal{S}}\Ket{c^{q,g}_I}}{\braket{c^{q,g}_K|c^{q,g}_K}}\,. \end{align}\tag{20}\] Although the soft function \(\boldsymbol{\mathcal{S}}\) has the same operator form in both the \(q\bar{q}\) and \(gg\) channels, the matrix elements in Eq. 20 differ between the partonic channels due to their distinct color structures [41]. To extract the NLO and NNLO contributions, we require \(\tilde{\boldsymbol{H}}^R_{q,g}\) up to NNLO. At LO, \[\begin{align} \Big[\tilde{\boldsymbol{H}}^{(0)R}_{q,g}\Big]_{IJ} = \Big[\tilde{\boldsymbol{H}}^{(0)\text{bare}}_{q,g}\Big]_{IJ} =\left(\tilde{d}^{(0)\text{bare}}_{J,q,g}\right)^{\dagger}\tilde{d}^{(0)\text{bare}}_{I,q,g}\,, \end{align}\] and the higher order terms \(\tilde{\boldsymbol{H}}^{(k)R}_{q,g}\) (\(k=1,2\)) follow from Eq. 17 after UV renormalization, which are given by \[\begin{align} \tilde{\boldsymbol{H}}^{(1)R}_{q,g} &= \left(\frac{\mu^2e^{\gamma_E}}{4\pi}\right)^{\epsilon}\tilde{\boldsymbol{H}}^{(1)\text{bare}}_{q,g} + Z^{(1)}_{\alpha_s}\tilde{\boldsymbol{H}}^{(0)\text{bare}}_{q,g}\,, \nonumber \\ \tilde{\boldsymbol{H}}^{(2)R}_{q,g} &= \left(\frac{\mu^2e^{\gamma_E}}{4\pi}\right)^{2\epsilon}\tilde{\boldsymbol{H}}^{(2)\text{bare}}_{q,g} + 2\left(\frac{\mu^2e^{\gamma_E}}{4\pi}\right)^{\epsilon}Z^{(1)}_{\alpha_s}\tilde{\boldsymbol{H}}^{(1)\text{bare}}_{q,g} + Z^{(2)}_{\alpha_s}\tilde{\boldsymbol{H}}^{(0)\text{bare}}_{q,g}\,. \end{align}\]

The second scheme is slightly more involved. We introduce a modified version of Eq. 19 : \[\begin{align} \label{eq:massive95scheme} \braket{\hat{\mathcal{M}}^{(0)R}_{q,g}|\hat{\mathcal{M}}^{R}_{q,g}}&=\mathcal{Z}_{[q,g]}^{(m|0)}(\epsilon,m_t,\mu) \, \mathcal{Z}_{[t]}^{(m|0)}(\epsilon,m_t,\mu) \nonumber \\ &\times \sum_{IJK}\Big[\boldsymbol{\mathcal{S}}(\epsilon,\{\tilde{p}\},m_t,\mu)\Big]_{KI}\Big[\hat{\boldsymbol{H}}^R_{q,g}(\epsilon,\{\tilde{p}\},\{p\},m_t,\mu)\Big]_{IJ}\braket{c_J|c_K}\,, \end{align}\tag{21}\] where the modified matrices \(\hat{\boldsymbol{H}}^R_{q,g}\) are defined up to NNLO as \[\begin{align} \Big[\hat{\boldsymbol{H}}^{(0)R}_{q,g}\Big]_{IJ} &= \Big[\boldsymbol{H}^{(0)R}_{q,g}\Big]_{IJ} =\left(d^{(0)\text{bare}}_{J,q,g}\right)^{\dagger}d^{(0)\text{bare}}_{I,q,g}\,,\nonumber \\ \hat{\boldsymbol{H}}^{(1)R}_{q,g} &= \left.\left(\frac{\mu^2e^{\gamma_E}}{4\pi}\right)^{\epsilon}\boldsymbol{H}^{(1)\text{bare}}_{q,g}\right|_{f^{(1)}_i\left[\vec{a}_j\right]\to \tilde{f}^{(1)}_i\left[\vec{a}_j\right]} + Z^{(1)}_{\alpha_s}\boldsymbol{H}^{(0)\text{bare}}_{q,g}\,, \nonumber \\ \hat{\boldsymbol{H}}^{(2)R}_{q,g} &= \left.\left(\frac{\mu^2e^{\gamma_E}}{4\pi}\right)^{2\epsilon}\boldsymbol{H}^{(2)\text{bare}}_{q,g}\right|_{f^{(2)}_i\left[\vec{a}_j\right]\to \tilde{f}^{(2)}_i\left[\vec{a}_j\right]} + \left.2\left(\frac{\mu^2e^{\gamma_E}}{4\pi}\right)^{\epsilon}Z^{(1)}_{\alpha_s}\boldsymbol{H}^{(1)\text{bare}}_{q,g}\right|_{f^{(1)}_i\left[\vec{a}_j\right]\to \tilde{f}^{(1)}_i\left[\vec{a}_j\right]} \nonumber \\ &+ Z^{(2)}_{\alpha_s}\boldsymbol{H}^{(0)\text{bare}}_{q,g}\,, \end{align}\] with \(\boldsymbol{H}^{(1)\text{bare}}_{q,g}\) and \(\boldsymbol{H}^{(2)\text{bare}}_{q,g}\) defined in Eqs. 12 and 14 , respectively. The massive integrands \(f^{(1)}_i\left[\vec{a}_j\right]\) and \(f^{(2)}_i\left[\vec{a}_j\right]\) are given in Eqs. 11 and 13 , while their massless counterparts \(\tilde{f}^{(k)}_i\left[\vec{a}_j\right]\) are defined in Eq. 18 . In other words, the modified matrices \(\hat{\boldsymbol{H}}^R_{q,g}\) are obtained by replacing all scalar integrals in the bare massive matrices \(\boldsymbol{H}^{\text{bare}}_{q,g}\) with their massless counterparts, and then performing UV renormalization as in the massless case. Thus, the only difference between Eqs. 21 and 19 arises from the distinction between \(\hat{\boldsymbol{H}}^R_{q,g}\) and \(\tilde{\boldsymbol{H}}^R_{q,g}\), which is power suppressed. In this setup, \(\hat{\boldsymbol{H}}^R_{q,g}\) retains full mass information from the external top-quark pair and partial mass information from internal top quark propagators. We therefore refer to this scheme as the “semi-massive scheme”.

4.2 Mapping between massive and massless phase-space points↩︎

The second ambiguity concerns the relation between the massive external momenta \(p_i\) and the corresponding massless ones \(\tilde{p}_i\) in Eq. 15 . Following the prescription used for \(t\bar{t}H\) production [42], we keep the directions of the 3-momenta identical between \(p_i\) and \(\tilde{p}_i\) and rescale their magnitudes (and energy components) such that the massless on-shell conditions and momentum conservation are both satisfied.

Specifically, since \(p_1\) and \(p_2\) are already massless, we simply set \(\tilde{p}_1=p_1\) and \(\tilde{p}_2=p_2\) resulting in \(\tilde{s}_{12}=s_{12}\). We parameterize the massive external momenta in the center-of-mass frame as \[\begin{align} \label{eq:momentamassive} p_1&=\frac{\sqrt{s_{12}}}{2}\left(1,0,0,1\right)\,,\nonumber\\ p_2&=\frac{\sqrt{s_{12}}}{2}\left(1,0,0,-1\right)\,,\nonumber\\ p_3&=\left(\sqrt{m_t^2+q_3^2},q_3\sin\theta_3\sin\phi_3,q_3\sin\theta_3\cos\phi_3,q_3\cos\theta_3\right)\,,\nonumber\\ p_4&=\left(\sqrt{m_t^2+q_4^2},q_4\sin\theta_4\sin\phi_4,q_4\sin\theta_4\cos\phi_4,q_4\cos\theta_4\right)\,,\nonumber\\ p_5&=q_5\left(1,\sin\theta_5,0,\cos\theta_5\right) \,, \end{align}\tag{22}\] where \(q_i\) is the magnitude of the spatial component of \(p_i\), \(\theta_i\) the polar angle, and \(\phi_i\) the azimuthal angle. For convenience and without loss of generality, we have set \(\phi_5=0\). Momentum conservation and on-shell conditions leave six independent parameters on the right-hand side of Eq. 22 , which we choose as \(s_{12}\), \(m_t\), \(q_5\), \(\theta_3\), \(\phi_3\), and \(\theta_5\). Their values satisfy the physical constraints \[\begin{align} &s_{12} \ge 4m_t^2\,, \quad 0\le q_5\le q_{5,\rm max}=\frac{s_{12}-4m_t^2}{2\sqrt{s_{12}}}\,,\nonumber\\ m_t >&0 \,, \quad 0\le \theta_3 \le \pi\,, \quad 0\le \phi_3 < 2\pi\,, \quad 0\le \theta_5 \le \pi \,. \end{align}\]

The massless momenta \(\tilde{p}_i\) are then obtained by rescaling \(\tilde{p}_{3,4,5}\) at fixed angles to enforce \(\tilde{p}_i^2=0\) and momentum conservation, following the same prescription as in \(t\bar{t}H\) production [42].

5 Numeric results↩︎

We are now ready to present numerical results for the two-loop amplitudes. For the input parameters, we take the top-quark mass \(m_t = \SI{173}{\GeV}\) and, by default, set the renormalization scale to \(\mu = m_t\). Other choices of \(\mu\) are of course possible and will be considered below.⁹

In addition, the amplitudes depend on the electric charges of the quarks, entering both through the external states and through closed internal quark loops. As a result, in the \(q\bar{q}\) channel, we must distinguish between up-type (e.g., \(u\bar{u}\)) and down-type (e.g., \(d\bar{d}\)) initial states. For contributions from closed light-quark loops with a photon attached, we likewise need to distinguish between up-type quarks with electric charge \(2/3\) and down-type quarks with charge \(-1/3\). In both cases the contributions from down-type light quarks can be obtained from those of up-type light quarks by including a factor of \(-1/2\) for each \(q\bar{q}\gamma\) vertex, without introducing additional computational complexity. For convenience, and without loss of generality, we therefore take the electric charge of all quarks to be \(Q_f = 2/3\) when presenting numerical results, which is expected to have only a mild impact on assessing the validity of our approximate approach.

To test the validity of the approximation, we first compare the approximate results with the exact ones. The exact results include the NLO squared amplitudes up to \(\mathcal{O}(\epsilon^2)\), and the IR poles of the NNLO squared amplitudes predicted by Eq. 10 . The NNLO squared amplitudes receive two types of contributions: the one-loop-squared amplitudes, which can be obtained exactly, and the interference between the two-loop and tree-level amplitudes. For clarity, we therefore drop the one-loop-squared contribution in the comparison and focus only on the interference terms.

We first show in Fig. 1 the numerical results for the squared amplitudes at \(\mathcal{O}(\epsilon^0)\) (NLO, left plots) and \(\mathcal{O}(\epsilon^{-1})\) (NNLO, right plots) in both the \(q\bar{q}\) (upper plots) and \(gg\) (lower plots) channels. The partonic center-of-mass energy \(\sqrt{s_{12}}\) is varied from \(\SI{1.5}{\TeV}\) to \(\SI{20}{\TeV}\), while the remaining independent kinematic variables are fixed as \(\theta_3 = 14\pi/29\), \(\phi_3=34\pi/29\), \(\theta_5=15\pi/29\), and \(q_5=20\,q_{5,\rm max}/29\). The black lines represent the exact results, while the red and orange curves denote the approximate results obtained using the two different schemes discussed in Sec. 4. As expected, the approximation improves with increasing energy for both NLO and NNLO. We observe that the semi-massive scheme (red curves) performs slightly better than the massless scheme (orange curves) at NLO, whereas at NNLO the opposite is true. This behavior can be traced to the fact that the semi-massive scheme incorporates part of the power-suppressed contributions (\(\mathcal{O}(m_t^2)\) and higher). However, this should be treated with caution, since beyond-LP terms can only be reliable once the unknown power-suppressed contributions from the scalar Feynman integrals are included. Overall, the massless scheme provides a reasonable approximation to the exact results down to \(\sqrt{s_{12}} \sim \SI{2}{\TeV}\) in all cases.

While the high-energy behavior is consistent with expectations, the different performance in the relatively low-energy region is more intriguing. From Fig. 1, we observe that at NLO the approximation behaves similarly in both the \(gg\) and \(q\bar{q}\) channels, whereas at NNLO it performs significantly better in the \(gg\) channel. In the \(q\bar{q}\) case, the relative deviation between the massless-scheme approximation and the exact result is much larger than the naive estimate of order \(m_t^2/|s_{ij}|\). Moreover, the difference between the semi-massive and massless schemes is also non-negligible. These effects can be attributed to accidental cancellations at LP in the small-mass limit, which making beyond-LP contributions more important than anticipated. A similar effect already appears at LO: the LO squared amplitude is a rational function of \(m_t^2\) and \(s_{ij}\) and can be straightforwardly expanded in the small-mass limit. We find cancellations at LP for the contributions from a subset of Feynman diagrams. Although these diagrams are not dominant at LO, their impact can be enhanced at NNLO due to the complicated two-loop effects. In contrast, in the \(gg\) channel, the relative difference between approximate and exact results is much smaller than what we expected, especially in the intermediate range \(\SI{2}{\TeV} \lesssim \sqrt{s_{12}} \lesssim \SI{5}{\TeV}\). This improvement originates from another accidental cancellation, which will become evident in the color decomposition of the squared amplitude defined in Eq. 23 in the following. As shown in Table [tab:ggnum1], the ratio of the approximate result in the massless scheme to the exact result for the squared amplitude as a whole is significantly smaller than the corresponding ratios for each individual color configuration in the \(\epsilon^{-1}\) column, owing to cancellations among different color coefficients.

We now turn to investigate the behavior of the approximate results as a function of the angle \(\theta_3\) between \(\vec{p}_3\) (the spatial momentum of the top quark) and \(\vec{p}_1\). The other parameters are fixed to \(\sqrt{s_{12}}= \SI{5}{\TeV}\), \(\phi_3=34\pi/29\), \(\theta_5=15\pi/29\) and \(q_5=20\,q_{5,\rm max}/29\). As expected, when \(\theta_3 \to 0\) or \(\pi\), the approximation becomes less reliable due to the smallness of either \(|s_{13}|\) or \(|s_{23}|\). Indeed, Fig. 2 shows that the approximate results deviate significantly from the exact ones at the two ends of the spectrum. Especially in the \(gg\) channel, the tree-level amplitudes contain propagators of the form \(1/(s_{13}-m_t^2)\) and \(1/(s_{23}-m_t^2)\), for which taking \(m_t \to 0\) introduces sizable differences. This explains the wild behavior of the orange and red curves in the \(gg\) channel near the two endpoints. Nevertheless, the massless scheme still provides a reasonable approximation across a large portion of the phase space in both the \(q\bar{q}\) and \(gg\) cases.

The situation is slightly different in the case of \(\theta_5\), the angle between \(\vec{p}_5\) (space configuration of the photon momentum) and \(\vec{p}_1\). Since the photon is color neutral, QCD corrections are expected to have only a mild dependence on its momentum, which is indeed observed in Fig. 3. Although collinear singularities arise at the two endpoints due to photon emission from the initial quark and antiquark (resulting in propagators \(1/s_{15}\) and \(1/s_{25}\)) in the \(q\bar{q}\) channel, the approximation in the massless scheme remains quite good across the whole range of \(\theta_5\), provided that \(\theta_3\) lies in the wide-angle region. By contrast, in the \(gg\) channel at NNLO the red curve displays weird behavior, indicating large numerical cancellations between the power-suppressed terms included in the semi-massive scheme and the unknown contributions from scalar integrals. Taken together with the previous observations, this indicates that the massless scheme performs more reliably overall than the semi-massive scheme. We therefore adopt the massless scheme as the default for presenting numerical results in what follows.

We now consider the NLO and NNLO squared amplitudes at various orders in \(\epsilon\) in the massless scheme. We define the following ratios of the approximate results to the exact ones at each order in \(\epsilon\): \[\begin{align} \bar{R}_{n,q/g}^{\text{NLO}}=\frac{\mathrm{Re}\braket{\bar{\mathcal{M}}^{(0)R}_{q/g}|\bar{\mathcal{M}}^{(1)R}_{q/g}}\bigg|_{\epsilon^n}}{\mathrm{Re}\braket{\mathcal{M}^{(0)R}_{q/g}|\mathcal{M}^{(1)R}_{q/g}}\bigg|_{\epsilon^n}}\,, \quad \bar{R}_{n,q/g}^{\text{NNLO}}=\frac{\mathrm{Re}\braket{\bar{\mathcal{M}}^{(0)R}_{q/g}|\bar{\mathcal{M}}^{(2)R}_{q/g}}\bigg|_{\epsilon^n}}{\mathrm{Re}\braket{\mathcal{M}^{(0)R}_{q/g}|\mathcal{M}^{(2)R}_{q/g}}\bigg|_{\epsilon^n}}\,, \end{align}\] where \(\ket{\bar{\mathcal{M}}^{(k)R}_{q/g}}\) are the amplitudes in the massless scheme, defined in Eq. 15 , while \(\ket{\mathcal{M}^{(k)R}_{q/g}}\) denote the exact massive amplitudes. \(\mathrm{Re}\braket{\bar{\mathcal{M}}^{(0)R}_{q/g}|\bar{\mathcal{M}}^{(k)R}_{q/g}}\) are evaluated using Eq. 19 in combination with Eq. 17 . The exact squared amplitudes \(\mathrm{Re}\braket{\mathcal{M}^{(0)R}_{q/g}|\mathcal{M}^{(k)R}_{q/g}}\) are known up to \(\epsilon^2\) at NLO and are computed up to \(\epsilon^{-1}\) at NNLO in this work.

We show \(\bar{R}_{n,q/g}^{\text{NLO}}\) and \(\bar{R}_{n,q/g}^{\text{NNLO}}\) as functions of \(\sqrt{s_{12}}\) in Fig. 4, using the same phase-space parameters as in Fig. 1. The approximation performs similarly across different coefficients in \(\epsilon\) for both channels, except for the \(\epsilon^{-1}\) term in the \(q\bar{q}\) channel at NNLO (already noted in Fig. 1), and for the leading \(1/\epsilon\) poles in the \(gg\) channel at both NLO and NNLO. As implied by Eq. 10 , the leading pole terms of the NLO and NNLO squared amplitudes are proportional to the LO ones as \(\epsilon \to 0\). We indeed observe an accidental cancellation at LP in the small-mass limit for the LO squared amplitudes in the \(gg\) channel at \(\mathcal{O}(\epsilon^0)\), while no such cancellation occurs at higher orders in \(\epsilon\). In Fig. 5, we show \(\bar{R}_{n,q,g}^{\text{NLO}}\) and \(\bar{R}_{n,q,g}^{\text{NNLO}}\) as functions of \(\theta_3\), using the same phase-space points as in Fig. 2. The approximation remains reliable across most values of \(\theta_3\) in the \(q\bar{q}\) channel, except near \(\theta_3 \to 0\), where the boosted limit is not valid due to small \(|s_{13}|\). In the \(gg\) channel, the approximation works very well in the central region but breaks down in the forward regions \(\theta_3 \to 0\) or \(\pi\), as expected.

The above results were obtained using the default scale choice \(\mu = m_t\). In Fig. 6, we show the dependence of the squared amplitudes on the renormalization scale at NLO and NNLO. The approximation remains stable across different values of \(\mu\), indicating that our formula correctly captures the scale dependence of the amplitudes.

We now present approximate predictions for the finite parts of the two-loop amplitudes for \(t\bar{t}\gamma\) production. Figure 7 shows the results as functions of \(\sqrt{s_{12}}\), \(\theta_3\) and \(\theta_5\), with phase-space points chosen as in Figs. 1, 2 and 3. Since the exact results for these finite terms are unknown, we provide approximate results from both the massless and semi-massive schemes. The two schemes agree well at high energies but begin to deviate as \(\sqrt{s_{12}}\) decreases. Based on our earlier findings, we conclude that the massless scheme yields a more reliable approximation at lower energies.

In addition to the squared amplitude as a whole, it is interesting to analyze its decomposition into color coefficients. The color structures of \(t\bar{t}\gamma\) production are identical to those of \(t\bar{t}\) and \(t\bar{t}H\) production. We therefore adopt the color decomposition of Refs. [37], [42], [73], [74]: \[\begin{align} \label{eq:color95dec} 2\,{\rm Re} \Braket{{\cal M}_{q}^{(0)}|{\cal M}_{q}^{(2)}} &= 2(N^2-1)\,\bigg( N^2 A^q + B^q + \frac{1}{N^2}\,C^q + N n_l\,D_l^q + N n_h\,D_h^q \nonumber \\ &+ \frac{n_l}{N}\,E_l^q + \frac{n_h}{N}\,E_h^q + n_l^2 F_l^q + n_l n_h\,F_{lh}^q + n_h^2 F_h^q \bigg) \,, \nonumber \\ 2\,{\rm Re} \Braket{{\cal M}_{g}^{(0)}|{\cal M}_{g}^{(2)}} &= (N^2-1)\,\bigg( N^3 A^g + N\,B^g + \frac{1}{N}\,C^g + \frac{1}{N^3}\,D^g \nonumber \\ &+ N^2 n_l\,E_l^g + N^2 n_h\,E_h^g + n_l\,F_l^g + n_h\,F_h^g + \frac{n_l}{N^2}\,G_l^g + \frac{n_h}{N^2}\,G_h^g \nonumber\\ &+ N n_l^2 H_l^g + N n_l n_h\,H_{lh}^g + N n_h^2 H_h^g + \frac{n_l^2}{N}\,I_l^g + \frac{n_l n_h}{N}\,I_{lh}^g + \frac{n_h^2}{N}\,I_h^g \bigg) \,. \end{align}\tag{23}\] In Tables [tab:qqnum1], [tab:ggnum1], [tab:qqnum2], [tab:ggnum2], [tab:qqnum3], and [tab:ggnum3], we list the numerical values of the color-decomposed two-loop squared amplitudes at different orders in \(\epsilon\). Results are shown at three representative phase-space points corresponding to relatively low, intermediate, and very high center-of-mass energies. We also perform a direct calculation of the \(n_f^2\) terms at NNLO, which contribute to the coefficients \(F^q\) in the \(q\bar{q}\) channel and \(H^g\), \(I^g\) in the \(gg\) channel. These results are also included in the tables. In the relatively low-energy region (Tables [tab:qqnum1] and [tab:ggnum1]), the relative differences between the massless-scheme results and the exact ones are significantly larger for \(F^q_h\) (\(H^g_h, I^g_h\)) than for \(F^q_l\) (\(H^g_l, I^g_l\)). In these cases, the closed quark loop effects factorize from the final-state top-quark pair, such that these coefficients share the same final-state contribution. Consequently, power corrections from the closed top-quark loop are more important than those from the external (anti-)top quark legs. This feature is interesting to explore further from the factorization perspective.¹⁰ A similar pattern is also observed for the terms proportional to \(n_f\), namely the coefficients \(D^q\) and \(E^q\) in the \(q\bar{q}\) channel, and \(E^g\), \(F^g\) and \(G^g\) in the \(gg\) channel. Overall, the massless scheme provides a reliable approximation to the exact results. These data may serve as a cross-check for a future exact evaluation of the two-loop squared amplitudes.

6 Conclusion↩︎

In this paper, we have calculated the two-loop amplitudes for \(t\bar{t}\gamma\) production at hadron colliders. First, for the infrared divergences, we employed the universal anomalous dimensions together with the one-loop massive amplitudes evaluated in dimensional regularization up to high order in \(\epsilon\) (see Eq. 10 ). Second, for the finite parts, we used the mass-factorization formula in the high-energy boosted limit (see Eq. 15 ), where scalar products of external momenta are much larger than the top-quark mass. In this framework, the massive amplitudes are expressed in terms of the corresponding massless ones multiplied by universal factors that capture collinear and soft dynamics. To this end, we generated the two-loop massive amplitudes to include part of the power corrections beyond LP, and computed the scalar integrals with vanishing parton masses using IBP reduction together with existing results for the massless master integrals. The massive amplitudes were then constructed in two schemes: the massless scheme, where all partons are treated as massless, and the semi-massive scheme, where only the scalar Feynman integrals are taken to be massless. In addition, we calculated exactly the \(n_f^2\) terms of the two-loop amplitudes.

We validated our approximations in the high-energy boosted limit from several perspectives. For the IR poles, we compared the predictions of both schemes with the exact results obtained in this work. For the finite remainders, we compared the approximate results with the exact ones for the \(n_f^2\) terms. We found that the massless scheme provides a better approximation across most of phase space, and we confirmed that our approach correctly reproduces the scale dependence of the exact results. Finally, we presented predictions for the finite parts of the two-loop squared amplitudes, which—together with the two-loop IR poles—constitute the main new results of this work. We also provided results for the squared amplitudes at three representative phase-space points, including their decomposition into color coefficients.

We find that power corrections in the \(q\bar{q}\) channel, especially in the relatively low-energy region, remain noticeable. Furthermore, the unknown power corrections from the scalar integrals are expected to significantly cancel the beyond-LP terms included in the semi-massive scheme. Improving the approximation therefore requires computing scalar integrals to higher powers in the high-energy limit, which can be achieved using strategies such as those described in Refs. [84]–[86]. Such results would also provide crucial input for studies of NLP factorization in the high-energy boosted limit, which we leave for future investigation. In addition, by combining our results with the real-emission contributions (double-real and real-virtual), one can compute the NNLO differential cross sections for \(t\bar{t}\gamma\) production in the high-energy boosted limit. Finally, by incorporating low-energy approximations such as threshold approximation and soft-photon approximation, and interpolating into the intermediate region, one would obtain a reasonable approximation across the whole phase space, enabling reliable predictions for differential cross sections from low to high energies.

This work was supported in part by the National Natural Science Foundation of China under Grant No. 12375097, 12535003, 12347103, and the ERC grant No. 101041109 (“BOSON”).

7 Renormalization constant↩︎

In this section, we present the renormalization constants needed in our work. Up to NNLO, the renormalization constant of QCD coupling constant \(\alpha_s\) is given by \[Z_{\alpha_s}=1-\left(\frac{\alpha_s}{4\pi}\right)\frac{\beta_0}{\epsilon}+\left(\frac{\alpha_s}{4\pi}\right)^2\left(\frac{\beta_0^2}{\epsilon^2}-\frac{\beta_1}{2\epsilon}\right) \,,\] where \[\beta_0=\frac{11}{3}C_A-\frac{4}{3}T_Fn_f,\quad \beta_1=\frac{34}{3}C_A^2-\frac{20}{3}C_AT_Fn_f-4C_FT_Fn_f.\] Up to NNLO, the top-quark mass renormalization constant in the on-shell scheme is given by: \[\begin{align} Z_m &= 1 + \frac{\alpha_s}{4\pi} C_F \left[ -\frac{3}{\epsilon} - \left(4+3L_\mu\right) - \epsilon\left(8+4L_\mu+\frac{\pi^2}{4}+\frac{3}{2}L_\mu^2\right) \right. \nonumber \\ & \left.+\epsilon^2\left(-L_{\mu}(8+\frac{\pi^2}{4}) -2L_{\mu}^2-\frac{L_{\mu}^3}{2} -48-\frac{\pi^2}{3}+\zeta_3 \right)\right] \nonumber \\ &+ \left(\frac{\alpha_s}{4\pi}\right)^2 \bigg\{C_A \left[\frac{11}{2 \epsilon ^2}-\frac{97}{12 \epsilon }+\frac{4 \pi ^2}{3}-\frac{1111}{24}+6 \zeta_3-\frac{11 L_{\mu }^2}{2} -4 \pi ^2 L_{\mu }-\frac{185 L_{\mu }}{6}\right]\nonumber \\ &+C_F \left[\frac{9}{2 \epsilon ^2}+\frac{36 L_{\mu }+45}{4\epsilon } +\frac{1}{8} \left(-96 \zeta_3+199-34 \pi ^2\right) + 9 L_{\mu }^2+8 \pi ^2 L_{\mu } +\frac{45 L_{\mu }}{2}\right]\nonumber \\ & + T_Fn_h \bigg(-\frac{2}{\epsilon ^2}+\frac{5}{3 \epsilon } +2 L_{\mu }^2+\frac{26 L_{\mu }}{3}-\frac{8 \pi ^2}{3}+\frac{143}{6}\bigg)\nonumber \\ &+T_Fn_l \bigg[-\frac{2}{\epsilon ^2}+\frac{5}{3 \epsilon } +\frac{1}{6} \left(71+8 \pi ^2\right)+2 L_{\mu }^2+\frac{26 L_{\mu }}{3}\bigg]\bigg\}+ \mathcal{O}(\alpha_s^3) \, , \end{align}\] where \(L_{\mu} = \ln\left(\mu^2/m_t^2\right)\). The on-shell wave-function renormalization constants are \[\begin{align} Z_{q} &= 1 + \left(\frac{\alpha_s}{4\pi}\right)^2 C_F \left(\frac{1}{2\epsilon} -\frac{5}{12} + L_{\mu}\right) + \mathcal{O}(\alpha_s^3) \, , \nonumber \\ Z_{g} &= 1 + \frac{\alpha_s}{4\pi} T_Fn_h \left[ -\frac{4}{3\epsilon} - \frac{4}{3}L_{\mu} - \epsilon\left(\frac{\pi^2}{9} + \frac{2}{3}L^2_{\mu} \right)+\epsilon^2\left( -\frac{\pi^2}{9} L_{\mu}-\frac{2}{9} L_{\mu}^3 +\frac{4}{9}\zeta_3\right)\right] \nonumber \\ &+ \left(\frac{\alpha_s}{4\pi}\right)^2 T_Fn_h \left\{C_A \left[\frac{35}{9 \epsilon ^2}-\frac{1}{\epsilon}\left(\frac{5}{2} - \frac{26 L_{\mu }}{9}\right)+\frac{13}{12}-5 L_{\mu }+\frac{13 \pi ^2}{54}+\frac{4 L_{\mu }^2}{9}\right] \right. \nonumber \\ &\left.-C_F \left(\frac{2}{\epsilon }+15+4 L_{\mu}\right)+n_h \left(\frac{8 L_{\mu }}{9 \epsilon}+\frac{2 \pi ^2}{27}+ \frac{4 L_{\mu }^2}{3}\right) \right. \nonumber \\ &\left.-n_l \left(\frac{8}{9 \epsilon ^2}+\frac{8 L_{\mu }}{9 \epsilon}+\frac{2 \pi ^2}{27}+\frac{4 L_{\mu }^2}{9}\right)\right\}+ \mathcal{O}(\alpha_s^3) \, , \nonumber \\ Z_{Q} &= 1 + \frac{\alpha_s}{4\pi} C_F \left[ -\frac{3}{\epsilon} - \left(4+3L_\mu\right) - \epsilon\left(8+4L_\mu+\frac{\pi^2}{4}+\frac{3}{2}L_\mu^2\right) \right. \nonumber \\ & \qquad\left.+\epsilon^2\left(-L_{\mu}(8+\frac{\pi^2}{4}) -2L_{\mu}^2-\frac{L_{\mu}^3}{2} -48-\frac{\pi^2}{3}+\zeta_3 \right) \right] \nonumber \\ & + \left(\frac{\alpha_s}{4\pi}\right)^2 C_F\left\{C_A \left(\frac{11}{2 \epsilon^2}-\frac{127}{12 \epsilon }-\frac{1705}{24}-\frac{215 L_{\mu }}{6}+5 \pi ^2-\frac{11 L_{\mu }^2}{2}-8 \pi ^2 \ln(2) + 12 \zeta _3\right) \right.\nonumber \\ & \left. +C_F \left[\frac{9}{2 \epsilon^2}+\frac{1}{\epsilon}\left(\frac{51}{4 }+9 L_{\mu }\right)+\frac{433}{8}+\frac{51 L_{\mu }}{2}-\frac{49 \pi ^2}{4}+9 L_{\mu }^2+16 \pi ^2 \log (2)-24\zeta _3\right]\right.\nonumber \\ & \left.+n_h\left[\frac{1}{2\epsilon}\left(1 + 4 L_{\mu}\right)+\frac{947}{36}+\frac{11 L_{\mu }}{3}-\frac{5 \pi ^2}{2}+3 L_{\mu }^2\right]\right.\nonumber \\ & \left.+n_l \left(-\frac{1}{\epsilon ^2}+\frac{11}{6 \epsilon }+\frac{113}{12}+\frac{19 L_{\mu}}{3}+\frac{2 \pi ^2}{3}+L_{\mu }^2\right)\right\}+ \mathcal{O}(\alpha_s^3) \, . \end{align}\]

References↩︎