Reanalysis of the top-quark pair production via the $e^{+}e^{-}$ annihilation near the threshold region up to N³LO QCD corrections

The production of top-quark pair in electron-positron annihilation [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18] holds substantial relevance for both the validation of the Standard Model (SM) [19] and the pursuit of new physics beyond it. As the heaviest species within the SM, a comprehensive probe into the diverse properties of the top quark is anticipated via the exploration of its threshold production. For instance, pinpoint theoretical estimations of the threshold production cross section for $e^{+}e^{-}\to t\bar{t}$ facilitate the ascertainment of the top quark’s exact mass, which also serves as a vital observable for future high-energy electron-positron colliders such as the International Linear Collider (ILC) [20], the Circular Electron Positron Collider (CEPC) [21] and the Future Circular Collider [22], and etc.

Recently, the full next-to-next-to-next-to-leading order (N³LO) QCD corrections to the heavy quark pair production $e^{+}e^{-}\to\gamma^{*}\to Q\bar{Q}$ has been completed in Ref.[23], which still exhibits significant renormalization scale dependence. Especially, the Coulomb-like interaction between the quark and the anti-quark triggers a pronounced enhancement of the cross section close to the production threshold, which has been incorporated to all orders in perturbation theory. To refine the current theoretical calculation, it is essential to address the mentioned renormalization scale dependence and the Coulomb singularity near the threshold. This can be accomplished by investigating renormalization scale-setting problem, estimating unknown higher-order (UHO) contributions, and resumming all of the Coulomb corrections. Such an endeavor will not only lead to a more robust theoretical calculation but also provide a more precise perturbative QCD (pQCD) prediction on the top-quark pair production cross section near the threshold region.

The fixed-order pQCD predictions for observables are usually assumed to suffer from an uncertainty in fixing the renormalization scale. This uncertainty in making fixed-order predictions occurs because one usually assumes an arbitrary renormalization scale together with an arbitrary range to ascertain its uncertainty. This simple treatment of the pQCD series and the renormalization scale causes the coefficients of the QCD running coupling ($\alpha_{s}$) at each order to be strongly dependent on the choice of both the renormalization scale and the renormalization scheme. The $\alpha_{s}$-running behavior is governed by the renormalization group equation (RGE),

$\displaystyle\beta(\alpha_{s})=\frac{{\rm d}\alpha_{s}(\mu^{2})}{{\rm d}\ln\mu^{2}}=-\sum_{i=0}^{\infty}\beta_{i}\alpha_{s}^{i+2}(\mu^{2}),$

(1)

where the $\{\beta_{i}\}$-functions have been computed up to five-loop level in modified minimal-subtraction scheme ($\overline{\rm MS}$-scheme) [24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. Consequently, the $\{\beta_{i}\}$-terms emerged in pQCD series can be reabsorbed into $\alpha_{s}$ to accurately fix the correct magnitude of $\alpha_{s}$. In respect of this fact, the Principle of Maximum Conformality (PMC) [34, 36, 37, 35, 38] has been proposed to address the scale-setting problem in fixed-order pQCD calculation. It eliminates conventional renormalization scale ambiguity by using the above RGE recursively so as to set the correct $\alpha_{s}$ value and to achieve a well-matched (maximumly-approaching) conformal series between $\alpha_{s}$ and its expansion coefficients, yielding naturally scheme-independent theoretical predictions at any fixed-order, cf. the recent reviews [39, 40, 41, 42]. This treatment is a kind of resummation of all the known type of RG-involved $\{\beta_{i}\}$-terms. More specifically, the PMC single-scale-setting approach (PMCs) [43, 44] determines a global effective coupling $\alpha_{s}(Q_{*}^{2})$ ($Q_{*}$ is called as the PMC scale which can be viewed as the effective momentum flow of the process) for any fixed-order prediction by using the RGE, which also provides a solid way to extend the well-known Brodsky-Lepage-Mackenzie (BLM) approach [45] to all orders ¹¹1A short review of the development of PMC from BLM and replies to some wrong comments on the PMC can be found in Ref.[46].. Using the PMCs, it has been demonstrated that the PMC prediction is independent to any choice of renormalization scheme and scale [47], being consistent with the fundamental renormalization group approaches [48, 49, 50, 51] and also the self-consistency requirements of the renormalization group such as reflectivity, symmetry and transitivity [52].

It has been found that Coulomb-type corrections for the heavy-quark pair ($Q\bar{Q}$) production via the electron-positron annihilation have sizable contributions in the threshold region [53, 54], which similar to the hadronic production case [55] shall also be enhanced by factors of $\pi$ and the resultant PMC scale can be relatively small when the heavy quark velocity $v\to 0$ [56]. Thus the terms which are proportional to various powers of $(\pi/v)$ need to be treated carefully. In Ref.[56], the total cross section of $e^{+}e^{-}\to\gamma^{*}\to Q\bar{Q}$ near the threshold region has been investigated up to next-to-next-to leading order (N²LO) level. At the N³LO level and higher, the question will be much more involved due to the much more complex Coulomb structures [9, 10, 11]. In this paper, we will begin by partitioning the N³LO-level total cross section of $e^{+}e^{-}\to\gamma^{*}\to Q\bar{Q}$ into two components, Coulomb and non-Coulomb corrections. Subsequently, we will employ the PMCs approach to deal with each components individually. Following that, we will perform a resummation of the divergent Coulomb corrections using the well-established Sommerfeld-Gamow-Sakharov (SGS) factor [57, 58, 13, 59]. Then as an explicit example, we will give the numerical results for the case of top-quark pair production. Finally, we will utilize the Bayesian analysis (BA) [60, 61, 62, 63] to estimate the contributions from the UHO-terms.

The remaining parts of paper are organized as follows. In Sec.II, we present the total cross section of $e^{+}e^{-}\to\gamma^{*}\to Q\bar{Q}$ near the threshold region up to N³LO-level, and explain how to deal with its Coulomb and non-Coulomb components separately. And then we derive the analytical expression of the Coulomb part by using the integer relation finding algorithm, e.g. the so-called PSLQ algorithm [64, 65]. The results for the N³LO-level Coulomb coefficients will be given in the Appendix. In Sec.III, we give a concise introduction to the PMCs method. Numerical results and discussions for the top-quark pair production are given in Sec.IV. Sec.V is reserved for a summary.

Up to N³LO-level QCD corrections, the total cross section of $e^{+}e^{-}\to\gamma^{*}\to Q\bar{Q}$ can be written as

$\displaystyle\sigma(s)=\sigma_{0}\left(1+r_{1}\alpha_{s}+r_{2}\alpha_{s}^{2}+r_{3}\alpha_{s}^{3}+{\cal O}(\alpha_{s}^{4})\right),$

(2)

where $\alpha_{s}\equiv\alpha_{s}^{(n_{l})}$ is QCD running coupling in $\overline{\rm MS}$-scheme with $n_{l}$ represents the active number of light flavors, and $\sigma_{0}$ is the LO cross section which takes the form

$\displaystyle\sigma_{0}=N_{C}\frac{4\pi\alpha^{2}}{3s}\frac{v(3-v^{2})}{2}e_{Q}^{2}.$

(3)

Here $N_{C}=3$ for $SU(3)_{C}$ color group, $\alpha$ is electromagnetic coupling constant, $v=\sqrt{1-4m_{Q}^{2}/s}$ is heavy quark velocity, $m_{Q}$ is the on-shell mass of the heavy quark $Q$, $s$ is the squared $e^{+}e^{-}$ collision energy, and $e_{Q}$ is the charge of the heavy quark $Q$ in units of the electric charge, which equals to $+2/3$ for $Q=t$. It is convenient to express the perturbative coefficients $r_{1}$, $r_{2}$ and $r_{3}$ as

	$\displaystyle r_{1}$	$\displaystyle=\frac{1}{v}r_{1,v}+r_{1,+},$		(4)
	$\displaystyle r_{2}$	$\displaystyle=\frac{1}{v^{2}}r_{2,v^{2}}+\frac{1}{v}r_{2,v}+r_{2,+},$		(5)
	$\displaystyle r_{3}$	$\displaystyle=\frac{1}{v^{2}}r_{3,v^{2}}+\frac{1}{v}r_{3,v}+r_{3,+},$		(6)

where $r_{i,+}$ represents the terms of non-negative powers of $v$ in $r_{i}$ with $i=(1,2,3)$, respectively. Clearly, in the vicinity of the two-particle threshold, i.e., as the limit $v\to 0$, the total cross section exhibits the well-known Coulomb singularity, which is dominated by the negative power term of $v$ in Eqs.(4,5,6).

The analytical results of $r_{1}$, $r_{2}$ and all $n_{l}$-terms in $r_{3}$ can be found in Ref.[12, 14], and the semi-numerical expression of $r_{3}$ can be found in Ref.[15, 16]. In Ref.[23], numerical results up to N³LO-level have been expressed as asymptotic expansions of $x=4m_{Q}^{2}/s$, which have been expanded up to $40_{\rm th}$-orders. Using those numerical results, we are able to reconstruct the analytical results of $r_{3}$ by using the PSLQ algorithm [64, 65]. The PSLQ algorithm can be used to reconstruct analytical expression for given high-precision numerical values. As an explanation, the inputs to the PSLQ algorithm are high-precision numerical value denoted as $y$, and a vector $x=(x_{1},\cdots,x_{n})$ consisting of real numbers that are analytically well-defined, such as $1$, $\ln 2$, $\pi$, ${\rm Li}_{n}(1/2)$, $\zeta_{n}$ and so on. Subsequently, the PSLQ algorithm determines a set of integer coefficients $m_{0},m_{1},\cdots,m_{n}$, such that the linear combination with these integer coefficients is close to zero within the required numerical precision, i.e., $m_{0}y+\sum_{i=1}^{n}m_{i}x_{i}=0$. Therefore, the analytical expression of $y$ can be written as $y=-\sum_{i=1}^{n}\frac{m_{i}}{m_{0}}x_{i}$. The PSLQ algorithm can be implemented using the Mathematica package PolyLogTools [66]. Our results of those coefficients at the scale $\mu=\sqrt{s}$ are given in the Appendix.

Note that the coefficients in Eq.(2) are parameterized in terms of $\alpha_{s}^{(n_{l})}$, which can be transformed from the results parameterized in terms of $\alpha_{s}^{(n_{l}+1)}$ by relating the coupling for the theory with $n_{l}+1$ flavors to the theory with $n_{l}$ flavors via the following decoupling relationship,

$\displaystyle\frac{\alpha_{s}^{(n_{l}+1)}(\mu^{2})}{\alpha_{s}^{(n_{l})}(\mu^{2})}=1+\frac{1}{3}T_{F}\ln\frac{\mu^{2}}{m_{h}^{2}}\frac{\alpha_{s}^{(n_{l})}(\mu^{2})}{\pi}+T_{F}\left[-\frac{2}{9}C_{A}+\frac{15}{16}C_{F}+\left(\frac{5}{12}C_{A}+\frac{1}{4}C_{F}\right)\ln\frac{\mu^{2}}{m_{h}^{2}}+\frac{1}{9}T_{F}\ln^{2}\frac{\mu^{2}}{m_{h}^{2}}\right]\left(\frac{\alpha_{s}^{(n_{l})}(\mu^{2})}{\pi}\right)^{2},$

(7)

where $m_{h}$ is the on-shell mass of the decoupled heavy quark, $C_{A}=3$, $C_{F}=4/3$ and $T_{F}=1/2$.

The Coulomb singularity plays a dominant role for achieving precise total cross section near the threshold region ($v\to 0$). Following the suggestion of Refs.[12, 13, 59], we schematically factorize the total cross section (2) as the product of Coulomb and non-Coulomb parts, i.e.,

$\displaystyle\sigma=\sigma_{0}\times\mathcal{R}_{\rm C}\times\mathcal{R}_{\rm NC},$

(8)

where $\mathcal{R}_{\rm C}$ and $\mathcal{R}_{\rm NC}$ represent the Coulomb part and the non-Coulomb part of the total cross section, respectively, which can be expressed as

	$\displaystyle\mathcal{R}_{\rm C}$	$\displaystyle=1+r_{1}^{\rm C}\alpha_{s}^{\rm V}(sv^{2})+r_{2}^{\rm C}\alpha_{s}^{\rm V,2}(sv^{2})+r_{3}^{\rm C}\alpha_{s}^{\rm V,3}(sv^{2}),$		(9)
	$\displaystyle\mathcal{R}_{\rm NC}$	$\displaystyle=1+r_{1}^{\rm NC}\alpha_{s}(s)+r_{2}^{\rm NC}\alpha_{s}^{2}(s)+r_{3}^{\rm NC}\alpha_{s}^{3}(s).$		(10)

Here the QCD coupling $\alpha_{s}^{\rm V}(\vec{q}^{2})$ has been introduced for describing the interaction of the non-relativistic heavy quark-antiquark pair, which is defined as the effective charge in the following Coulomb-like potential [67, 68, 69, 70, 71, 72, 73, 74]:

$\displaystyle V(\vec{q}^{2})=-4\pi C_{F}\frac{\alpha_{s}^{\rm V}(\vec{q}^{2})}{\vec{q}^{2}},$

(11)

where $\alpha_{s}^{\rm V}(\vec{q}^{2})$ absorbs all the higher-order QCD corrections, which is related to the $\overline{\rm MS}$-scheme coupling via the following way

$\displaystyle\frac{\alpha_{s}^{\rm V}(\mu^{2})}{\alpha_{s}(\mu^{2})}=1+a_{1}\alpha_{s}(\mu^{2})+a_{2}\alpha_{s}^{2}(\mu^{2})+\mathcal{O}(\alpha_{s}^{3}),$

(12)

where the expansion coefficients $a_{1}$ and $a_{2}$ have been calculated in Ref.[74], and their specific forms are as follows:

	$\displaystyle a_{1}=$	$\displaystyle\,\frac{1}{4\pi}\left(\frac{31}{9}C_{A}-\frac{20}{9}T_{F}n_{l}\right),$		(13)
	$\displaystyle a_{2}=$	$\displaystyle\,\frac{1}{(4\pi)^{2}}\bigg{[}\left(\frac{4343}{162}+4\pi^{2}-\frac{\pi^{2}}{4}+\frac{22}{3}\zeta_{3}\right)C_{A}^{2}$
		$\displaystyle\,-\left(\frac{1798}{81}+\frac{56}{3}\zeta_{3}\right)C_{A}T_{F}n_{l}-\left(\frac{55}{3}-16\zeta_{3}\right)C_{F}T_{F}n_{l}$
		$\displaystyle\,+\left(\frac{20}{9}T_{F}n_{l}\right)^{2}\bigg{]},$		(14)

where $\zeta_{n}=\zeta(n)$ is Riemann’s Zeta function. Then we can obtain

	$\displaystyle r_{1}^{\rm C}$	$\displaystyle=C_{F}\frac{\pi}{2v},$		(15)
	$\displaystyle r_{2}^{\rm C}$	$\displaystyle=C_{F}^{2}\frac{\pi^{2}}{12v^{2}},$		(16)
	$\displaystyle r_{3}^{\rm C}$	$\displaystyle=C_{F}\left(\frac{\pi^{3}}{3v}\beta_{0}^{2}(n_{l})-C_{F}\frac{2\zeta_{3}}{v^{2}}\beta_{0}(n_{l})\right),$		(17)

and the coefficients $r_{i}^{\rm NC}$ can be obtained by expanding $\sigma(s)/(\sigma_{0}\mathcal{R}_{\rm C})$ over $\alpha_{s}$ and $v$.

In the previous literature [43, 44], it has been found that the PMCs approach provides an all-orders single-scale-setting method to fix the conventional renormalization scale ambiguity, which adopts the standard RGE and the general QCD degeneracy relations among different orders [75] to transform the RG-involved $n_{f}$-series into the $\{\beta_{i}\}$-series. In this subsection, we give a simple introduction of PMCs for self-consistency.

If the pQCD approximant $\rho(Q)$ starts and ends at the $n_{\rm th}$ and $(n+p-1)_{\rm th}$ order of $\alpha_{s}$, one has

$\displaystyle\rho(Q)=\sum_{k=1}^{p}c_{k}\alpha_{s}^{n+k-1}(\mu_{R}^{2}),$

(18)

where $\mu_{R}$ is an arbitrary initial choice of renormalization scale, $Q$ represents the kinetic scale at which the observable is measured or the typical momentum flow of the reaction. The perturbative coefficients $c_{k}$ can be rewritten as the following forms by using the general QCD degeneracy relations:

	$\displaystyle c_{1}$	$\displaystyle=c_{1,0},$		(19)
	$\displaystyle c_{2}$	$\displaystyle=c_{2,0}+n\beta_{0}c_{2,1},$		(20)
	$\displaystyle c_{3}$	$\displaystyle=c_{3,0}+(n+1)\beta_{0}c_{3,1}+\frac{n(n+1)}{2}\beta_{0}^{2}c_{3,2}+n\beta_{1}c_{2,1},$		(21)
		$\displaystyle\cdots.$

Specifically, the fixed-order pQCD series $\rho$ can be further rewritten as the following compact form:

$\displaystyle\rho(Q)=\sum_{k=1}^{p}\hat{c}_{k,0}\alpha_{s}^{n+k-1}(\mu_{R}^{2})+\sum_{k=1}^{p-1}\left[(n+k-1)\beta(\alpha_{s})\alpha_{s}^{n+k-2}(\mu_{R}^{2})\sum_{i=1}^{p-k}(-1)^{i}\Delta_{n,k}^{(i-1)}\sum_{j=0}^{i}C_{i}^{j}\hat{c}_{k+i-j,i-j}L_{\mu_{R}}^{j}\right],$

(22)

where $L_{\mu_{R}}\equiv\ln(\mu_{R}^{2}/Q^{2})$, and the combinatorial coefficients are $C_{i}^{j}=j!/i!(j-i)!$. $\hat{c}_{i,0}=c_{i,0}$ is the conformal coefficients, $\hat{c}_{i,j\geq 1}=c_{i,j}|_{\mu_{R}=Q}$ is non-conformal coefficients, $c_{i,j}$ and $\hat{c}_{i,j}$ can be connected as the relation $c_{i,j}=\sum_{k=0}^{j}C_{j}^{k}\hat{c}_{i-k,j-k}L_{\mu_{R}}^{k}$. $\Delta_{n,k}^{(i-1)}$ is related to the running of the RGE, which is defined as

$\displaystyle\Delta_{n,k}^{(i-1)}(\mu_{R})=\frac{1}{(n+k-1)\beta\alpha_{s}^{n+k-2}}\frac{1}{i!}\frac{{\rm d}^{i}\alpha_{s}^{n+k-1}}{({\rm d}\ln\mu^{2})^{i}}\bigg{|}_{\mu=\mu_{R}},$

(23)

we thus have

	$\displaystyle\Delta_{n,k}^{(0)}$	$\displaystyle=1,$		(24)
	$\displaystyle\Delta_{n,k}^{(1)}$	$\displaystyle=-\frac{1}{2!}\sum_{i=0}(n+k+i)\beta_{i}\alpha_{s}^{i+1},$		(25)
	$\displaystyle\Delta_{n,k}^{(2)}$	$\displaystyle=\frac{1}{3!}\sum_{j=0}\sum_{i=0}^{j}(n+k+i)(n+k+j+1)\beta_{i}\beta_{j-i}\alpha_{s}^{j+2},$		(26)
		$\displaystyle\cdots.$

According to the PMCs procedure, all non-conformal terms in Eq.(22) should be absorbed into $\alpha_{s}(Q_{*}^{2})$, where $Q_{*}$ is a global effective scale and its value can be determined by requiring all non-conformal terms vanish in Eq.(22), i.e., one can obtain PMC scale $Q_{*}$ by solving the following equation:

	$\displaystyle\sum_{k=1}^{\ell+1}$	$\displaystyle\,\bigg{[}(n+k-1)\alpha_{s}^{k-1}(Q_{*}^{2})\sum_{i=1}^{\ell-k+2}(-1)^{i}\Delta_{n,k}^{(i-1)}(Q_{*})$
		$\displaystyle\,\times\sum_{j=0}^{i}C_{i}^{j}\hat{c}_{k+i-j,i-j}L_{Q_{*}}^{j}\bigg{]}=0,$		(27)

where $\ell$ indicates that the PMC scale $Q_{*}$ can be fixed up to N^ℓLL accuracy. Since the $\{\beta_{i}\}$-terms in the perturbative series (22) have been resummed to determine the accurate value of $\alpha_{s}$, and the resultant PMC series is free of divergent renormalon terms and becomes scheme- and scale-independent pQCD conformal series [47], e.g.

$\displaystyle\rho|_{\rm PMC}=\sum_{k=1}^{p}\hat{c}_{k,0}\alpha_{s}^{n+k-1}(Q_{*}^{2}).$

(28)

Practically, to achieve the analytic expression, the PMC scale $Q_{*}$ is expressed as the form of $\ln\left(Q_{*}^{2}/Q^{2}\right)$ by expanding Eq.(III) in terms of $\alpha_{s}(Q_{*}^{2})$ or $\alpha_{s}(Q^{2})$ up to the required N^ℓLL accuracy [43, 44]. In this paper, we will directly solve Eq.(III) numerically so as to derive more accurate PMC scales $Q_{*,\rm C}$ and $Q_{*,\rm NC}$ ²²2Due to double suppression from both the $\alpha_{s}$-suppression and the exponential suppression, the usual analytic treatment is accurate enough. For the present case, the N^ℓLL-level expansion series of $\ln\left(Q_{*}^{2}/Q^{2}\right)$ does not convergent enough in the threshold region; It is thus helpful to keep all known divergent $1/v$-terms so as to achieve a more accurate prediction..

We are now ready to deal with the Coulomb and non-Coulomb corrections within the PMC framework.

For the Coulomb part $\mathcal{R}_{\rm C}$ (9), after applying the PMC, we obtain the resultant PMC series, i.e.

	$\displaystyle\mathcal{R}_{\rm C}|_{\rm PMC}$	$\displaystyle=$	$\displaystyle 1+\sum_{i}\hat{r}_{i,0}^{\rm C}\alpha_{s}^{{\rm V},i}(Q_{*,\rm C}^{2}),$		(29)
		$\displaystyle=$	$\displaystyle 1+\frac{\pi}{2v}C_{F}\alpha_{s}^{\rm V}(Q_{*,\rm C}^{2})+\frac{\pi^{2}}{12v^{2}}C_{F}^{2}\alpha_{s}^{\rm V,2}(Q_{*,\rm C}^{2})$
			$\displaystyle+{\cal O}(\alpha_{s}^{\rm V,4}),$

whose first three conformal coefficients $\hat{r}_{i,0}^{\rm C}$ can be read from the coefficients $r_{i}^{\rm C}$ given in Eqs.(15,16,17) by applying the degeneracy relations (19-21) with the help of the replacements $c_{i}\to r_{i}^{\rm C}$. It is interesting to find that since $\hat{r}_{3,0}^{\rm C}=0$, the N³LO Coulomb coefficient $r_{3}^{\rm C}$ (17) is exactly non-conformal, which can be resummed through the RGE and determines an effective $\alpha_{s}$, giving the required effective momentum flow near the threshold region (corresponding to the PMC scale $Q_{*,C}$). That is, the resultant PMC Coulomb series can be resummed as a Sommerfeld-Gamow-Sakharov (SGS) factor [12, 13, 59] ³³3The SGS factor (30) can be expanded over $X$ as $1+{X}/{2}+{X^{2}}/{12}-{X^{4}}/{720}+{\cal O}(X^{5})$, in which the $X^{3}$-coefficient is exactly zero. As shown by Eq.(17), the coefficient $r_{3}^{C}$ is non-zero, thus the N³LO-level Coulomb part of the series (9) cannot be resummed into a simple SGS factor. Thus we will not provide the resummed Coulomb results for the predictions under conventional scale-setting approach.:

$$\mathcal{R}_{\rm C}|_{\rm PMC}^{\rm resum}=\frac{X}{1-e^{-X}},\;\;X=\frac{\pi}{v}C_{F}\alpha_{s}^{\rm V}(Q_{*,\rm C}^{2}).$$

(30)

The $1/v$ in the numerator can be canceled by the $v$ factor in the LO cross section $\sigma_{0}$, thus the Coulomb singularity near the threshold can be eliminated.

Similarly, for the non-Coulomb part (9), its PMC conformal series up to N³LO can be achieved and be written in the following form:

$\displaystyle\mathcal{R}_{\rm NC}|_{\rm PMC}=1+\sum_{i=1}^{3}\hat{r}_{i,0}^{\rm NC}\alpha_{s}^{i}(Q_{*,\rm NC}),$

(31)

where conformal coefficients $\hat{r}_{i,0}^{\rm NC}$ can be obtained in the same way as that of $\hat{r}_{i,0}^{\rm C}$.

It is worth mentioning that, even though the Coulomb divergence has been eliminated, $\ln v$-terms still remain in the $v^{0}$ terms of $\mathcal{R}_{\rm NC}$. However, when compared to the terms with negative powers of $v$, these terms are relatively more benign, and therefore, further treatment of those terms is not addressed in this article.

In the paper, we have presented an improved analysis for the total cross-section of $e^{+}e^{-}\to\gamma^{*}\to t\bar{t}$ near the threshold region up to N³LO QCD corrections. Analytical expression for the N³LO-level Coulomb-terms has been achieved by employing the PSLQ algorithm, which agrees well with the recently given numerical values in Ref.[23]. The PMCs approach has been applied to eliminate the large renormalization scale uncertainties for the Coulomb and non-Coulomb terms separately. After that, the Coulomb-terms have been resummed as a Sommerfeld-Gamow-Sakharov factor, which greatly improves the reliability and precision of total cross-section near the threshold region.

Our results confirm the correctness and the importance of the PMC procedures for achieving scheme-and-scale invariant pQCD predictions, which can also be ensured by the commensurate scale relations among different schemes [79, 80] and agrees with the renormalization group invariance [81, 82, 83, 84]. The elimination of the uncertainty in setting the renormalization scale for pQCD predictions using the PMCs, together with the reliable estimate for the UHO contributions obtained using the Bayesian analysis, greatly increases the precision of collider tests of the Standard Model and thus the sensitivity to new phenomena.

For the sake of completeness and convenience, we give the NLO and N²LO coefficients $r_{1}$ and $r_{2}$ at the scale $\mu=\sqrt{s}$ in the following,

	$\displaystyle r_{1}=$	$\displaystyle\,\frac{C_{F}}{\pi}\bigg{\{}\frac{1+v^{2}}{v^{2}}\left[4{\rm Li}_{2}\left(\frac{1-v}{1+v}\right)+2{\rm Li}_{2}\left(-\frac{1-v}{1+v}\right)+\ln\left(\frac{1-v}{1+v}\right)\left(\ln\left(\frac{2}{1+v}\right)+2\ln\left(\frac{2v}{1+v}\right)\right)\right]$
		$\displaystyle\qquad-2\left(\ln\left(\frac{2}{1+v}\right)+2\ln\left(\frac{2v}{1+v}\right)\right)-\frac{(1-v)(33-39v-17v^{2}+7v^{3})}{8v(3-v^{2})}\ln\left(\frac{1-v}{1+v}\right)+\frac{3(5-3v^{2})}{4(3-v^{2})}\bigg{\}}$		(38)
	$\displaystyle=$	$\displaystyle\,\frac{C_{F}}{\pi}\left[\frac{\pi^{2}}{2v}-4+\frac{\pi^{2}}{2}v+\left(-\frac{86}{9}+8\ln 2+\frac{16}{3}\ln v\right)v^{2}+\mathcal{O}(v^{4})\right],$		(39)
	$\displaystyle r_{2}=$	$\displaystyle\,C_{F}\bigg{\{}\frac{\pi^{2}C_{F}}{12v^{2}}+\frac{1}{v}\left[-2C_{F}+C_{A}\left(\frac{31}{72}-\frac{11}{12}\ln v\right)+T_{F}n_{l}\left(-\frac{5}{18}+\frac{1}{3}\ln v\right)\right]-\left(C_{A}+\frac{2}{3}C_{F}\right)\ln v$
		$\displaystyle\qquad+C_{F}\left(-\frac{35}{18}+\frac{39}{4\pi^{2}}+\frac{\pi^{2}}{6}+\frac{4}{3}\ln 2-\frac{1}{\pi^{2}}\zeta_{3}\right)+C_{A}\left(\frac{179}{72}-\frac{151}{36\pi^{2}}-\left(\frac{8}{3}+\frac{22}{3\pi^{2}}\right)\ln 2-\frac{13}{2\pi^{2}}\zeta_{3}\right)$
		$\displaystyle\qquad+T_{F}n_{l}\left(\frac{11}{9\pi^{2}}+\frac{8}{3\pi^{2}}\ln 2\right)+T_{F}n_{h}\left(-\frac{4}{9}+\frac{44}{9\pi^{2}}\right)$
		$\displaystyle\qquad+\left[\frac{839}{216}+\frac{173}{18}\ln 2+\frac{101}{36}\ln v+T_{F}n_{l}\left(-\frac{13}{18}+\frac{1}{3}\ln v\right)+T_{F}n_{h}\left(\frac{1}{6}+\frac{1}{3}\ln 2\right)\right]v+\mathcal{O}(v^{2})\bigg{\}},$		(40)

where ${\rm Li}_{n}(z)=\sum_{k=1}^{\infty}z^{k}/k^{n}$ is polylogarithm function. As for the N³LO coefficient $r_{3}$, the analytic expression of its coefficient $r_{3,+}$ with non-negative power terms of $v$ is very tedious and hard to give its analytic form, so we directly adopt its numerical form given by Ref.[23].

As for the N³LO-level coefficients $r_{3,v}$ and $r_{3,v^{2}}$ with negative power terms $1/v$ and $1/v^{2}$ that are helpful for fixing the PMC scale of the Coulomb part and for getting the SGS factor, we adopt the PSLQ algorithm with the help of Mathematica package PolyLogTools [66] so as to derive their analytic form, e.g.

	$\displaystyle r_{3,v^{2}}=$	$\displaystyle\,C_{F}\left[-\frac{\pi}{3}C_{F}^{2}+C_{F}C_{A}\left(\frac{31\pi}{216}-\frac{11}{6\pi}\zeta_{3}-\frac{11\pi}{36}\ln v\right)+C_{F}T_{F}n_{l}\left(-\frac{5\pi}{54}+\frac{2}{3\pi}\zeta_{3}+\frac{\pi}{9}\ln v\right)\right],$		(41)
	$\displaystyle r_{3,v}=$	$\displaystyle\,C_{F}\bigg{\{}C_{F}^{2}\left(\frac{39}{8\pi}-\frac{35\pi}{36}+\pi\ln 2-\frac{\zeta_{3}}{2\pi}-\frac{\pi}{3}\ln(2v)\right)+C_{A}^{2}\left(\frac{4343}{5184\pi}+\frac{175\pi}{432}-\frac{\pi^{3}}{128}+\frac{11}{48\pi}\zeta_{3}-\frac{247}{108\pi}\ln v+\frac{121}{72\pi}\ln^{2}v\right)$
		$\displaystyle\qquad+C_{A}C_{F}\left[-\frac{275}{72\pi}+\frac{179\pi}{144}-\left(\frac{4\pi}{3}+\frac{11}{3\pi}\right)\ln 2-\frac{13}{4\pi}\zeta_{3}+\left(\frac{11}{3\pi}-\frac{\pi}{2}\right)\ln v\right]+C_{F}T_{F}n_{h}\left(\frac{22}{9\pi}-\frac{2\pi}{9}\right)$
		$\displaystyle\qquad+T_{F}n_{l}\left[C_{F}\left(\frac{331}{288\pi}+\frac{\zeta_{3}}{2\pi}+\frac{4}{3\pi}\ln 2-\frac{13}{12\pi}\ln v\right)+C_{A}\left(-\frac{899}{1296\pi}-\frac{11\pi}{54}-\frac{7\zeta_{3}}{12\pi}+\frac{217}{108\pi}\ln v-\frac{11}{9\pi}\ln^{2}v\right)\right]$
		$\displaystyle\qquad+T_{F}^{2}n_{l}^{2}\left[\frac{25}{162\pi}+\frac{\pi}{27}-\frac{10}{27\pi}\ln v+\frac{2}{9\pi}\ln^{2}v\right]\bigg{\}}.$		(42)

It should be pointed out that the correctness of the PSLQ algorithm can be further tested by comparing the analytical results of $r_{1}$ and $r_{2}$ given by Ref.[12, 14] with the numerical results given in Ref.[23]; and by comparing with the above $r_{3}$ with the semi-numerical expression of $r_{3}$ given in Refs.[15, 16].