# Heavy-quark-pair production at lepton colliders at NNNLO in QCD

Xiang Chen<EMAIL_ADDRESS>School of Physics, Peking University, Beijing
100871, China  Xin Guan<EMAIL_ADDRESS>School of Physics, Peking
University, Beijing 100871, China  Chuan-Qi He<EMAIL_ADDRESS>School of
Physics, Peking University, Beijing 100871, China  Xiao Liu
<EMAIL_ADDRESS>Rudolf Peierls Centre for Theoretical Physics,
Clarendon Laboratory, Parks Road, Oxford OX1 3PU, UK  Yan-Qing Ma
<EMAIL_ADDRESS>School of Physics, Peking University, Beijing 100871, China
Center for High Energy Physics, Peking University, Beijing 100871, China

###### Abstract

We compute the total cross section of heavy-quark-pair production in
$e^{+}e^{-}$ annihilation at the next-to-next-to-next-to-leading order (NNNLO)
in Quantum Chromodynamics. The result is expressed as a piecewise function
defined by several deeply expanded power series. The result significantly
reduces the theoretical uncertainty. For a typical collision energy of 500
GeV, the scale dependence has been reduced from $0.94\%$ at the next-to-next-
to-leading order (NNLO) to $0.11\%$ at the NNNLO, which meets the request by
future lepton colliders.

Introduction. —  As the heaviest particle in the Standard Model of particle
physics, the top quark plays a crucial role in both precision Electroweak (EW)
physics and physics beyond the Standard Model [1]. Currently, the top quark is
mainly studied by hadron colliders, which, for example, determine top-quark
mass with uncertainties about 500-600 MeV [2, 3].

Study of the top-quark pair production at future lepton colliders is important
to further pin down properties of the top quark [4, 5, 6, 7, 8, 9, 10], due to
the substantially cleaner hadronic environment aspects related to hadronic
initial state radiation. For example, at the International Linear Collider
[10], uncertainties of top-quark mass can be reduced to 50 MeV by measuring
top-quark pair production near threshold. While by measuring the cross section
and forward-backward asymmetry for the top-quark pair production at 500 GeV,
form factors for the top quark couplings to the photon and the Z boson can be
determined with relative uncertainties smaller than 0.3%. To exploit the full
potential of future colliders, significant theoretical effort is required,
where one of the most important things is the computation of cross section of
top-pair production at higher orders in perturbation theory.

Full next-to-leading order (NLO) Chromodynamics (QCD) correction was first
computed in Ref. [11] and NLO EW effects together with NLO QCD correction was
provided in Ref. [12]. Next-to-next-to-leading order (NNLO) QCD correction to
the total cross section was first obtained in Ref. [13] by
Pad$\acute{\text{e}}$ approximation based on the results of threshold
expansions, high-energy expansions and low-energy expansions. Direct
calculation of the cross section and differential distributions at NNLO QCD
with full top-mass dependence was provided in Refs. [14, 15, 16]. Recently, it
was re-computed purely numerically using the Local-Unitarity method [17]. At
the next-to-next-to-next-to-leading order (NNNLO) level, an exploration study
of the total cross section by Pad$\acute{\text{e}}$ approximation can be found
in Ref. [18] which unfortunately cannot provide an reliable result in physical
region. For top-pair production cross section near threshold, NNNLO QCD
correction has been achieved in Refs. [19, 20]. In Refs. [21, 22], massive
form factors at three loops was obtained, which are building blocks for
complete NNNLO QCD corrections.

In this Letter, we present a full computation of NNNLO QCD correction to the
total cross section of $e^{+}e^{-}\to\gamma^{*}/Z^{*}\to t\bar{t}$ process
mediated by an off-shell photon ($\gamma^{*}$) and/or an off-shell $Z$ boson
($Z^{*}$), valid in the whole physical region. Our computation is achieved in
a semi-numerical way taking advantage of newly proposed methods and released
tools. All integrals including phase-space integrals and loop integrals are
expressed as piecewise functions represented by several deeply expanded power
series in the physical region. As a result, the final cross section can be
used precisely and efficiently. As a by product, we also obtain NNNLO QCD
corrections for the production cross sections of $c\bar{c}$ and $b\bar{b}$.

Calculation. —  We generate the Feynman amplitudes with qgraf [23], with some
sample Feynman diagrams shown in Fig 1. We use our in-house Mathematica
package to deal with Lorentz and color algebras, and express all squared
amplitudes as linear combinations of scalar integrals, which have been mapped
to predefined integral families respectively. Phase-space integrals are
treated in the same manner as loop integrals in our calculation taking
advantage of the reverse unitarity [24, 25, 26]. The coefficients of integrals
are some simple rational functions of the squared center-of-mass energy $s$,
squared heavy quark mass $m^{2}$, dimensional regulator $\epsilon=(4-D)/2$,
number of light fermions $n_{l}$ and $N_{c}$ of the $SU(N_{c})$ gauge group.

(a)

(b)

(c)

(d)

Figure 1:  Some sample Feynman diagrams for heavy quark pair production at
NNNLO in QCD.

Integration-by-parts (IBP) reduction [27] is used to express integrals in each
family as linear combinations of a minimal set of so-called master integrals.
We realize this by a two-step strategy. At the first step, we use LiteRed [28]
and FiniteFlow [29] to generate and solve the system of IBP identities based
on Laporta’s algorithm [30] over finite field. A small number (typically 50)
of numerical samplings are sufficient for us to construct the block-triangular
relations for target integrals proposed in Ref. [31, 32]. At the second step,
we use these much more efficient block-triangular relations to obtain large
amount (typically 1000) of samplings to eventually reconstruct the reduction
coefficients. In this work, the two-step strategy typically reduces the
computational time by one to two orders, depending on families, compared with
the reduction without using block-triangular relations. We note that this is
the first time the block-triangular relations are used to tackle a physical
problem.

The master integrals are computed using differential equations [33] based on
power series expansion [34, 35]. The differential equations of master
integrals with respect to $x=4m^{2}/s$ are constructed using aforementioned
IBP reduction. To fix the boundary conditions, say at an arbitrary regular
point $x=4/23$, we utilize the auxiliary mass flow method [36, 37, 38, 39]
implemented in AMFlow [40]. With these in hand, we are able to construct a
piecewise function for each master integral using the differential equation
solver in AMFlow, as 5 deeply expanded power series at the following points

$\displaystyle x\to\left\\{0,\frac{1}{4},\frac{1}{2},\frac{3}{4},1\right\\},$
(1)

where $0$, $1/4$ and $1$ are singular points and $1/2$ and $3/4$ are regular
points. The expansion at each $x_{0}$ in Eq. (1) is valid at least in
$x\in(x_{0}-1/8,x_{0}+1/8)$, thus the physical region $x\in(0,1)$ is fully
covered.

It is well known that, in the dimensional regularization, the anticommutation
relation $\\{\gamma^{\mu},\gamma_{5}\\}=0$ and the cyclicity of Dirac trace
cannot be satisfied simultaneously. In practice, we prefer to keep the
anticommutation relation $\\{\gamma^{\mu},\gamma_{5}\\}=0$ not only to
simplify the computation, but also because using a non-anticommuting
$\gamma_{5}$ with the other Dirac matrices $\gamma^{\mu}$ leads to “spurious
anomalies” which violate chiral symmetry and hence gauge invariance [41]. In
this case, one must impose the relevant Ward-Takahashi and Slavnov-Taylor
identities order by order to keep the renormalizability of a perturbative
theory, which will make practical calculations much more different and tedious
obviously. Thus we choose the so-called naïve-$\gamma_{5}$ scheme so that
anticommutation relations hold. Then a new problem arises: where should we
read a fermion loop? Due to the lack of the cyclicity of the trace, reading a
fermion loop from a different point may give different final results. This
ambiguity can be fixed by the KKS-scheme [42, 43, 44], where the final result
is defined as the average of all possible reading points starting from the
head and the tail of an axial-vector three-point fermion-gauge subgraph that
contains the maximal one-particle-irreducible non-singlet-type loop correction
[44, 45].

As we are working in the bare perturbation theory, we need to replace bare
quantities by renormalized ones before physical results can be obtained. QCD
coupling is renormalized in the $\overline{\text{MS}}$ scheme, and heavy quark
mass and fields are renormalized in the on-shell scheme. Both coupling
renormalization and field renormalization can be realized by simply
multiplying the bare results by proper powers of
$Z_{\alpha_{s}}^{\overline{\text{MS}}}$ and $Z_{2}^{\text{OS}}$. Heavy quark
mass renormalization is more complicated because it is not multiplicatively
renormalizable. We replace bare mass by renormalized mass via
$m^{b}=m^{\text{OS}}+\delta m$, and keep on-shell condition of external heavy
quark momentum $p^{2}=(m^{\text{OS}})^{2}$ untouched. In spirit of
perturbation theory, $\delta m$ is a small quantity and can be expanded to any
desired order. Especially, the expansion for master integrals is achieved by
setting up and solving differential equations with respect to $\delta m$.
After the expansion of $\delta m$, we can substitute renormalization
constants, which can be found in Ref. [46, 47, 48], to obtain final
renormalized cross section. We note that performing renormalization in this
way makes our calculation very systematic, and at the same time introduces
negligible efforts against the calculation of bare cross section at the
highest order in $\alpha_{s}$.

Another key technique in our calculation is that the master integrals are
calculated with numerical values of $\epsilon$ and the $\epsilon$ dependence
is only reconstructed at the cross section level, as proposed in Ref. [39,
40]. The advantage of this technique is that we do not need to manipulate
Laurent expansions of $\epsilon$ during the intermediate stages of
calculations, which significantly reduces the computational time.

Results. —  By combining everything together, we end up with final result at
the NNNLO level, which is free of ultraviolet and infrared divergences as
expected from the Kinoshita-Lee-Nauenberg theorem [49, 50]. Our result of
NNNLO cross section is expressed as a piecewise function of $x=4m^{2}/s$
represented by 5 power series, all expanded to 40 orders, allowing us to
obtain at least 10 correct digits in physical region. The expressions can be
found in a computer-readable ancillary file attached to this Letter.

To provide numerical result for top-quark-pair production, we choose top-quark
mass as $m=173.34$ GeV [51], and set all other quarks as massless.
Electromagnetic coupling is chosen as a fixed value $\alpha=1/132.2$, and
strong coupling $\alpha_{s}(\mu)$ is running as a function of renormalization
scale $\mu$, which is computed using the RunDec package [52, 53] with input
value $\alpha_{s}^{n_{f}=5}(m_{Z}=91.1876\text{ GeV})=0.1181$.

Figure 2:  Total cross section for $t\bar{t}$ production.

Fig. 2 shows our result for the total cross section of $t\bar{t}$ production,
where LO, NLO, NNLO and NNNLO are shown in black, red, blue and brown
respectively. The middle lines of each band correspond to the choice of
$\mu=\sqrt{s}$ for the renormalization scale, and the upper and lower lines
correspond to the scale variations between $\mu=\sqrt{s}/2$ and
$\mu=2\sqrt{s}$. It can be found that the NNNLO correction significantly
reduces the scale dependence. The dashed, dotted, and dot-dashed lines
correspond to $\gamma-\gamma$, $Z-Z$, $\gamma-Z$ and $Z-\gamma$ contributions
at NNNLO, respectively, since the intermediate off-shell boson of the
amplitude and complex conjugate of which could be a photon ($\gamma$) or $Z$
boson ($Z$).

Figure 3:  $\Delta_{1}$, $\Delta_{2}$ and $\Delta_{3}$ as functions of the
center-of-mass energy.  Figure 4:  Scale dependence of $\Delta_{1}$,
$\Delta_{1}+\Delta_{2}$ and $\Delta_{1}+\Delta_{2}+\Delta_{3}$ for $t\bar{t}$
production at $\sqrt{s}=500$ GeV.

The total cross section for $t\bar{t}$ production can also be expressed in the
form

$\sigma_{\text{NNNLO}}=\sigma_{\text{LO}}(1+\Delta_{1}+\Delta_{2}+\Delta_{3}),$
(2)

where the order $\alpha_{s}$, $\alpha_{s}^{2}$ and $\alpha_{s}^{3}$
corrections $\Delta_{1}$, $\Delta_{2}$ and $\Delta_{3}$ are displayed in Fig.
3 as functions of the center-of-mass energy, respectively. The smallness of
$\Delta_{3}$ confirms a good convergence of perturbative expansion with
respect to $\alpha_{s}$. In Fig. 4, we further show the reduction of the scale
dependence after including $O(\alpha_{s}^{3})$ correction. We vary the
renormalization scale between $\mu=0.3\sqrt{s}$ and $\mu=3\sqrt{s}$. It turns
out that the scale dependence has been reduced from $0.94\%$ at NNLO to
$0.11\%$ at NNNLO for a collision energy of 500 GeV, which meets the precision
requested by, e.g., International Linear Collider [10]. However, the NNNLO
result still shows a large scale dependence near the top-quark pair production
threshold due to Coulomb interaction and requires resummation for further
improvement [54, 19, 20].

Figure 5:  Total cross section for $b\bar{b}$ production.  Figure 6:  Total
cross section for $c\bar{c}$ production.

Our results for $b\bar{b}$ and $c\bar{c}$ total cross section are shown in
Figs. 5 and 6, where on-shell heavy-quark mass is chosen as $m_{b}=4.8$ GeV
and $m_{c}=1.5$ GeV, respectively. We find that theoretical uncertainties are
also reduced, and thus they are useful to describe corresponding processes
when experimental precision is improved. Note that here we only consider the
bottom-quark and charm-quark pair production mediated by an off-shell photon
($\gamma^{*}$) because the mass of $Z$ boson is much larger than the energy
scale where we are interested in. Also, we only take account of non-singlet
contributions in the bottom and charm quark case because the singlet
contributions are very small, which we can see from the calculation of top-
quark pair production. For example, at $\sqrt{s}=500$ GeV, the singlet
contributions amount to $-5.8\%$ of the total third order correction
$\Delta_{3}$ defined in Eq 2. Here we define “singlet” and “non-singlet” at
the level of Feynman diagrams which correspond to that the external current
does not couple/couples directly to the fermion line of the final-state heavy
quarks.

Verification. —  Besides free of ultraviolet and infrared divergences, our
results have been checked in several other aspects, which we describe in the
following.

Our result for $t\bar{t}$ production at NNLO agrees with that obtained in
Refs. [14, 16]. The three-loop contribution of our NNNLO result agrees with
the latest results for massive form factors [21, 22]. All master integrals
involved in our calculation, obtained by solving differential equations with
respect to $x$ with boundary conditions given at $x=4/23$, have been checked
by AMFlow at another phase space point $x=15/16$.

The most rigorous check is that we calculate the cross section in another
totally different way and get a perfect agreement. We calculate the forward
scattering amplitude of a virtual photon at 4-loop level without cut, and then
relate its imaginary part to the desired cross section by virtue of optical
theorem. Note that using optical theorem could induce massless final states
without any heavy quark, thus we must subtract these contributions from the
final results. Similarly, the contribution from four-heavy-quark production
needs to be subtracted when we consider the region above its threshold, i.e.
$x<1/4$.

Additionally, to ensure the appropriateness of the choice of $\gamma_{5}$
scheme, we conducted a cross check by calculating the singlet part of axial-
axial contribution at NNNLO using the so-called Larin scheme [55, 56]. Our
focus primarily lies on the axial-axial contribution since there is no
$\gamma_{5}$ in vector-vector contribution and the vector-axial contribution
is zero. Furthermore, for the non-singlet part of axial-axial contribution, it
is unambiguous using anticommuting $\gamma_{5}$ scheme because there will
always be even number of $\gamma_{5}$ in the Dirac trace which allows them to
contract without any issues, taking advantage of the anticommutation relation
$\\{\gamma^{\mu},\gamma_{5}\\}=0$ and $\gamma_{5}^{2}=1$. Hence we only need
to consider the singlet part of axial-axial contribution. We found that the
results obtained using the KKS scheme agree with those obtained using the
Larin scheme, after finite renormalization was conducted in the later scheme.

Summary. —  In this Letter we provide the first complete NNNLO perturbative
correction to heavy-quark-pair production cross section at lepton colliders.
Our results are obtained by combining several powerful techniques, such as
block-triangular relations and numerical sampling of $\epsilon$. The
calculation are done in two totally different ways for cross check: one way is
to calculate virtual corrections and real corrections separately, and the
other way is based on optical theorem to calculate forward scattering
amplitude. Our result significantly reduces theoretical uncertainty and meet
the request by the future lepton colliders. Besides, our strategy is
applicable for many other processes, which will be helpful to push the
corresponding perturbative calculation to an unprecedented order.

###### Acknowledgements.

Acknowledgments. —  We thank L. Chen for many useful communications and
discussions. The work was supported in part by the National Natural Science
Foundation of China (Grants No. 11875071, No. 11975029), the National Key
Research and Development Program of China under Contracts No. 2020YFA0406400,
and the High-performance Computing Platform of Peking University. The research
of XL was also supported by the ERC Starting Grant 804394 HipQCD and by the UK
Science and Technology Facilities Council (STFC) under grant ST/T000864/1.
JaxoDraw [57] was used to generate Feynman diagrams.

## References

  * [1] K. Agashe et al., Report of the Topical Group on Top quark physics and heavy flavor production for Snowmass 2021, [arXiv:2209.11267] [InSPIRE].
  * [2] ATLAS , M. Aaboud et al., Measurement of the top quark mass in the $t\bar{t}\rightarrow$ lepton+jets channel from $\sqrt{s}=8$ TeV ATLAS data and combination with previous results, Eur. Phys. J. C 79 (2019) 290 [arXiv:1810.01772] [InSPIRE].
  * [3] CMS , V. Khachatryan et al., Measurement of the top quark mass using proton-proton data at ${\sqrt{(s)}}$ = 7 and 8 TeV, Phys. Rev. D 93 (2016) 072004 [arXiv:1509.04044] [InSPIRE].
  * [4] The International Linear Collider Technical Design Report - Volume 2: Physics, [arXiv:1306.6352] [InSPIRE].
  * [5] P. Bambade et al., The International Linear Collider: A Global Project, [arXiv:1903.01629] [InSPIRE].
  * [6] FCC , A. Abada et al., FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2, Eur. Phys. J. ST 228 (2019) 261–623 [InSPIRE].
  * [7] CEPC Study Group , M. Dong et al., CEPC Conceptual Design Report: Volume 2 - Physics & Detector, [arXiv:1811.10545] [InSPIRE].
  * [8] M. Vos et al., Top physics at high-energy lepton colliders, [arXiv:1604.08122] [InSPIRE].
  * [9] CLICdp , H. Abramowicz et al., Top-Quark Physics at the CLIC Electron-Positron Linear Collider, JHEP 11 (2019) 003 [arXiv:1807.02441] [InSPIRE].
  * [10] ILD Concept Group , H. Abramowicz et al., International Large Detector: Interim Design Report, [arXiv:2003.01116] [InSPIRE].
  * [11] J. Jersak, E. Laermann, and P. M. Zerwas, Electroweak Production of Heavy Quarks in e+ e- Annihilation, Phys. Rev. D 25 (1982) 1218 [InSPIRE]. [Erratum: Phys.Rev.D 36, 310 (1987)].
  * [12] W. Beenakker, S. C. van der Marck, and W. Hollik, e+ e- annihilation into heavy fermion pairs at high-energy colliders, Nucl. Phys. B 365 (1991) 24–78 [InSPIRE].
  * [13] K. G. Chetyrkin, J. H. Kuhn, and M. Steinhauser, Three loop polarization function and O (alpha-s**2) corrections to the production of heavy quarks, Nucl. Phys. B 482 (1996) 213–240 [hep-ph/9606230] [InSPIRE].
  * [14] J. Gao and H. X. Zhu, Electroweak prodution of top-quark pairs in e+e- annihilation at NNLO in QCD: the vector contributions, Phys. Rev. D 90 (2014) 114022 [arXiv:1408.5150] [InSPIRE].
  * [15] J. Gao and H. X. Zhu, Top Quark Forward-Backward Asymmetry in $e^{+}e^{-}$ Annihilation at Next-to-Next-to-Leading Order in QCD, Phys. Rev. Lett. 113 (2014) 262001 [arXiv:1410.3165] [InSPIRE].
  * [16] L. Chen, O. Dekkers, D. Heisler, W. Bernreuther, and Z.-G. Si, Top-quark pair production at next-to-next-to-leading order QCD in electron positron collisions, JHEP 12 (2016) 098 [arXiv:1610.07897] [InSPIRE].
  * [17] Z. Capatti, V. Hirschi, and B. Ruijl, Local Unitarity: cutting raised propagators and localising renormalisation, [arXiv:2203.11038] [InSPIRE].
  * [18] Y. Kiyo, A. Maier, P. Maierhofer, and P. Marquard, Reconstruction of heavy quark current correlators at O(alpha(s)**3), Nucl. Phys. B 823 (2009) 269–287 [arXiv:0907.2120] [InSPIRE].
  * [19] M. Beneke, Y. Kiyo, P. Marquard, A. Penin, J. Piclum, and M. Steinhauser, Next-to-Next-to-Next-to-Leading Order QCD Prediction for the Top Antitop $S$-Wave Pair Production Cross Section Near Threshold in $e^{+}e^{-}$ Annihilation, Phys. Rev. Lett. 115 (2015) 192001 [arXiv:1506.06864] [InSPIRE].
  * [20] M. Beneke, Y. Kiyo, A. Maier, and J. Piclum, Near-threshold production of heavy quarks with $\tt{QQbar\\_threshold}$, Comput. Phys. Commun. 209 (2016) 96–115 [arXiv:1605.03010] [InSPIRE].
  * [21] M. Fael, F. Lange, K. Schönwald, and M. Steinhauser, Massive Vector Form Factors to Three Loops, Phys. Rev. Lett. 128 (2022) 172003 [arXiv:2202.05276] [InSPIRE].
  * [22] M. Fael, F. Lange, K. Schönwald, and M. Steinhauser, Singlet and nonsinglet three-loop massive form factors, Phys. Rev. D 106 (2022) 034029 [arXiv:2207.00027] [InSPIRE].
  * [23] P. Nogueira, Automatic Feynman graph generation, J. Comput. Phys. 105 (1993) 279–289 [InSPIRE].
  * [24] C. Anastasiou and K. Melnikov, Higgs boson production at hadron colliders in NNLO QCD, Nucl. Phys. B 646 (2002) 220–256 [hep-ph/0207004] [InSPIRE].
  * [25] C. Anastasiou, L. J. Dixon, and K. Melnikov, NLO Higgs boson rapidity distributions at hadron colliders, Nucl. Phys. B Proc. Suppl. 116 (2003) 193–197 [hep-ph/0211141] [InSPIRE].
  * [26] C. Anastasiou, L. J. Dixon, K. Melnikov, and F. Petriello, Dilepton rapidity distribution in the Drell-Yan process at NNLO in QCD, Phys. Rev. Lett. 91 (2003) 182002 [hep-ph/0306192] [InSPIRE].
  * [27] K. G. Chetyrkin and F. V. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl. Phys. B192 (1981) 159–204 [InSPIRE].
  * [28] R. N. Lee, LiteRed 1.4: a powerful tool for reduction of multiloop integrals, J. Phys. Conf. Ser. 523 (2014) 012059 [arXiv:1310.1145] [InSPIRE].
  * [29] T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs, JHEP 07 (2019) 031 [arXiv:1905.08019] [InSPIRE].
  * [30] S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations, Int. J. Mod. Phys. A 15 (2000) 5087–5159 [hep-ph/0102033] [InSPIRE].
  * [31] X. Liu and Y.-Q. Ma, Determining arbitrary Feynman integrals by vacuum integrals, Phys. Rev. D 99 (2019) 071501 [arXiv:1801.10523] [InSPIRE].
  * [32] X. Guan, X. Liu, and Y.-Q. Ma, Complete reduction of two-loop five-light-parton scattering amplitudes, Chin. Phys. C 44 (2020) 9 [arXiv:1912.09294] [InSPIRE].
  * [33] A. V. Kotikov, Differential equations method: New technique for massive Feynman diagrams calculation, Phys. Lett. B254 (1991) 158–164 [InSPIRE].
  * [34] M. Caffo, H. Czyz, M. Gunia, and E. Remiddi, BOKASUN: A Fast and precise numerical program to calculate the Master Integrals of the two-loop sunrise diagrams, Comput. Phys. Commun. 180 (2009) 427–430 [arXiv:0807.1959] [InSPIRE].
  * [35] M. Czakon, Tops from Light Quarks: Full Mass Dependence at Two-Loops in QCD, Phys. Lett. B664 (2008) 307–314 [arXiv:0803.1400] [InSPIRE].
  * [36] X. Liu, Y.-Q. Ma, and C.-Y. Wang, A Systematic and Efficient Method to Compute Multi-loop Master Integrals, Phys. Lett. B779 (2018) 353–357 [arXiv:1711.09572] [InSPIRE].
  * [37] X. Liu, Y.-Q. Ma, W. Tao, and P. Zhang, Calculation of Feynman loop integration and phase-space integration via auxiliary mass flow, Chin. Phys. C 45 (2021) 013115 [arXiv:2009.07987] [InSPIRE].
  * [38] X. Liu and Y.-Q. Ma, Multiloop corrections for collider processes using auxiliary mass flow, Phys. Rev. D 105 (2022) L051503 [arXiv:2107.01864] [InSPIRE].
  * [39] Z.-F. Liu and Y.-Q. Ma, Feynman integrals are completely determined by linear algebra, [arXiv:2201.11637] [InSPIRE].
  * [40] X. Liu and Y.-Q. Ma, AMFlow: a Mathematica package for Feynman integrals computation via Auxiliary Mass Flow, [arXiv:2201.11669] [InSPIRE].
  * [41] F. Jegerlehner, Facts of life with gamma(5), Eur. Phys. J. C 18 (2001) 673–679 [hep-th/0005255] [InSPIRE].
  * [42] D. Kreimer, The $\gamma$(5) Problem and Anomalies: A Clifford Algebra Approach, Phys. Lett. B 237 (1990) 59–62 [InSPIRE].
  * [43] J. G. Korner, D. Kreimer, and K. Schilcher, A Practicable gamma(5) scheme in dimensional regularization, Z. Phys. C 54 (1992) 503–512 [InSPIRE].
  * [44] D. Kreimer, The Role of gamma(5) in dimensional regularization, [hep-ph/9401354] [InSPIRE].
  * [45] L. Chen, An observation on Feynman diagrams with axial anomalous subgraphs in dimensional regularization with an anticommuting $\gamma_{5}$, [arXiv:2304.13814] [InSPIRE].
  * [46] K. Melnikov and T. van Ritbergen, The Three loop on-shell renormalization of QCD and QED, Nucl. Phys. B 591 (2000) 515–546 [hep-ph/0005131] [InSPIRE].
  * [47] M. Czakon, A. Mitov, and S. Moch, Heavy-quark production in gluon fusion at two loops in QCD, Nucl. Phys. B 798 (2008) 210–250 [arXiv:0707.4139] [InSPIRE].
  * [48] P. Bärnreuther, M. Czakon, and P. Fiedler, Virtual amplitudes and threshold behaviour of hadronic top-quark pair-production cross sections, JHEP 02 (2014) 078 [arXiv:1312.6279] [InSPIRE].
  * [49] T. Kinoshita, Mass singularities of Feynman amplitudes, J.Math.Phys. 3 (1962) 650–677 [InSPIRE].
  * [50] T. Lee and M. Nauenberg, Degenerate Systems and Mass Singularities, Phys.Rev. 133 (1964) B1549–B1562 [InSPIRE].
  * [51] Particle Data Group , R. L. Workman et al., Review of Particle Physics, PTEP 2022 (2022) 083C01 [InSPIRE].
  * [52] K. G. Chetyrkin, J. H. Kuhn, and M. Steinhauser, RunDec: A Mathematica package for running and decoupling of the strong coupling and quark masses, Comput. Phys. Commun. 133 (2000) 43–65 [hep-ph/0004189] [InSPIRE].
  * [53] F. Herren and M. Steinhauser, Version 3 of RunDec and CRunDec, Comput. Phys. Commun. 224 (2018) 333–345 [arXiv:1703.03751] [InSPIRE].
  * [54] A. H. Hoang and M. Stahlhofen, The Top-Antitop Threshold at the ILC: NNLL QCD Uncertainties, JHEP 05 (2014) 121 [arXiv:1309.6323] [InSPIRE].
  * [55] S. A. Larin and J. A. M. Vermaseren, The alpha-s**3 corrections to the Bjorken sum rule for polarized electroproduction and to the Gross-Llewellyn Smith sum rule, Phys. Lett. B 259 (1991) 345–352 [InSPIRE].
  * [56] S. A. Larin, The Renormalization of the axial anomaly in dimensional regularization, Phys. Lett. B 303 (1993) 113–118 [hep-ph/9302240] [InSPIRE].
  * [57] D. Binosi and L. Theussl, Jaxodraw: A graphical user interface for drawing feynman diagrams, Computer Physics Communications 161 (2004) 76–86 [InSPIRE]. https://www.sciencedirect.com/science/article/pii/S0010465504002115.