Semi-analytical calculation of gluon fragmentation into ${^{1}\hspace{-0.6mm}S_{0}^{[1,8]}}$ quarkonia at next-to-leading order

We calculate the NLO corrections for the gluon fragmentation functions to a heavy quark-antiquark pair in ${^{1}\hspace{-0.6mm}S_{0}^{[1]}}$ or ${^{1}\hspace{-0.6mm}S_{0}^{[8]}}$ state within NRQCD factorization. We use integration-by-parts reduction to reduce the original expression to simpler master integrals (MIs), and then set up differential equations for these MIs. After calculating the boundary conditions, MIs can be obtained by solving the differential equations numerically. Our results are expressed in terms of asymptotic expansions at singular points of $z$ (light-cone momentum fraction carried by the quark-antiquark pair), which can not only give FFs results with very high precision at any value of $z$, but also provide fully analytical structure at these singularities. We find that the NLO corrections are significant, with K-factors larger than 2 in most regions. The NLO corrections may have important impact on heavy quarkonia (e.g. $\eta_c$ and $J/\psi$) production at the LHC.


Introduction
Study of heavy quarkonium production is important to understand both perturbative and non-perturbative physics in QCD. Currently, the most widely used theory for quarkonium production is the non-relativistic QCD (NRQCD) factorization [1]. Although many important processes have been calculated to next-to-leading order in α s expansion , there are still some notable difficulties in quarkonium production within the NRQCD framework (see, e.g. [26]). To further explore the quarkonium production mechanism, it may be better to study quarkonium production at high transverse momentum p T region, where long-distance interactions between quarkonium and initial-state particles are suppressed and thus factorization is easier to hold.

JHEP04(2019)116
The inclusive production differential cross section of a specific hadron H at high p T can be calculated in collinear factorization [27], dσ A+B→H(p T )+X = i dσ A+B→i(p T /z)+X ⊗ D i→H (z, µ) + O(1/p 2 T ) , (1.1) where i sums over all quarks and gluons, z is the light-cone momentum fraction carried by H with respect to the parent parton i, and A and B are colliding particles whose effect should be further factorized into partons if they are hadrons. dσ A+B→i(p T /z)+X are perturbatively calculable hard parts, while D i→H (z, µ) are non-perturbative but universal fragmentation functions (FFs) describing the probability of parton hadronizing to H with momentum fraction z. For quarkonium production, O(1/p 2 T ) contributions can be further factorized to double parton FFs [28][29][30][31][32]. In both single parton FFs and double parton FFs, there is a collinear factorization scale µ dependence, and this dependence can be canceled between hard parts and FFs perturbatively order by order, leaving physical differential cross section to be independent of the scale. The evolution of single parton FFs with respect to µ are controlled by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation [33][34][35], and similar evolution equations for double parton FFs are calculated in [29]. With these evolution equations, the only unknown information for FFs are their values at a chosen factorization scale µ = µ f .
When µ f is close to the quarkonium mass m H , it is natural to calculate FFs via NRQCD factorization. For single parton FFs that will be considered in this paper, we have where d i→QQ(n) represent the perturbative calculable short-distance coefficients (SDCs) to produce a heavy quark-antiquark pair QQ with quantum number n, and Ō H n are normalized long-distance matrix elements (LDMEs). 1 The quantum number is usually expressed in spectroscopic notation n = 2S+1 L [c] J , with c = 1, 8 respectively for color-singlet state or color-octet state. According to velocity scaling rule [1], Ō H n is usually suppressed if L is too large. Therefore, the most important states for phenomenological purpose are S-wave and P -wave states. Because LDMEs are supposed to be process independent, they can be determined by fitting experimental data, while SDCs need to be calculated perturbatively.

JHEP04(2019)116
SDCs of g → QQ( 1 S [1] 0 ) + X, which involves not only tree-level diagrams but also one-loop diagrams. Numerical results for this process have been calculated in ref. [54]. Considering the complicity of the calculation, an independent check by another group is badly needed.
As 1 S [1] 0 is the dominant Fock state for η c,b , the FF g → QQ( 1 S [1] 0 ) + X is important to study η c,b production at high transverse momentum at LHC [55]. At the LHC, we have even much more data of J/ψ production at high transverse momentum. Theoretical studies [13,17,56,57] show that 1 S [8] 0 channel may be crucial to explain the J/ψ data. To calculate 1 S [8] 0 contribution precisely, we need to calculate the FF g → QQ( 1 S [8] 0 ) + X to at least NLO.
In this paper, we aim to calculate NLO SDCs of FFs of g → QQ( 1 S [1] 0 ) + X and g → QQ( 1 S [8] 0 ) + X to high precision using similar methods in our previous paper [52]. With sufficient numerical precision, analytical results can in principle be extracted by using PSLQ algorithm. The rest of the paper is organized as following. In section 2, we first introduce the definition of SDCs of quarkonium FFs, including projection operators and Feynman rules related to gauge link, and then give the LO results. NLO corrections include real emission Feynman diagrams and one-loop Feynman diagrams, the calculation of them will be presented in section 3 and section 4, respectively. In the calculation, we use integration-by-part (IBP) reduction [58][59][60][61][62] to express both real contributions and virtual contributions in terms of linear combination of a small set of simpler integrals, which are usually called master integrals (MIs). High precision MIs can be obtained by solving differential equations of MIs numerically. Renormalization will be presented in section 5. After renormalization, the obtained SDCs are free of ultraviolet (UV) and infrared (IR) divergences. Final results and discussions will be given in section 6. We find that our results for g → QQ( 1 S [1] 0 )+X seem to be different from that calculated in ref. [54], while our results for g → QQ( 1 S [8] 0 ) + X are new. Finally, high precision results and some technical details will be given in appendices.

Definitions
The definition of FF from a gluon to a hadron (quarkonium) is given by Collins and Soper [63], where G µν is the gluon field-strength operator, P and P c are respectively the momenta of the produced hadron H and the initial-state fragmenting gluon g, and z = P + /P + c is the ratio of momenta along the "+" direction. It is convenient to choose the frame in which the hadron has zero transverse momentum, P = zP + c , m 2 H /(2zP + c ), 0 ⊥ , with P 2 = 2P + P − = m 2 H . The projection operator P H(P ) is defined by where X sums over all unobserved particles. The gauge link E(x − ) is an eikonal operator that involves a path-ordered exponential of gluon field operators along a light-like path, where g s = √ 4πα s is the QCD coupling constant and A µ (x) is the matrix-valued gluon field in the adjoint representation: . From this definition, we can derive Feynman rules related to gauge link, which are shown in figure 1, where n = (0, 1 − , 0 ⊥ ), K and P denote momenta, µ and ν denote Lorentz indexes, and a, b and c denote color indices.
With these Feynman rules, we can obtain the amplitude of all Feynman diagrams denoted as M λ Q λQλ 0 λ i (P, k i , m Q ), where λ Q and λQ are respectively spins of produced onshell heavy quark and heavy antiquark, λ 0 and λ i (i = 1, 2, . . . ) are spins of the initial-state virtual gluon and final-state unobserved light particles, respectively, k i are the momenta of final-state light particles, and m Q is the heavy quark mass. For the processes of gluon fragmenting to S-wave quarkonium, the relative momentum between the QQ pair can be chosen as 0 directly at the lowest order in velocity expansion, and thus it does not appear in the amplitude. If we project the free QQ pair to specific states 1 S where Γ c , Γ 5 are the projection operators defined as

JHEP04(2019)116
where P 2 = M 2 = 4m 2 Q . By summing over spin and color of initial-state and final-state particles, we get the squared amplitude Then the SDCs for gluon fragmenting to spin-singlet S-wave quarkonium can be written as where S is the symmetry factor for final-state particles.
To be convenient, we extract the dependence on m Q explicitly by rescaling momenta in the delta function in eq. (2.8) by M , Thus the phase space in eq. (2.8) changes to where n is the number of final-state light particles, and dΦ is similar to dΦ by changing all momenta to the dimensionless ones. If we further denotê we get a similar relation as that in eq. (2.7), which means that the same SDCs can be obtained by replacing all momenta by their corresponding rescaled ones. In the rest of the paper, we will only use the rescaled momenta, but omitting the "ˆ" for simplicity.

LO SDCs
The Feynman diagrams of gluon fragmenting into 1 S where a equals 0 or 1/2, k is the momentum of the emitted gluon with k + = (1 − z)P + /z and k − = k 2 ⊥ /(2k + ), and (2.14) These integrals can be performed easily.
Then we get LO SDCs: where µ r is the renormalization scale, d [1] LO and d [8] LO respectively denote SDCs of gluon fragmenting into 1 S [1] 0 and 1 S [8] 0 states, and The color-singlet result and color-octet result are consistent with refs. [54] and [45], respectively.

Reduction to MIs
Real NLO corrections to FFs of g → QQ( 1 S SDCs can be expressed as linear combinations of integrals of the form ) + gg. The other diagrams can be obtained by permuting the heavy quark and anti-quark or the two emitted gluons. where a i are integers, k 1 and k 2 are momenta of the final-state light particles, the phase space dΦ real is defined in eq. (2.8) with S = 2, and 2) In eq. (3.1), we safely ignore infinitesimal imaginary parts in denominators because E i (i = 1, · · · 11) are positive-definite and that SDCs are well regularized by dimensional regularization. The later condition implies that only the region where all E i (i = 1, · · · 11) are not too small can contribute to the phase space integration. Note that, for qq pair emission, although the symmetry factor should be 1, we can also express the SDCs as linear combinations of integrals in eq. (3.1).
To take advantage of multi-loop techniques, we express delta functions by propagator denominators,

JHEP04(2019)116
We replace the three delta functions in eq. (3.1) following the above rule, and denote Then each phase space integral in eq. (3.1) is translated to 8 loop integrals, with either positive or negative infinitesimal imaginary parts in new denominators.
If we forget about infinitesimal imaginary parts in denominators for the moment, we need to deal with loop integrals with integers a i , which can be expressed in terms of corresponding simpler MIs by using IBP reduction [58][59][60][61][62]. MIs are also the same kind of integrals, but usually with smaller a i . Note that, we can always choose MIs with powers of E 12 , E 13 and E 14 being no larger than 1. For MIs with integrand involving 1 E 12 , we can replace the denominator by δ + (k 2 1 ) considering the relation eq. (3.3), while for MIs with integrand E −a 12

12
(a 12 ≤ 0) we can set it to zero. Similar replacement can be done for E 13 and E 14 . Therefore, all MIs for loop integration are changed back to corresponding MIs for phase space integration defined in eq. (3.1). Once these MIs are also calculated, we obtain final results of real corrections.
In the above procedure, we actually assume that IBP reduction relations are independent of infinitesimal imaginary parts in denominators. This assumption, unfortunately, does not always hold. If one or more integrals cannot be fully regularized by dimensional regularization, one may get wrong final results. In the appendix A, we will discuss this problem in more details, and then propose a solution. Eventually, the above procedure is justified with a small modification.

Calculation of MIs
To calculate these MIs, we use differential equations (DEs) method [64][65][66][67][68][69][70][71][72][73][74][75][76], which has also been used in our previous paper [52] to calculate SDCs of g → QQ( 3 S 1 ) + X. We get 95 MIs using the IBP reduction program FIRE5 [62], without using the symmetry rules. We set up DEs by first differentiating these MIs I k (k = 1, . . . , 95) with respect to z, and then reducing the resulted integrals to MIs again by using IBP reduction, which results in where I represents the vector of MIs I k , and A is a 95 × 95 matrix whose elements are rational functions of z and . Having the DEs, we also need boundary conditions of I k to fully determine these MIs. We choose the boundary at z → 1, and calculate the boundary conditions in appendix C.
With boundary conditions, we can solve the DEs to obtain MIs at any value of z. One possible choice is to solve the DEs analytically, which can be done by transforming DEs into canonical form (or -form) [68,69]. In this way, we successfully express MIs in terms of Goncharov polylogarithms (GPLs) [77]. All obtained GPLs have weights at most three, and they can be expressed in terms of logarithms and classical polylogarithms Li n (z), (n ≤ 3) [78]. Even though, the obtained analytical expression is too long to present in this paper. Furthermore, for virtual correction, boundary conditions are hard to calculate analytically.
Another choice is to solve DEs numerically, which is a well-studied mathematical problem. DEs can help to do asymptotic expansions of MIs around any point z = z 0 . Because Feynman integrals have Feynman parametric representation, their asymptotic expansions have the form (see e.g. ref. [79]) where s is a linear function of , n s is an integer determined by s, I s i j k ( ) are functions of , and the radius of convergence of the summation over j is usually determined by the nearest singular point. For the special case when z 0 is an analytical point, we have s = n s = 0. When z 0 is a singular point, different values of s and i correspond to different regions of MIs, which are independent of each other. Therefore, each region satisfies the same DEs as the original MIs, and the DEs can generate recurrence relations to express I s i j k ( ) in terms of I s 0 0 k ( ) for each fixed s, i and . It implies that, when calculating boundary conditions in appendix C, we only need to calculate I s 0 0 k ( ) for each region. In practice, as we are only interested in MIs up to a fixed order in expansion, we will do a Laurent expansion of in both (z − z 0 ) s and I s i j k ( ). As it is clear, singular points play important role in the procedure of solving DEs numerically. There are 12 singular points in the DEs (3.6) for real corrections, which are located at z = 0, 1/2, ±1, ±2, ±4, ±2i, 1 ± i, as shown in figure 5. singularities are far enough from the physical region. Among these three singularities, the point z = 1/2 is in fact a removable singularity. However, as we will discuss in appendix B, this singularity determines the radius of convergence of the asymptotic expansion at z = 0 and 1. We thus estimate values of MIs in regions 0 ∼ 1/4, 1/4 ∼ 3/4 and 3/4 ∼ 1 respectively by the asymptotic expansions of MIs at z = 0, 1/2 and 1. For example, if we want to obtain values in the physical region with precision about 15 digits, we should calculate the expansion in eq. (3.7) with j to as large as 50.

Virtual NLO corrections
Some diagrams that contributed to virtual NLO corrections to FFs of g → QQ( 1 S [1,8] 0 ) + X are shown in figure 6. The other diagrams are either self-energy diagrams for external legs (including initial virtual gluon), or they can be obtained by permuting the heavy quark and anti-quark.
SDCs of the virtual corrections can be expressed as linear combination of integrals of the form 1) where z = P + /(k + + P + ), a i are integers, k is the momentum of the final-state gluon, l is JHEP04(2019)116 the loop momenta, and (4.2) Similar to real corrections, by replacing δ functions using eq. (3.3), integrals in eq. (4.1) can be reduced to corresponding simpler MIs, the number of which is 66. DEs for these MIs can also be set up.
As for real corrections, asymptotic expansion of virtual-correction MIs at any point z = z 0 can be obtained in the form of eq. (3.7) with the help of DEs, once we have boundary conditions for these DEs. The DEs have 6 singularities in the complex-z plane, z = 0, ±1, 2, 2(± √ 2−1), as shown in figure 5. For the physical region, the relevant poles are is a removable singularity. We will discuss in appendix B that this removable singularity does not affect the radius of convergence of asymptotic expansion at z = 1, although it can decrease the precision if we estimate values for z < 2( √ 2−1) from the asymptotic expansion at z = 1. The later problem has no impact if boundary conditions can be calculated to sufficient high precision, which is indeed the case as we will explain later. Therefore, virtual-correction MIs in regions 0 ∼ 1/4, 1/4 ∼ 3/4 and 3/4 ∼ 1 can be respectively estimated by the asymptotic expansions of MIs at z = 0, 1/2 and 1, where we introduce an expansion at a non-singular point z = 1/2 so that the combination of real corrections and virtual corrections can be expressed by a single piecewise function.
Finally, let us discuss how to obtain boundary conditions for DEs of virtual-correction MIs. We find that, if we choose boundary conditions at z → 1, calculation of these MIs either analytically or numerically to high precision is very hard. The method proposed in refs. [75,76] provides a way to calculate MIs numerically to very high precision at any non-singular point z, which we will explain in appendix D. With this method, we can not only provide boundary conditions for DEs, but also do a self-consistent check. To this purpose, we use this method to calculate MIs at two points, say z = z 1 and z 2 . With results at z = z 1 as boundary conditions, the DEs can give prediction for MIs at z = z 2 , and the later values can be compared with the values obtained by this method. In our work, we have done this self-consistent check, and find a perfect agreement.

Renormalization
Bare quantities of fields Ψ b and A b , coupling constant g sb , and heavy quark mass m Qb are related to corresponding renormalized ones by the renormalization constants δ 2 , δ 3 , δ g and δ m ,

JHEP04(2019)116
In this paper, we choose MS renormalization scheme for the coupling constant, and choose on-shell renormalization scheme for gluon field, heavy quark field and heavy quark mass. It is convenient to rescale the renormalization constants as following Summing over all counter terms, we obtain where |M LO | 2 is the squared amplitude at LO in α s , and dΦ Born is defined in eq. (2.14). This integral can be calculated easily. Besides, we need to renormalize the operator defining the FF. In MS scheme, the counter term gives where µ f is the factorization scale, d [1/8] LO (z) are given respectively in eq. (2.15) and eq. (2.16), and the Altarelli-Parisi splitting function P gg (z) is

JHEP04(2019)116
where b 0 is given below eq. (5.3), d LO (z) is given in eq. (2.18), and d [1/8] 3) The coefficients A k ij , B k j , C k ij can be evaluated numerically to very high precision, then analytical results can be obtained by fitting numerical results using PSLQ algorithm. For example, with 20-digit precision, we get In practice, however, numerical results with high precision will be sufficient. In appendix E, we present these coefficients up to j = 40 with 14 digits for each coefficient. With these numerical results, we can calculate d [1/8] NLO (z) to more than 12-digit precision for any value of z. To obtain results with 160-digit precision at any value of z, we attach an ancillary file with these coefficients calculated up to j = 530.

Numerical results
To see the effects of NLO corrections, we choose parameters the same as that in ref. [54], with m b = 4.75 GeV, N c = 3, n f = 4, and α s (µ r = 2m b ) = 0.181. In figure 7, we plot the curves of LO FFs and LO+NLO FFs with µ r = µ f = 2m b . To show color-singlet FFs and color-octet FFs in the same figure, we introduce overall factors c [1] = 6m 3 and c [8] = 96m 3 /5 for them, respectively. We find that our result of NLO FF of g → QQ( 1 S [1] 0 ) + X has some differences from that obtained in ref. [54], especially when z → 0. With µ r = µ f = 2m b , we also provide K-factors (the ratio of LO+NLO over LO) of some special values of z in table 1, where we find that K-factors are very significant for most values of z.
As shown in eq. (6.3), NLO FFs are negative and divergent at both z = 0 and z = 1, with leading behavior 1/z at z → 0 and ln 2 (1 − z) at z → 1. Thus total fragmenting probability obtained by integrating NLO FFs over z from 0 to 1 are divergent. The divergence at z → 0 is not a big problem because physical cross sections are obtained by convolving 0 ) and g → bb( 1 S [8] 0 ) at LO and NLO. The dotted line is for d [1] LO (z)×(6m 3 b ) or d [8] LO (z)×(96m 3 b /5), the solid line is for d [1] LO (z)+d [1] NLO (z) × (6m 3 b ) and the dashed line is for d [8] LO (z) + d [8] NLO (z) × (96m 3 b /5), with scale choices µ r = µ f = 2m b . The superscript [1] or [8] respectively denotes the color-singlet or color-octet states of bb. FFs with partonic hard parts, which behave as z n in the small z region with n usually larger than 4. 2 Therefore, the small z region has negligible contributions to physical cross sections. The divergence at z → 1 region, on the other hand, means that perturbative calculation in this region is not good. We leave the study of resummation of FFs in the z → 1 region for future works.

A IBP reduction with unregularized rapidity divergence
If all integrals are well regularized by dimensional regularization, IBP reduction relations should be independent of the infinitesimal imaginary parts iη, which means that coefficients of the relations are the same no matter a denominator is E j + iη or E j − iη. This is the reason why we ignore the infinitesimal imaginary parts when using IBP reduction.
However, in this paper we encounter some integrals that can not be regularized by dimensional regularization only. There is a MI in the calculation of real correction , (A.2) where we integrated out k − 1 , k − 2 and k + 2 , denoted k + 1 = (1−z)z 1 P + c , and did the replacement It is clear now that the integration over z 1 is divergent at z 1 = 0 and it can not be regularized by dimensional regularization. This divergence is usually called rapidity divergence, and it is in fact well-known that it cannot be regularized by dimensional regularization.
which does not change the integral because of δ functions in the definition of dΦ real . If we then simply replace these δ functions by propagator denominators using eq. (3.3) and perform IBP reduction of original expression by ignoring the infinitesimal imaginary parts, we find two equal loop integrals, Because they are equal to each other, we can choose either the former or the later as our MI. However, on the other hand, once we replace propagator denominators back to δ functions, the eq. (A.6) will vanish as it lacks of E 12 , while the eq. (A.5) will be changed to MI in eq. (A.1). Therefore, the final results are ambiguous.
To get unambiguous results both in the reduction step and in the calculation of MIs, we in principle need all involved integrals to be well regularized. We thus introduce an additional regulator besides spacetime dimension D = 4 − 2 , and take the limit of this new regulator to zero before take the limit of → 0. In this way, divergences that are regularized by dimensional regularization will not be affected. A possible choice of the new regulator is gluon mass in the phase space integration, which means we use instead of dΦ. Note that, although gluon mass should also be introduced in Feynman amplitudes to be self-consistent, it is easy to show that only the gluon masses in phase space integration have non-vanishing effect. With this regulator, we find all involved integrals in our calculation are well regularized, and thus the IBP reduction do not introduce any ambiguity. After the IBP reduction and then take the limit m g → 0 in any place as far as the operation does not result in unregularized integrals, m g still presents in four MIs besides the other four MIs obtained by changing E 4 to E 5 .
As an example, we calculate the first MI in eq. (A.8) in the limit of m g → 0. After we integrate out k − 1 , k − 2 , k + 2 , k 2⊥ and k 1⊥ sequentially, we get

JHEP04(2019)116
where t = (1 − z)/z. Because of dimensional regularization, only the region z 1 ∼ m 2 g survives in the limit of m g → 0. So we set z 1 = m 2 g y and take limit of m g → 0, then we get It tells us that, after introducing and then removing the gluon mass regulator, the MI in eq. (A.1) is eventually well regularized by dimensional regularization. We can similarly calculate the other three MIs in eq. (A.8), and find that they can be obtained from the first MI by multiplying factors 2 , 1 and z/(1 − z), respectively. Before describing how to apply the above method to our problem, let us first examine the following integral which is well regularized by dimensional regularization, 3 and thus we can calculate it numerically without introducing any other regulator. On the other hand, again without introducing any other regulator, we use IBP naively to reduce it to MIs. We find the reduced result is unique. All MIs obtained here except the one in eq. (A.1) are well regularized by dimensional regularization, which can be easily evaluated. Then if we replace the unregularized MI by eq. (A.10), we find the numerical result of eq. (A.11) agrees with the value calculated by applying IBP reduction. This test tells us two things. The first is that our gluon mass regulator can indeed give correct result. We thus take eq. (A.10) as the value of MIs defined in eq. (A.1), and similarly for other unregularized MIs. The second is that, if the original express is well regularized by dimensional regularization, using IBP naively may have no problem. Based on the above lessons, we divide our original express two parts. The first part is well regularized by dimensional regularization, which is then reduced to MIs by using IBP naively. We check this part numerically and find good agreement between results before and after the IBP reduction. As the second part is unregularized, we introduce the above gluon mass regulator before applying IBP. After inserting the values of MIs, we find the second part in our decomposition eventually vanishes.

B Removable singularities and their effects
In this work, we encounter some removable singularities. Some of them determine the convergence radius of asymptotic expansion at some points, and others only decrease the precision of higher order coefficients in the asymptotic expansion. In the following discussion, to be definite we discuss the case where there is a removable singularity at z = 1/2 and there is a non-removable singularity at z = 1. We will do asymptotic expansion at z = 0.
Let us first discuss the case where there are more than one analytical structure at z = 0, and the singularity at z = 1/2 is removable only after the summation of contributions from JHEP04(2019)116 all structures. Here is an example, where z = 1/2 is indeed a removable singularity. When we do the asymptotic expansion at z = 0, we get two series, where one comes from the analytical part and the other one comes from the part proportional to ln z. As z = 1/2 is a non-removable singularity of each of the two parts, the convergence radius of each of the series is 1/2. Although z = 1/2 becomes a removable singularity in the summation of two parts, the convergence radius of the asymptotic expansion at z = 0 is still 1/2. The reason is that we have no way to reorganize the two series to a single series so that it is convergent everywhere in 1/2 < |z| < 1. In our calculation, MIs in real corrections have this kind of removable singularity at z = 1/2, which determine convergence radius of asymptotic expansion at both z = 0 and z = 1. Now let us discuss the case where z = 1/2 is a removable singularity for each nonvanishing analytical structure at z = 0. A special case is that there is only one nonvanishing analytical structure, for example where 1/2 is a removable singularity while 1 is a branch point. If we denote we have a n = 2 n+1 ln 2 − based on which we can calculate the convergence radius: lim n→∞ an a n+1 = 1. We thus find that the singularity z = 1/2 does not affect the convergence radius at z = 0. However, we will show that this singularity has other effects. To this purpose, we note that g(z) satisfies the following DE, with initial condition g(0) = a 0 = 2 ln 2. The DE can generate the recursion relation which determines higher order coefficients in the expansion. However, when we solve DEs numerically, the initial condition can have only finite precision. If we denote the absolute error of a 0 as λ, the absolute error of a n calculated from the recursion relation is 2 n λ. At the point z = x, the contributed error from a n is (2x) n λ, which is much larger than λ if x > 1/2 and n is large. If we reduce this error by truncating the expansion to small n, then

JHEP04(2019)116
there will be a large systematic error at the order of x n . The best accuracy at z = x that one can obtain is to choose a truncation n so that (2x) n λ ∼ x n , which gives n ∼ log 2 λ −1 and x n ∼ λ − log 2 x . For example, for the point x = √ 2/2, the smallest absolute error that we can get is λ 1/2 , which is larger than the absolute error λ at x = 0.
In virtual corrections, if we do asymptotic expansion at z = 1, the removable singularity at z = 2( √ 2 − 1) belongs to the second type, and the convergence radius is determined by the singularity at z = 0. Let us denote 1 − 2( √ 2 − 1) = a −1 , then the best accuracy at z = 1−x that we can obtained is determined by (ax) n λ ∼ x n , which gives n ∼ log a λ −1 and x n ∼ λ − log a x . In this work, we want to estimate the value at z = 3/4 from the expansion at z = 1, which gives the best accuracy about λ 0.786 . If we need the accuracy to be about 10 −15 , we find λ ∼ 10 −19 and n ∼ 25, which means that we need initial condition for the expansion at z = 1 to have four more significant digits.

C Boundary conditions of MIs in real corrections
MIs in real corrections have the form where d ∈ {8, 9, 10, 11}, a, b, c ∈ {1, . . . , 7} and dΦ real is defined in eq. (3.1). We calculate most MIs in the limit z → 1 in this appendix, with the other MIs which are not regularized by dimensional regularization calculated in appendix A. From eq. (3.1), δ k 1 · n + k 2 · n − 1−z z P · n together with conditions k + 1 > 0 and k + 2 > 0 requires that k + 1 and k + 2 must be at the order of 1 − z when z → 1. Otherwise, if taking k + 1 (1 − z)P + as an example, the integral will be proportional to which equals 0 in dimensional regularization. Introducing the parametrization k + 1 = (1 − z)z 1 P + c and k + 2 = (1 − z)(1 − z 1 )P + c , dΦ real becomes

JHEP04(2019)116
In the limit of z → 1, E i given in eq. (3.2) becomê where λ = 1 − z. As λ → 0, each MI has at most four non-vanishing regions, To obtain boundary conditions, we only need to calculate the leading contribution in each region, which is proportional to λ n with n = −2, −3 or −4. The calculation is a little different depending on whether E 1 presents in eq. (C.1). Without E 1 , there is no cross term k 1⊥ · k 2⊥ in the limit λ → 0. We thus first rescale momenta by and then integrate out k 2⊥ . After that, the integration over k 1⊥ is very simple unless it has the form with a(z 1 ) = 0, 1. For this kind of integrals, we integrated out k 1⊥ after Feynman parametrization. Finally, the integration over z 1 and Feynman parameters can be calculated analytically with the help of sector decomposition refs. [80,81], which can isolate mixed divergences from parameter integrals. There are two widely used programs that can do the sector decomposition, SecDec [82][83][84][85] and FIESTA [86][87][88][89]. We use SecDec in this paper.
If E 1 presents, we first rescale momenta by

JHEP04(2019)116
and then do the the replacement k 1⊥ → k 1⊥ + k 2⊥ , which changesÊ 1 to k 2 1⊥ and moves the cross term k 1⊥ ·k 2⊥ to other denominators. To proceed, we introduce Feynman parametrization and integrate out k 2⊥ . Then the integration of k 1⊥ has the form which can be easily integrated out. Finally, integration of Feynman parameters can be worked out with the help of sector decomposition. All analytical results of MIs calculated here have been checked by numerical results computed by SecDec, and good agreement is found.

D Calculation of MIs in virtual corrections
MIs for virtual corrections in eq. (4.1) can be expressed as where F i are defined in eq. (4.2) with k 2 = 0, k − = k 2 ⊥ /(2k + ) and k + = (1 − z)P + /z. We apply the method proposed in refs. [75,76] to calculate these MIs at any regular point z = z 0 . To this purpose, we change F i to F i + iη for i = 1, 2 to obtain new MIs. We can set up DEs of the new MIs by first differentiating them with respect to η and then reducing the obtained expressions to the new MIs using IBP reduction. If we also know boundary conditions of the new MIs at a special value of η, we can solve the DEs numerically to obtain the new MIs at η = 0 + with very high precision, which are nothing but our desired old MIs.
The boundary that we choose is at η → ∞. To calculate the boundary conditions, we first perform Feynman parameterization and then shift l to remove cross terms. The obtained results are proportional to (D.2) where b, c, d are functions of z and the Feynman parameters x 1 , . . . , x n , and a is a function of z. As η → ∞, there are only two regions for this integral, The leading term of the first region gives which can be easily integrated out. The leading term of the second region gives

JHEP04(2019)116
the integrand of which is proportional to b n 1 −1+ after integrating out k ⊥ and l. Though b is a function of Feynman parameters and z, the dependence of z can be factorized out. So the integration over Feynman parameters can be easily performed.

E Coefficients
In this appendix, we give the coefficients defined in eq. (6.3). The coefficients of asymptotic expansion at z = 0 with different powers of ln(z) are shown respectively in tables 3-5. The coefficients of asymptotic expansion at z = 1/2 are shown in table 6. The coefficients of asymptotic expansion at z = 1 with different powers of ln(1 − z) are shown respectively in tables 7-10. To obtain results with 160-digit precision at any value of z, we attach an ancillary file with these coefficients calculated up to j = 530.
Note added. While this paper was being finalized, two related preprints appeared [90,91]. In ref. [90], the authors calculated NLO corrections for FF of g → QQ( 1 S [8] 0 ) + X using FKS subtraction method; while in ref. [91], the authors calculated NLO corrections for FFs of g → QQ( 1 S [1] 0 ) + X and g → QQ( 1 S [8] 0 ) + X using sector decomposition. Our highprecision results agree with K-factors obtained in these two works within their estimated errors.

JHEP04(2019)116
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.