Singlet and non-singlet three-loop massive form factors

We consider Quantum Chromodynamics with external vector, axial-vector, scalar and pseudo-scalar currents and compute three-loop corrections to the corresponding vertex function taking into account massive quarks. We consider all non-singlet contributions as well as those singlet contributions where the external current couples to a massive quark loop. We apply a semi-numerical method which is based on expansions around singular and regular kinematical points. They are matched at intermediate values of the squared partonic center-of-mass energy $s$ which allows to cover the whole kinematic range for negative and positive values of $s$. Our method permits a systematic increase of the precision by varying the expansion depth and the choice of the intermediate matching points. In our current set-up we have at least seven significant digits for the finite contribution of all form factors. We present our results as a combination of series expansions and interpolation functions which allows for a straightforward use in practical applications.


Introduction
Form factors are important building blocks for a number of processes. For example, they constitute the virtual corrections for lepton pair production via the Drell-Yan process, for Higgs boson production in gluon fusion and for Higgs boson decay into heavy quarks. For many processes there are precise experimental results such that higher-order corrections have to be included in the theory predictions. This is particularly true for QCD corrections, the topic of this paper.
There is a rapid increase in complexity when increasing the number of loops for a given process and often one has to rely on approximations valid in certain regions of phase space, in case one wants to perform analytic calculations. Alternatively, it is possible to use numerical methods, which, however, are less flexible in practical applications of the results.
In this paper we consider three-loop corrections to massive form factors in QCD. We use a method which leads to compact expansions of the bare three-loop expressions around kinematical points with high-precision numerical coefficients. Their numerical evaluation is fast and their combination covers the whole kinematic range. The counterterms are known in analytic form and allow for a flexible construction of finite expressions in different renormalization schemes.
We define the external vector, axial-vector, scalar and pseudo-scalar currents via j v µ =ψγ µ ψ , j a µ =ψγ µ γ 5 ψ , j s = mψψ , j p = imψγ 5 ψ . (1) The heavy quark mass m has been introduced in the scalar and pseudo-scalar currents for convenience such that all currents have vanishing anomalous dimensions [1].
Here the momentum q 1 (q 2 ) is incoming (outgoing) and q = q 1 − q 2 is the outgoing momentum at the current j δ with q 2 = s. The external quarks are on-shell, i.e. q 2 1 = q 2 2 = m 2 , and we have σ µν = i[γ µ , γ ν ]/2. We note that in all cases the colour structure is a simple Kronecker delta in the fundamental colour indices of the external quarks and not written out explicitly.
We distinguish singlet and non-singlet form factors. In the first case the external current couples to a closed quark loop which is connected to the quark line in the final state via gluons. In the non-singlet case the external current couples directly to the finalstate quarks. In this paper we consider the complete non-singlet contributions as well as those singlet contributions where the closed quark loop has the same mass as the external quarks. While we consider all four currents in the non-singlet case, we only consider the vector, scalar and pseudo-scalar currents in the singlet case. The proper treatment of the axial-vector contribution requires the singlet contributions with a massless closed quark loop, which is postponed to the future.
We define the perturbative expansion of the scalar form factors as where at lowest order we have F = 0. The three-loop non-singlet form factors can be decomposed into the ten colour factors where C F = T F (N 2 C − 1)/N C and C A = 2T F N C are the quadratic Casimir operators of the SU(N C ) gauge group in the fundamental and adjoint representation, respectively, n l is the number of massless quark flavors, and T F = 1/2. For convenience we introduce n h = 1 for closed quark loops which have the same mass as the external quarks.
The singlet form factors start at two-loop order. Due to Furry's theorem the vector form factor is non-zero only at three-loop order, but there are non-zero results for F s, (2) sing and F p, (2) sing . 1 For the three-loop singlet form factors we have the colour factors where (d abc ) 2 = (N 2 C − 1)(N 2 C − 4)/(16N C ) arises from Feynman diagrams in which the closed fermion loop is connected to the external quarks by three gluons, see Fig. 1(d). This is the only colour structure present in F Two-loop corrections to the vector form factors have been computed about 20 years ago, first in the context of QED [2,3] and later also for QCD [4] (see also Ref. [5] for the fermionic contributions). Several groups have provided cross checks and computed higher order terms in [6][7][8][9][10]. Higher order perturbative contributions in the high-energy limit of the form factors have been predicted in Refs. [6,8,11] using renormalization group equations. Two-loop axial-vector, scalar and pseudo-scalar contributions have been computed in Refs. [9,[12][13][14].
Analytic three-loop corrections to the form factors have first been considered in the large-N C limit [7,[15][16][17] where only planar integrals have to be computed. The first non-planar contributions appear in the light-fermion contributions which are available from Ref. [10]. In Ref. [18] all contributions with a closed heavy quark loop have been considered and around 2000 expansion terms around s = 0 have been computed. Let us also mention that all-order corrections to massive form factors in the large-β 0 limit have been considered in Ref. [19], where β 0 is the one-loop correction to the QCD beta function.
In Ref. [20] we computed first complete results for the vector form factors F v 1 and F v 2 taking into account all colour factors and covering the whole s range. In this work we present details of the computational method and extend the results to the axial-vector, scalar and pseudo-scalar currents. We present results for the non-singlet contributions, where the external current couples to the quarks in the final state. In this case it is possible to use anti-commuting γ 5 . Furthermore, we consider the singlet contributions, where the external current couples to a closed massive quark loop. For the treatment of γ 5 in the pseudo-scalar singlet case we follow Ref. [21]. Sample Feynman diagrams for the heavy-quark form factors are shown in Fig. 1.
The remainder of this paper is organized as follows: In Section 2 we describe the setup used for the computation of the amplitudes for the form factors. In Section 3 we describe our approach for the construction of the approximations of the master integrals. Section 4 is dedicated to the singlet form factors and the renormalization and infrared subtraction is discussed in Section 5. Our results are presented in Section 6. We conclude in Section 7. Supplementary material is relegated to the Appendix. In Appendix A we provide results for the three-loop on-shell integrals with special emphasis on the higher order terms in the expansion, which we needed as boundary conditions to determine the master integrals. Appendix B contains the collection of the plots which show the accuracy of the pole cancellation and in Appendix C we present semi-analytic expansions around s = 0, 4m 2 and ∞.

Reduction to master integrals and choice of basis
In this Section we provide details for the non-singlet form factors; the discussion for the singlet form factors is postponed to Section 4.
We generate the amplitudes with qgraf [39] and use q2e and exp [40][41][42] to rewrite the output to FORM [43] notation and map each diagram to a predefined integral family. In this way we can express the form factors as linear combinations of scalar Feynman integrals with twelve indices where nine correspond to the exponents of propagators and the remaining three to the exponents of irreducible numerators. In total we have 34 different integral families.
Before performing the reduction of the Feynman integrals contributing to the form factors, we improve the basis by getting rid of those denominators in the coefficients of the master integrals which are multivariate polynomials. Since the ratio s/m 2 is the only kinematic variable, it is possible to find a basis such that the denominators completely factorize into univariate polynomials of either s/m 2 or d [44,45]. Moreover, it turns out that we can choose a basis such that all polynomials of s/m 2 in the denominators are linear polynomials raised to some power. We call this basis good in the following. To construct this basis we first reduce all integrals up to the top-level sector and with up to two dots for every integral family individually with the help of Kira [46,47] employing Fermat [48]. As initial basis we simply take the default master integrals found with the integral ordering 6 of Kira, i.e. the sectors are primarily ordered by the number of lines and dots are preferred over scalar products. These reduction tables then serve as input to search for a good basis for every family with the help of an improved version of ImproveMaster.m developed in Ref. [44]. 2 The construction of a good basis takes about three hours for the most expensive families and almost all of that runtime is spent for the reduction.
The actual reduction of the integrals for the form factors is also performed with Kira employing Fermat. First, we reduce the integrals for every family individually to the good basis of this family. The most expensive families run for about a week on eight cores and require about 200 GiB of memory. The resulting reduction tables, especially the denominators, are simplified by using Mathematica to factorize the expressions which is not done by Kira and Fermat. Secondly, we employ symmetries between the families to reduce the number of master integrals. To this end we sort the families by ascending number of master integrals and use Kira to reduce the good master integrals for every family to those of easier families if possible. This reduces the number of master integrals from 3131 to 422. 136 of these master integrals are needed for the fermionic contributions with at least one closed massive or massless quark loop. This step takes about one day and can be done in parallel to the individual reductions of all families since the non-minimal good basis is already known.
We set up the differential equations for the master integrals by differentiating the master integrals with respect to s/m 2 with the help of LiteRed [49,50] and reducing the resulting integrals again to master integrals with Kira. Also this step is performed on a per-family basis and the already known symmetry tables are only inserted in the very last step. The construction of the differential equations takes at most a few hours per family. The only complication at this point is that a few blocks on the diagonal of the system of differential equations contain unfavorable pole structures that prohibit a solution of the master integrals as series in . However, this problem can be solved by searching the available reduction tables for integrals in the affected sectors with a pole in front of one of the master integrals of the same sector. Inverting the relation, i.e. switching integral and master integral, thus introduces a positive power of and improves on the pole structure [51]. By repeating this step all sectors can be made solvable. For our problem this can be achieved without spoiling the good properties of the basis. Let us emphasize that we explicitly did not try to get rid of all poles in our reduction tables because this would most certainly spoil the good properties.

Method
In order to (numerically) solve the master integrals we apply the method of analytical expansions and numerical matching introduced in Ref. [52]. Let us briefly summarize the basic idea: (a) After establishing a system of differential equations for the master integrals, we calculate initial values at a kinematic point where the integrals simplify.
(b) We construct an analytic series expansion around this kinematic point by inserting a suitable ansatz into the system of differential equations, re-expanding in and around the special kinematic point, establishing a system of linear equations for the expansion coefficients and solving it in terms of a small set of boundary coefficients. We use the information from step (a) to fix these coefficients analytically.
(c) Construct an analytic expansion around a neighboring kinematic point. Here we cannot fix the boundary constants analytically, but we can evaluate both expansions at a kinematic point where the radii of convergence of both expansions overlap and use these values as numerical initial conditions.
(d) Repeat this procedure until the whole kinematic interval of interest is mapped out with partially overlapping series expansions.
In the following we give details on how we applied this method to the calculation of the massive form factors. We concentrate on the non-singlet form factors; details for the singlet form factors are provided in Section 4.

(a) Calculation of initial values
For our initial expansion we choose the value s = 0. At this point the master integrals simplify to on-shell propagator integrals which are well studied in the literature [24,53,54]. However, one obstacle in obtaining the initial values for our system of differential equations is the appearance of high spurious poles in the dimensional regulator in the physical amplitude and the differential equation. This requires the knowledge of some master integrals beyond the order given in Ref. [24]. The same problem was also encountered in Ref. [18], where a subset of integrals was already obtained to sufficiently high order in . We extended the calculation of the remaining integrals to the order in needed in our approach. Details on the calculation and the results are given in Appendix A. We note that due to the spurious poles of O(1/ 8 ) we need some integrals up to transcendental weight 9, however, in the physical amplitude all contributions with weight higher than 5 cancel.

(b) Construction of series expansions
In order to construct series expansions around specific values of s we insert a suitable ansatz into the differential equation, re-expand in and around the chosen value of s and construct a system of linear equations for the expansion coefficients of the ansatz by equating coefficients. The suitable ansatz can be found from physical considerations. We only have three kinematic points with non-analytical expansions. They correspond to the two-particle threshold s = 4m 2 , the four-particle threshold s = 16m 2 and the high-energy expansion s → ±∞. For these three expansion points we have to use a power-log ansatz.
Since the threshold expansions are governed by the particle velocities, we use the variables z = √ 4 − s and z = √ 16 − s for the two-and four-particle thresholds respectively. The ansatz for a master integral I n is then given by For the present three-loop problem we have o = 3 and we choose k = 10 as a conservative lower bound. In principle both lower values depend on the master integral and can be determined beforehand. In practice, we choose them uniformly for all master integrals for simplicity. The high-energy expansion can be performed in the variable z = 1/(−s), but one has to allow for Sudakov-like double logarithms in the ansatz. The upper summation variable of the sum over m in Eq. (6) has therefore to be modified to i + 2o. All other expansions are given by simple Taylor expansions in the variable z = s 0 − s, with s 0 the value we want to expand around. Thus, the sum over m can be dropped and k = 0.
The resulting system of linear equations, which does not contain any variable anymore, is solved by Kira together with Fermat or FireFly [55,56]. For the use with FireFly we use a special version of Kira which allows the efficient use of finite field methods and rational reconstruction even for systems without variables. We make sure to prefer boundary constants corresponding to low powers of z and in the ansatz when solving the system.

(c) Numerical matching of neighboring expansions
Once we found the symbolic expansion around s = 0 we can use the previously calculated initial values to fix the remaining boundary constants. This results in an analytic expansion of the master integrals around s = 0, where we calculated 76 terms in the expansion.
In a next step we construct a symbolic expansion around s = 1. For s = 1 the boundary constants cannot be easily fixed analytically. However, since the radii of convergence for the expansion around s = 0 and s = 1 overlap, we can use the numerical evaluation of the expansion around s = 0 at s = 1/2 to obtain numerical boundary conditions for the expansion around s = 1.
In practice, we equate the expansions around s = 0 and s = 1 after setting s = 1/2 to obtain a linear system of equations for the remaining boundary constants in the s = 1 expansion. For analytical values at s = 1/2 this system has a unique solution. However, due to the finite numerical precision this is not the case anymore and we have to be careful in solving the system. We proceed in the following way: • In order to get numerically more stable results we rationalize the numbers in the system with a preset accuracy.
• We sort the equations with ascending number of appearing boundary constants and start by solving the equations one by one and insert the solutions back into the remaining equations. After inserting we rationalize again to the preset accuracy in order to avoid bigger and bigger rational numbers. Furthermore, we set numbers with absolute values smaller than a threshold to zero. This is important, since sometimes boundary conditions in the system are multiplied by very small numbers of the order of the preset accuracy. If we accidentally solve for these constants we introduce large values into the system of equations which render the solution of the system unstable.
• Sometimes after reinserting relations, which have been found before, back into the system, equations with no boundary constant remain. In an analytic matching these numbers would be identically zero. We use the absolute value of these remaining numbers to judge the accuracy of the matching.
To evaluate the stability of the final solution we solve the systems twice: once with rationalized coefficients with an accuracy of 500 digits and setting all numbers smaller than 10 −100 to zero and once with rationalized coefficients with an accuracy of 100 digits and setting all numbers smaller than 10 −50 to zero. We find agreement of both runs with the expected accuracy of ∼ 100 digits.

(d) Mapping out the kinematics
In previous calculations of the form factors the variable x defined by has been used. This choice of variable maps the special points of the physical amplitude corresponding to low energy (s = 0), production threshold of two heavy quarks (s = 4m 2 ) and the high-energy limit (s = ±∞) to the points x = 1, 0, −1 and, up to two-loop order, allows to express the final result in terms of harmonic polylogarithms. However, in our approach it turns out that the variable x is not a good choice. One problem we encounter is that the production threshold of four heavy quarks, which starts contributing at threeloop order, is mapped to x = 4 √ 3 − 7 ≈ −0.072. As a consequence the expansion around x = 0 has a very small radius of convergence and large numerical coefficients appear in the expansions, making the numerical solution less stable. In addition, we also have to construct expansions around the value of the new threshold. In principle this can be achieved by defining the square root as a symbol before passing the system of equations to Kira. However, this complicates the solution of the system considerably. We therefore choose to do the expansions in the variable s.
To obtain the results for the non-singlet 3 form factors over the whole real axis we construct expansions at the values of starting from our initial expansion at s = 0. 4 The numerical matching between two neighboring expansions is always done close to the center of the interval. Some of these values are chosen to correspond to additional, unphysical poles in the differential equation which we observe at s m 2 = −4, −2, −1, −1/2, 1/2, 1, 2, 3, 16/3 .
However, it turns out that these additional poles are all spurious and do not spoil the convergence of the series expansions. Thus, expansions around s/m 2 = 4, 16, ∞ are in principle sufficient to construct the form factors at arbitrary values of s. However, we find that this approach suffers from a slow convergence of the individual series expansions close to the singular points, which means that very deep expansions are required to get decent accuracy at the matching points. Since very deep expansions are expensive in our approach we trade the expansion depth with more expansions at intermediate points with a default expansion depth of 51 terms. These intermediate points are all Taylor expansions and therefore quite inexpensive to calculate.
For series expansions close to a singular point we employ Möbius transformations, which have already been discussed in Ref. [57]. Assume, we want to expand around the point x k and there are singular points of the differential equations at x k−1 and x k+1 with x k−1 < x k < x k+1 . Naively the radius of convergence is limited by the distance to the closer singular point. However, the variable transformation maps the points x k−1 , x k , x k+1 to −1, 0, 1. The radius of convergence of the series expansion is therefore extended into the direction of the farther singularity although the convergence at the boundaries can be quite slow. For example we construct the expansion around s/m 2 = 7/2 in the difference x = 7/2m 2 − s, but afterwards reexpress the result in While the expansion in the variable x converges for s ∈ (3, 4), the convergence for the expansion in y is increased to s ∈ (−∞, 4), making the evaluation at the matching point s = 15/4 much more precise.
Let us comment on the resources needed to construct the series expansions in our approach. To set up the system of linear equations for the construction of the series expansions we use Mathematica and can trivially parallelize this step over all master integrals. For the most complicated master integrals this step takes around 1 h on a single core. The solution of the resulting system for a simple Taylor expansion over one prime field takes around 6 h and requires about 50 GiB of memory. In total around 50 evaluations are needed which can all run in parallel. The power-log expansions are more expensive to compute. Due to the Sudakov-like double logarithms the high-energy expansion is the most involved. Here the solution of the linear system of equations over one prime field takes around 10 d with a memory requirement of 250 GiB, but again the evaluation over 50 prime fields is enough to reconstruct the full result.
The idea to utilize the associated system of differential equations to obtain expansions or numerical evaluations of the master integrals has received a lot of attention especially in the recent past [57][58][59][60][61][62][63], and there are even public packages available which implement different algorithms [64][65][66][67]. However, most of them have restrictions on the structure of the system of differential equations, e.g. it has to be triangular for SeaSyde [67] or Fuchsian without resonances for DESS [64], or on the solution, e.g. only formal power series [60]. Some methods also lack the proof that they can handle complicated physical problems with a few hundred master integrals. Here, we show that the method laid out above can be used to compute a nontrivial quantity, namely the massive form factors at three-loop order in QCD, and obtain precise numerical results over the whole parameter space.

Singlet contributions to the form factors
In this Section we discuss the calculation of the singlet contribution to the three-loop massive QCD form factors, i.e. contributions of the type shown in Fig. 1(d) where the current couples to a closed heavy-quark loop. We present results for the vector, scalar and pseudo-scalar case. The axial-vector contribution is anomalous and needs the inclusion of singlet diagrams where the current couples to a closed light-quark loop. Since many of the technical details of the calculation closely follow the non-singlet case, we only discuss them briefly. We put special emphasis on the computation of the boundary conditions in the limit s → 0 since the calculation is different from the one in the non-singlet case.
Results for the singlet contribution are shown in Section 6 together with the non-singlet case.
To generate the amplitude in the singlet case we need 17 different integral families. The resulting list of scalar integrals is again reduced with Kira and Fermat where we find 316 master integrals after symmetrization over all families. In this case we not only make sure to reduce to a good basis, but in addition reduce the number of spurious poles in . This is achieved with an improved version of ImproveMasters.m. Afterwards we again establish a closed system of differential equations with the help of LiteRed with subsequent reductions using Kira. Overall the reduction of the singlet diagrams is significantly simpler than the one in the non-singlet case.
In principle the whole method discussed in Sec. 3 can be directly applied to the singlet diagrams. However, the calculation of the boundary conditions at s = 0 is more involved, since the singlet diagrams can possess massless cuts. Hence, a simple Taylor expansion around the limit s = 0 is not sufficient, but we have to perform an asymptotic expansion. We realize the asymptotic expansion by the method of regions [68,69], where the relevant regions can be found with the program asy.m [70].
Applying asy.m to the singlet master integrals we find that there are regions with three different scalings: with y = −s/m 2 . Region 1 corresponds to the hard region, where we can naively expand the integrands around s = 0 and apply the same procedure as for the non-singlet diagrams to obtain the boundary conditions.
The expansions in the two soft regions 2 and 3 do not permit a straightforward diagrammatic expansion. In order to calculate the boundary conditions in these regions we perform the expansion on the level of the integral representation in terms of α parameters and use direct integration methods to solve the occurring integrals. Let us exemplify the calculation of the boundary conditions for the two soft regions for the most involved master integral J we encounter. The corresponding graph is shown in Fig. 2.
The Symanzik polynomials U and F for this diagram are given by The graph corresponding to master integral J. Full (dashed) lines correspond to massive (massless) scalar propagators. The dot represents the influx of the external momentum q.
and the master integral can be written as With the help of asy.m we find the possible scalings of the α-parameters to be {0, 0, 0, 0, 0, 0} or {−1, −2, −1, −2, 0, −2}. The first scaling corresponds to the hard region, which we can calculate with the same method as in the non-singlet case. We obtain The second scaling corresponds to region 3, which can be seen by applying the scaling relations to the integrand in Eq. (14). In this region the U and F polynomials reduce to To obtain our boundary condition in region 3 we therefore have to solve the integral The α parameters α 1 , α 2 and α 5 can be integrated using the formulae together with suitable rescalings of the integration variables. In the end we are left with .
The remaining integrals cannot be further reduced by applying the formulae in Eqs. (19) and (20). We find it convenient to integrate Eq. (21) with the help of the program HyperInt [71], since the divergent pieces for → 0 are already factored out into the Γ functions and the integrand can be simply expanded in . After performing the integrals the result is given by Region 2 does not contribute to the master integral J, so we have J 2 = 0. The complete boundary condition can be obtained by adding up the contributions of all regions: Note that for regions 2 and 3 for all other boundary conditions the formulae in Eqs. (19) and (20) are sufficient to obtain the boundary in terms of Γ functions and therefore exact in .
After all boundary conditions for all three regions are found we again obtain analytic series expansions around s = 0. In the singlet case we calculate symbolic expansions at the values of and match numerically again starting from our initial expansion at s = 0. The singular points at s/m 2 = 4, 16, ∞ are shared among the singlet and non-singlet case. For the singlet diagrams also the point s = 0 is a singular point and we need to perform a powerlog expansion. Let us mention that in the singlet case the differential equation has the additional unphysical poles

Ultraviolet renormalization and infrared subtraction
The form factors develop both ultraviolet (UV) and infrared (IR) divergences. To account for the UV divergences we have to renormalize the quark mass, the strong coupling constant and the wave function of the external quarks. We denote the UV-renormalized form factors by F UVren which still have poles in of IR nature. They are taken care of with the help of an IR renormalization factor in minimal subtraction, Z, which is constructed from the QCD beta function and the cusp anomalous dimension (see, e.g., Ref. [7]). It is convenient to consider log Z which is given by with cusp are the expansion coefficients of the cusp anomalous dimension defined through The cusp anomalous dimension in QCD is available to three loops from Refs. [72][73][74][75].
Using Z, we can define a finite form factor via Spelling out this equations leads to where F (1+2+CT) contains the exact one-and two-loop expressions (including higher orders in ) and the three-loop counterterm contribution. The quantities Z and F (1+2+CT) are expressed in terms of HPLs, which depend on the quantity is the bare three-loop contribution, the main result of this paper. We construct numerical approximations valid in the whole s/m 2 range.
The one-and two-loop results for the form factor are available to order 2 and 1 , respectively, see, e.g., Ref. [15], which is necessary to obtain the finite contributions at order α 3 s . One has to take into account the following counterterms, which are well established in the literature, see, e.g., Refs. [76][77][78][79]: • α s in the MS scheme up to two loops. This is a multiplicative renormalization.
Since the one-and two-loop form factors and the counterterms are ξ independent also the induced three-loop terms are ξ independent. We express our final results in terms of α (n l ) s , the strong coupling constant with n l active quark flavours.
• Wave-function renormalization in the on-shell scheme to three loops. This is a multiplicative renormalization. Note that Z OS 2 has a ξ dependence at three loops in the colour factors C F C 2 A and C F C A T F n h . • Quark-mass renormalization in the on-shell scheme to two loops. Since we have for the external (anti)quark momenta q 2 1 = q 2 2 = m 2 , a multiplicative renormalization in the bare one-and two-loop expressions is not possible. Instead one has to generate and compute one-and two-loop counterterm diagrams with explicit mass insertions. While the mass counterterms themselves are ξ independent, linear and quadratic ξ terms are generated from the mass insertions nonetheless.
• The anomalous dimensions of the vector and axial-vector currents vanish and thus no renormalization is needed. The anomalous dimensions of the scalar and pseudoscalar currents only vanish because of the additional factor m in Eq. (1), which has to be renormalized. We choose to renormalize this factor in the MS scheme.
The MS renormalization constants only contain pole parts whereas for the on-shell quantities also higher order coefficients are needed since the one-and two-loop form factors develop 1/ and 1/ 2 poles, respectively.
We use the results from Ref. [15] and obtain exact results for the three-loop counterterm contributions. We expand the expression for s → 0 and check the cancellation in this limit analytically.
In principle we have to expand the counterterm contributions also around all other values of s for which we perform expansions. However, this can be quite tedious since we have to expand the iterated integrals around quite involved arguments. Instead we check the pole cancellation numerically using the approximated bare three-loop form factors, the exact results for the counterterm contributions and the exact result for the cusp anomalous dimension. We observe that the poles cancel with a relative precision of at least nine digits in the whole s range. This is discussed in more detail in Subsection 6.1.
Note that the singlet contributions start at two loops. Apart from that the renormalization proceeds in the same way as for the non-singlet case. Since the singlet contributions to the vector form factors vanish at two loops due to Furry's theorem, its three-loop amplitude is already finite and does not need to be renormalized.

Results
In this section we present our results for the three-loop corrections of the massive form factors. First, we discuss how we estimate the uncertainty of our numerical results and which cross checks we have performed to validate them. We then discuss the static and high-energy limits as well as the two-particle threshold and present the leading terms of these expansions. Furthermore, we comment on the behaviour of the four-particle threshold. Finally, we show plots for the form factors over the whole kinematic range. This is accompanied by a numeric package to evaluate the form factors as the main result of this paper. Unless stated otherwise we choose for the renormalization scale µ 2 = m 2 .

Estimation of the accuracy and cross checks
We estimate the precision of our result from the numerical pole cancellations of the renormalized and infrared-subtracted form factors: At each random sample point, for every colour factor, and for every order in , we add the numerical bare results and the numerical evaluations of the counterterms as well as Z defined in Eq. (29) and divide by the absolute value of the counterterms and Z: This corresponds to the precision of the pole terms. Real and imaginary parts are checked separately.
For illustration we show the cancellations for the three non-fermionic colour structures Fig. 3. Figures for the non-fermionic colour structures of the five remaining form factors can be found in Appendix B. Since the precision of the fermionic colour factors is much better, we refrain from showing the associated plots. In general we observe a progression in the orders of , i.e. the 1/ −3 poles cancel with the highest precision and we lose some digits with every higher order. Sometimes, this general progression is violated, usually if the value of the colour factor changes sign and crosses  zero. These zero crossing are visible in the plots, because the precision slowly decreases and then slowly increases again, see for example the region around s ≈ −20m 2 for the real part of the colour factor C 3 Fig. 3(e).
Let us now analyze the four physical regions s < 0, 0 ≤ s < 4m 2 , 4m 2 ≤ s < 16m 2 , and 16m 2 ≤ s separately in more detail. In each region we provide the minimal precision over all colour factors and form factors. Since the form factors are very small close to zero crossings, we also provide the minimal precision when removing all points for which the size of the coefficient is smaller than 2.5 % of the average in the region. In addition, we discard all points close to the Coulomb singularity, i.e. 3.95m 2 ≤ s ≤ 4.05m 2 .
The poles generally cancel with at least 15 digits for the 1/ −3 poles, at least 10 digits for the 1/ −2 poles, and at least 10 digits for the 1/ poles. Removing the points close to the zero crossings this improves to 15, 12, and 11 digits, respectively.
2. The region 0 ≤ s < 4m 2 is the most precise one: For the 1000 random sample points, the poles cancel with at least 17, 14, 12 digits for the 1/ −3 , 1/ −2 , 1/ poles, respectively. Removing the points close to zero crossings with the threshold chosen above does not improve the precision. However, these reported worst pole cancellations all belong to the power-log expansion around the two-particle threshold at s = 4m 2 . For 0 ≤ s < 3.75m 2 our precision is well beyond 20 digits as can be seen in the figures.
3. The region 4m 2 ≤ s < 16m 2 between the thresholds at s = 4m 2 and s = 16m 2 is least precise: For the 1000 random sample points, the poles of the real part only cancel with at least 11, 11, 8 digits for the 1/ −3 , 1/ −2 , 1/ poles, respectively, and the poles for the imaginary parts with 12, 6, 6 digits. Removing the points close to zero crossings mildly improves the precision of the real part to 13, 12, 9 digits. The imaginary part on the other hand significantly improves to 14, 11, 9 digits.
4. The region 16m 2 ≤ s becomes more precise again, since it is matched from +∞. For the 250 random sample points between 16m 2 ≤ s ≤ 100m 2 , the poles cancel with at least 14, 13, 8 digits for real parts of the 1/ −3 , 1/ −2 , 1/ poles, respectively, and with 15, 11, 8 digits for the imaginary parts. This improves to 14, 13, 10 digits for the real part and to 15, 11, 9 digits for the imaginary part when removing the points close to zero crossings.
Extrapolating these numbers to the finite terms, we expect that our result is correct up to at least 7 digits away from the zero crossings, with a much better precision for most colour factors and form factors over most parts of the real axis.
We have performed the calculation of the form factors for general QCD gauge parameter ξ and have checked that ξ cancels in the renormalized form factors. Note that the mass counterterm contributions depend on ξ which cancels against the bare three-loop expressions. We have checked the cancellation numerically and observe that the coefficient in front of ξ is of order 10 −18 or smaller in most of the phase space.
After specifying to the large-N C limit via C F → N C /2 and C A → N C we can compare the N 3 C terms against the exact results from Ref. [7,10]. In this limit only about 90 planar master integrals contribute and we observe a significantly increased precision of our result. In fact, in the whole s/m 2 region we can reproduce the exact result with at least 14 digits. with the exception very close to the singularity at s = 4. For example, for s = 3.9m 2 and s = 4.1m 2 we have an agreement of about 12 digits.
We observe similar results for the light-fermion colour factors C 2 F n l , C A C F n l , C F n h n l and C F n 2 l , which we compare against the exact results from Ref. [10], and the n 2 h terms where the exact expressions can be found in Ref. [18].
For the singlet contribution there are no 1/ 2 and 1/ 3 poles from the counterterms. For the 1/ pole we have a much higher precision as can be seen in Fig. 4 for the colour sing . We observe that the poles cancel with at least 20 digits, and with several more digits for most regions of the phase space. The cancellation for the other two colour factors with a second fermion loop is even more precise. Due to the high precision we refrain from showing plots for these two colour factors and F p,f, (3) sing . As a conservative estimate we claim a precision of at least 10 digits for the 0 terms of the singlet form factors. Since contributions to F v sing only start at three-loop order, we cannot check the pole cancellation in this case. However, we do not observe any difference in the behavior of the master integrals which only contribute in this case and, thus, assume the same precision.

Analytic and numeric expansions
Next we present expansions around the special kinematic points s = 0 and s = ±∞. For s = 4m 2 we construct expansions for the cross sections and decay rates, respectively. Such expansions around s = 0, 4m 2 and ∞ may serve as input for approximation procedures as those based on Padé approximations, see, e.g., Refs. [80][81][82].
6.2.1 Static limit: s → 0 In the static limit we construct an analytic expansion including 5 s 66 from the boundary values at s = 0. We restrict ourselves again to the five colour structures which are not known analytically. All colour factors as well as higher orders in the expansion are available in the ancillary file accompanying this paper [83]. For illustration we show the results for F v 1 and F s in the main part of the paper and relegate the remaining four form factors to Appendix C. For the non-singlet form factors the expansions up to s/m 2 are given by where l 2 = log(2), a 4 = Li 4 (1/2) and ζ n is Riemann's zeta function evaluated at n. For s = 0 our results for F v 2 (0) and F a 1 (0) agree with Refs. [84] and [85], respectively. Note that for s = 0 the cusp anomalous dimension in Eq. (26) vanishes and we have Z = 1.
The n h -singlet contribution to the scalar form factor reads with l √ −s/m = log( √ −s/m). This logarithm as well as the expansion in powers of √ −s/m, instead of s/m 2 as for the non-singlet contributions show above, originate from the massless cuts through singlet diagrams discussed in Section 4. The results for the vector and pseudo-scalar form factors can be found in Appendix C.

High-energy expansion: s → −∞
Also for the high-energy expansion we focus on the colour factors C 3 F T F n h , and C F C A T F n h and refer to the literature [10,15] for the remaining fermionic contributions. The high-energy expansions of the non-singlet form factors up to m 2 /s read F F a,f,(3) l and n 2 h terms, (40) with l s = log(m 2 /(−s − iδ)).
The leading logarithmic contributions of order α n s log 2n (m 2 /s) are given by the Sudakov exponent [86,87] , F s,f, (3) , and F p,f,(3) . Our numerical results (shown above with lower precision) are sufficiently precise to reconstruct the analytic coefficient (41) Similarly, we can reconstruct the analytic coefficients for the leading logarithms of the first mass corrections and find The latter agree with Refs. [88][89][90] where the results in Eq. (42) have been obtained using an involved asymptotic expansion of the three-loop vertex diagrams. Moreover, we confirm that there are only subleading contributions from the non-singlet diagrams to the remaining form factors. While This result disagrees with Ref. [91]. However, we can make both results agree by modifying Eq. (2.14) in Ref. [91] The correctness of our result has been confirmed by the authors of Ref. [91].
Finally, we show the reconstructed analytic coefficients for the remaining leading and first subleading logarithms for the first two terms in the high-energy expansion for all currents: Apart form the leading and subleading logarithms discussed above, our approach provides the whole tower of logarithms and also higher order contributions in m 2 /(−s). We estimate the accuracy of the non-logarithmic terms in Eqs. (35)-(40) to ten digits. For the subleading terms the accuracy decreases. Note, however, that we use the s → ∞ expansion only for |s/m 2 | 45 and that 1/45 3 ≈ O(10 −5 ).
For the scalar and pseudo-scalar singlet contributions we obtain the following results for the leading logarithmic contributions of the power-suppressed term which is in agreement with Ref. [90].

Threshold expansions: s → 4m 2 and s → 16m 2
Let us next discuss the two-and four-particle thresholds at s = 4m 2 and s = 16m 2 . Close to the two-particle threshold F 1 develops the famous Coulomb singularity with negative powers in the velocity of the produced quarks, β = 1 − 4m 2 /s, up to third order multiplied by log(β) terms. In this limit real radiation is suppressed by three powers of β and it is thus possible to construct physical quantities from the square of the form factors. For the four currents under consideration we follow Ref. [15] and define These quantities form building blocks for, e.g., cross sections of heavy quark production in electron-positron annihilation or decay rates for scalar or pseudo-scalar Higgs bosons (see also Ref. [15]). For reference we provide the (exact) leading order results which are given by where we adapt the notation from Eq. (3). We parametrize the QCD corrections to R δ with the quantities ∆ δ,(i) which we introduce as with K v = 3/2, K a = K s = K p = 1, n v = n p = 1 and n a = n s = 3. For convenience we set µ = m. In contrast to Ref. [15] we present results parametrized in terms of α (n l ) s (and not α (n l +1) s ). Furthermore, for the scalar and pseudo-scalar current we keep the factor m in the definition of the currents (see Eq. (1)) in the MS scheme and refrain from transformation to the on-shell scheme. Note that this is the natural choice for Higgs decays where the factor m takes over the role of the Yukawa couplings.
The three leading terms in β for the four currents read with l 2β = log(2β). Note that the fermionic contributions are suppressed by one additional power of β and, thus, we can show the term β 0 for them in contrast to the non-fermionic contributions. The non-fermionic part of ∆ v, (3) has already been shown in Ref. [20]. In Eqs. (50) to (53) we only show five digits for each coefficient, however, our results for ∆ δ, (3) contain more significant digits. For example, in the vector and pseudo-scalar case our numerical results reproduce the analytic expressions from Ref. [92] (see also Refs. [82,93]) with at least 13 digits accuracy. The light-fermion contributions can be compared with the analytic results of Ref. [15] and agreement is found for 19 digits. Similarly, after specifying to the large-N C limit we can reproduce the first 14 digits of Ref. [15] for all four currents.
The four-particle thresholds are much less pronounced in the results for the form factors.
For the individual master integrals we observe behaviours of the form (β 4 ) k × log l (β 4 ) with k = 3 and l = 5 where β 4 = 1 − 16m 2 /s. However, all form factors have a smooth behaviour for s → 16m 2 . In fact, in all cases there are no log(β 4 ) terms in our expansions for β 4 → 0. Furthermore, we observe the first non-analytic terms, i.e. terms where β 4 is raised to an odd power, at order (β 4 ) 7 for the axial-vector and scalar form factor and (β 4 ) 9 in the vector and pseudo-scalar case. This statement is true both for the non-singlet and singlet form factors. Such a high suppression can partly be explained by the fact that the massive four-particle phase-space, which is one of our master integrals, already provides a factor (β 4 ) 7 . Due to the more divergent behaviour of the master integrals it is nevertheless necessary to perform a careful matching to s = 16m 2 , both from above and below. This means we have to go quite close to s = 16m 2 using Taylor expansions around regular points (in our case we chose s = 15m 2 and s = 17m 2 ). Furthermore, we constructed 100 terms for the expansion in β 4 .

Numerical results for finite three-loop form factors
For illustration we show in Figs. 5, 6 and 7 the finite non-singlet form factors (see Eq. (30)) for negative s, for 0 < s < 4m 2 and above the two-particle threshold, respectively. Only in the latter case the imaginary parts are different from zero. We restrict ourselves to the nonfermionic colour factors and the contributions containing a closed heavy quark loop. In total we present results for the colour factors C 3 F v, f, (3)  F a, f, (3) F a, f, F s, f, (3) F p, f, (3)  F v, f, (3) F F a, f, (3) F a, f, F s, f, (3) F p, f, (3)  F v, f, (3) F F a, f, (3) F a, f, F s, f, (3) F p, f, (3) for the scalar and pseudo-scalar case they have a similar order of magnitude as the other colour structures.
In Fig. 6 one can clearly see the Coulomb singularities for the non-fermionic contributions close to s = 4m 2 . In the n h contributions the closed heavy-quark loop regularizes the 1/β behaviour and leads to a finite limit for s → 4m 2 , see also Section 6.2.3. It is interesting to note that for s → 0 we have F 1,sing (0) = const whereas the scalar and pseudo-scalar currents behave as log 2 (−s/m 2 ). The logarithms, which are present in higher expansion terms for all form factors and colour factors, are the reason for the imaginary parts for s > 0.
The behaviour around s = 4 is smoother than for the non-singlet form factors. In particular, there are no 1/β 3 or 1/β 2 singularities. The two vector form factors have a smooth limit and only the colour factor C 2 F T F n h of the scalar and pseudo-scalar form factors develop 1/β and log(β)/β singularities, both for the real and imaginary parts. The other three colour factors have a finite limit for s → 4.
In the high-energy region all form factors vanish except F

Conclusions
In this paper we compute three-loop corrections to massive quark form factors with external vector, axial-vector, scalar and pseudo-scalar currents and provide precise numerical results in the whole s/m 2 range. We apply the method developed in Ref. [52] to obtain expansions around regular and singular points. It is based on differential equations which are used to construct the expansions. Two neighboring expansions are numerically matched at an intermediate value of s.   We consider both the non-singlet and the singlet contributions, with the restriction that the external currents couple to a closed massive quark loop in the latter case . The main difference between the two contributions is the computation of the boundary conditions at s = 0. While we have simple Taylor expansions in the non-singlet case, it is necessary to perform an asymptotic expansion for the singlet contributions. Our expansions around s = 0 are analytic; the expansions around the other s values have precise numeric coefficients. In some cases the precision is sufficient to reconstruct coefficients analytically as, e.g., for leading and sub-leading logarithmic contributions in the high-energy limit. We also provide expansions close to threshold and make the Coulomb singularity explicit.
Our results can be downloaded from the website [94] where an easy-to-use routine is provided which provides numerical results for all six colour factors, both for the singlet and non-singlet contributions. We provide results for each individual colour factor which makes it straightforward to specify to QED.     F p, f, (3)

A Three-loop on-shell integrals to higher orders in
For the expansion around s = 0 we can analytically fix the boundary conditions, since in this case all integrals reduce to massive three-loop on-shell propagators. The most recent calculation of these integrals is given in Ref. [24] and extends to weight 7, formally enough for four-loop calculations. However, since we encounter spurious poles in our reductions of the full system, this is not enough to fix all boundary constants. Some integrals have to be calculated up to weight 9.
A subset of the integrals can be found in the appendix of Ref. [18]. Translating to the notation of Ref. [24] we can fix I 12 = G 45 , I 11 = G 46 , and I 9 = G 52 . But we still need the expansions of the integrals G 43 , G 51 , G ( ) 53 , G 61 , G ( ) 63 , G 64 , G 66 up to weight 9 (integrals with a ( ) actually only to weight 8). In the following we will describe the steps of the calculation and provide results for all master integrals which cannot be expressed in terms of Γ-functions up to weight 9.
In a first step we use the Mathematica package Summertime [97] to calculate the needed integrals up to the necessary order in with 5000 digits accuracy, for the integral G 8 we only obtain with 3000 digits accuracy. This is the input we use to apply the PSLQ algorithm [98] to fix the analytic form of the expansions. In order to apply this algorithm we need a basis of constants. In all the orders obtained by now, the basis of constants given by harmonic polylogarithms evaluated at argument 1 was sufficient. We use the basis given in Ref. [99]. This basis is convenient since the harmonic polylogarithms can be calculated to (in principle) arbitrary precision with ginac. For the usage in the PSLQ algorithm we have calculated them with 7500 digits accuracy. Since we are not dealing with integrals of uniform transcendentality one has to use all products of constants up to the desired weight in the ansatz. Therefore, the number of unknown constants from weight 1 to weight 9 is given by 2,4,7,12,20,33,54,88,143. After the successful reconstruction we change to the notation of SUMMER [100]. This has the advantage that most of the transcendental constants are known to Mathematica. The additional constants are given by: − 5376a 4 ζ 3 + 21504a 5 ζ 3 − 280π 2 l 2 ζ 3 + 2408 15   , c.f. Eq. (31). Note that the 1/ 3 pole of the colour factor C 2 A C F is zero.  , c.f. Eq. (31). Note that the 1/ 3 pole of the colour factor C 2 A C F is zero. C Analytic results for s → 0 In the following we present analytic expansion s → 0 of three-loop term for the non-singlet form factors F v 2 , F a 1 , F a 2 and F p . The results for F v 1 and F s can be found in Section 6.2.1. Our results read