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

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 deﬁned by several deeply expanded power series. The result sig-niﬁcantly 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 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, it is imperative to undertake advanced computations for the cross-section and differential distributions of top-pair production at higher orders within the framework of perturbation theory.
Full next-to-leading order (NLO) QCD correction was first computed in Ref. [11] and NLO electroweak 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] using Padé 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][17].At the next-to-next-to-next-to-leading order (NNNLO) level, partial total cross-section was ob-tained utilizing Padé approximation [18,19], specifically excluding the singlet contribution like diagrams such as Fig. 1 (b).For top-pair production cross-section near threshold, NNNLO QCD correction has been achieved in Refs.[20,21].In Refs.[22,23], massive form factors at three loops was obtained, which are building blocks for complete NNNLO QCD corrections.Despite significant theoretical advancements, a comprehensive NNNLO QCD calculation for the total cross-section is notably absent.Moreover, tackling the complexities of differential observables remains a substantial challenge, with no single study addressing this in the existing literature.
In this Letter, we offer a comprehensive computation of the NNNLO QCD corrections to the total cross-section in the process e + e − → γ * /Z * → t t.This encompasses scenarios involving an off-shell photon (γ * ) and/or an off-shell Z boson (Z * ), and the validity extends across the entire physical region.Additionally, we introduce a first computation of the differential distribution up to NNNLO in QCD, exemplified by the invariant mass distribution of the top-quark pair.Our computations are achieved in a semi-analytical 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 total cross-section and invariant mass distribution can be used precisely and efficiently.Employing identical methodologies, we further derive NNNLO QCD corrections for the production cross-sections of cc and b b.
Calculation.-We take the total cross-section as an example to illustrate our computation methods.We generate the Feynman amplitudes with qgraf [24], with some sample Feynman diagrams shown in Fig. 1.Note that the four-top final state contributions, which can be arXiv:2209.14259v3[hep-ph] 11 Mar 2024 measured separately, are not included in our calculation.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 [25][26][27].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 ϵ = (4 − D)/2, number of light fermions n l and N c of the SU (N c ) gauge group.
FIG. 1: Representative Feynman diagrams at NNNLO.Integration-by-parts (IBP) reduction [28] 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 [29] and FiniteFlow [30] to generate and solve the system of IBP identities based on Laporta's algorithm [31] 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 Refs.[32,33].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 blocktriangular 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 [34] based on power series expansion [35,36].The differential equations of master integrals with respect to x = 4m 2 /s are constructed using the aforementioned IBP reduction.To fix the boundary conditions, say at an arbitrary regular point x = 4/23, we utilize the auxiliary mass flow method [37][38][39][40] implemented in AMFlow [41].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 x → 0, 1 4 , 1 2 , 3 4 , 1 , where 0, 1/4 and 1 are singular points and 1/2 and 3/4 are regular points.The expansion at each point, say x 0 , offers a rapidly converging estimation, particularly within the range x ∈ (x 0 − 1/8, x 0 + 1/8).Consequently, the entire physical region x ∈ (0, 1) is comprehensively covered.
It is widely recognized that within dimensional regularization, the simultaneous satisfaction of the anticommutation relation {γ µ , γ 5 } = 0 and the cyclicity of the Dirac trace poses a challenge.In practical applications, the maintenance of the anticommutation relation is not only favored for computational simplification but is also valued for its potential to safeguard chiral symmetry and, consequently, gauge invariance [42].This approach necessitates a specific prescription for reading a fermion loop, which we choose the KKS scheme [43][44][45].In the KKS scheme, the ultimate result is defined as the average of all conceivable reading points, originating from both the head and tail of an axial-vector three-point fermion-gauge subgraph that encompasses the maximal one-particleirreducible non-singlet-type loop correction [45,46].In this context, the terms "singlet" and "non-singlet" are defined within the framework of Feynman diagrams, indicating whether the external current directly couples to a fermion loop within the relevant part of the Feynman diagram.This distinction is exemplified in Fig. 1 (b) and 1(a, c, d), respectively.For our calculation within the Standard Model, which is anomaly-free, no additional finite renormalization is required to restore symmetries when using the KKS scheme, at least to the order considered here [45,46].
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 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 MS αs and Z OS 2 .Heavy quark mass renormalization is more complicated because it is not multiplicatively renormalizable.We replace bare mass by renormalized mass via m b = m OS + δm, and keep on-shell condition of external heavy quark momentum p 2 = (m OS ) 2 untouched.In the spirit of perturbation theory, δm is a small quantity and can be expanded to any desired order.The expansion can be readily implemented either at the integrand level or, alternatively, at the integral level by formulating and solving differential equations with respect to δm.After the expansion of δm, we can substitute renormalization constants, which can be found in Refs.[47][48][49], to obtain the 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 α s .
Another key technique in our calculation is that the master integrals are calculated with numerical values of ϵ and the ϵ dependence is only reconstructed at the crosssection level, as proposed in Refs.[40,41].The advantage of this technique is that we do not need to manipulate Laurent expansions of ϵ during the intermediate stages of calculations, which significantly reduces the computational time.
Finally, we discuss how to compute differential crosssections based on the previously outlined methods.Using the invariant mass distribution of the heavy-quark pair as an example, the differential cross-section can be derived by inserting a delta function δ Q Q into the final state phase space integration.This integration can be approached similarly to the total cross-section calculation, noting that the delta function can be expressed as cut propagators using reverse unitarity [25][26][27].We highlight that the boundary conditions for differential equations of master integrals in differential cross-section computation can be determined by aligning them with integrals from the total cross-section.Specifically, we solve the differential equation and formulate master integrals as piecewise functions of M Q Q, represented by 13 deeply expanded power series with undetermined coefficients.By integrating over M Q Q, we collapse the cut propagator and arrive at integrals that belong to the integral families of the total cross-section, which are already known.Consequently, we can ascertain both the unknown coefficients and the piecewise functions.
Thanks to all these strategies mentioned above, the computational resources utilized in this work amount to less than 10 5 CPU core hours in total. Results.
-By combining everything together, we end up with final results at the NNNLO level, which are free of ultraviolet and infrared divergences as expected from the Kinoshita-Lee-Nauenberg theorem [50,51].Our result of the NNNLO total cross-section is expressed as a piecewise function of x = 4m 2 /s represented by 5 power series.Expanding these series to 40 orders enables us to achieve at least 10 correct digits in the physical region, with a relative error of approximately 2 −40 ∼ 10 −12 .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 = 172.69GeV [52] and set all other quarks as massless.Electromagnetic coupling is chosen as a fixed value α = 1/132.2,and strong coupling α s (µ) is running as a function of renormalization scale µ, which is computed using the RunDec package [53,54] with input value α Other electroweak parameters are chosen from Ref. [52].
Fig. 2 shows our result for the total cross-section of t t production, where LO, NLO, NNLO and NNNLO are shown in black, blue, green and red respectively, and the It can be found that the NNNLO correction significantly reduces the scale dependence.However, the NNNLO result near the production threshold, say for √ s < 370 GeV, still suffers from large uncertainty due to Coulomb interaction.Perturbative calculations become unreliable in this region, necessitating the application of resummation for further improvement [20,21,55].
The total cross-section can also be expressed in the form where the order α s , α 2 s and α 3 s corrections ∆ 1 , ∆ 2 and ∆ 3 are displayed in the lower panel of Fig. 2 as functions of the center-of-mass energy, respectively.The renormalization scale is set to √ s.The smallness of ∆ 3 confirms a good convergence of perturbative expansion with respect to α s .In Fig. 3, we further show the reduction of the scale dependence after including O(α 3 s ) correction by varying the renormalization scale in a larger range.It is found that, for a collision energy of 500 GeV, the scale dependence has been reduced from 0.60% at NNLO to 0.06% at NNNLO by varying the scale between µ = √ s/2 and µ = 2 √ s, which meets the precision requested by, e.g., International Linear Collider [10].The situation is analogous for b b or cc production cases, as illustrated in Fig. 4. In these instances, an onshell heavy-quark mass is chosen as m b = 4.8 GeV and m c = 1.5 GeV, respectively.Except the threshold region where perturbation is unreliable, theoretical uncertainties are significantly reduced, making them valuable for describing corresponding processes, like the so-called R ratio.It is important to note that, in our calculation, quarks with masses heavier than the considered heavy quark are excluded, while those with masses lighter than the considered quark are treated as massless.Additionally, our consideration is limited to diagrams mediated by an off-shell photon (γ * ) since the mass of the Z boson significantly exceeds the energy scale of interest here.As a first instance of NNNLO differential crosssections, Fig. 5 illustrates the invariant mass distribution of the top-quark pair for a collision energy of 500 GeV.The distribution increases as M t t becomes larger and becomes singular when M t t approaches its maximum value of √ s, where QCD radiations must be soft.In this region, the fixed-order prediction breaks down, prompting us to take an average over a 20 GeV bin.Although fixed-order perturbation theory also breaks down when M t t ∼ 2m t , this region is negligible due to the significant suppression by phase space considerations.On the whole, the pure O(α 3 s ) correction introduces a significant modification to the lower order result, falling within the range of approximately 9% to 13% for M t t ∈ [370, 480] GeV and about −1% for the bin at the end point.The scale uncertainties have been mitigated to approximately 5% across the majority of the distribution range.This represents a considerable impact for phenomenological study.
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 t total cross-section and invariant mass distribution at NNLO agrees with that obtained in Refs.[14,16].Specifically, we have replicated Fig. 3, Fig. 7, and the left plot of Fig. 10 in Ref. [14], and reproduced Fig. 1, Fig. 2, and the right plot of Fig. 5 in Ref. [16].Our three-loop contribution agrees with Refs.[22,23] to 7 digits provided there.Our non-singlet results at NNNLO agree with [18,19] within their uncertainty.1All master integrals involved in total cross-section, 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 with at least 10 digits precision.By integrating the invariant mass differential distribution at NNNLO, we found perfect agreement with total cross-section.
For the vector current contributions, we also calculate the forward scattering amplitude of a virtual photon at 4loop level without cut, and then relate its imaginary part to the desired cross-section by virtue of optical theorem.The results agree perfectly with that obtained using default methods.Note that using optical theorem could induce undesired final states, including states with no heavy quarks and states with four heavy quarks, which are subtracted by calculating virtual corrections and real radiations using the default methods.
To validate the appropriateness of the chosen γ 5 scheme, we performed a cross-check by computing the singlet part of the axial-axial contribution at NNNLO using the Larin scheme [56,57].Our focus was specifically on the axial-axial contribution, as the vector-vector contribution does not involve γ 5 , and the vector-axial contribution is zero.Moreover, for the non-singlet portion of the axial-axial contribution, the use of the anticommuting γ 5 scheme is unambiguous, because there will always be an even number of γ 5 in the Dirac trace.Consequently, our consideration is confined to the singlet part of the axial-axial contribution.Our findings indicate that the results obtained using the KKS scheme align with those obtained using the Larin scheme, following a finite renormalization in the latter scheme [56][57][58].
Summary.-In this Letter, we present the complete NNNLO QCD correction to the cross-section of heavy-quark pair production at lepton colliders.Additionally, we provide the first differential cross-section, specifically the invariant mass distribution, at the same order.Using top-quark pair production at a collision energy of 500 GeV as a benchmark, we observe a correction of approximately 0.1% for the total cross-section, while the correction for the invariant mass distribution ranges around 10% for most values away from the end point.The NNNLO QCD correction largely reduces the scale dependence.These results hold significant value for the precision testing of heavy-quark pair production.
The derivation of our results involves the integration of various powerful techniques, like reverse unitarity, finite field, block-triangular relations, and numerical sampling of ϵ.Our strategic approach is not only tailored to the specific process studied here but is also versatile enough to be applied to numerous other processes, including differential observables.This versatility holds promise for advancing perturbative calculations to unprecedented orders in various contexts.

FIG. 2 :
FIG.2: Total cross-section for t t production.Refer to the text for details.

FIG. 4 :
FIG. 4: Total cross-section for b b and cc production.