Life of \Pi

The heavy-quark contribution to the polarisation function $\Pi$ at higher perturbative orders is presently only known approximately. We scrutinise the accuracy of state-of-the-art approximations at three- and four-loop order. At three loops, we present for the first time a result with arbitrary numerical precision for general kinematics and compare to the best Pad\'e estimate. At four loops, we calculate the fourth (inverse) moment of the non-singlet heavy-quark vacuum polarisation in order to test the prediction for this moment based on Pad\'e approximation.


Introduction
Vacuum polarisation is one of the earliest and phenomenologically most important predictions of quantum electrodynamics (QED). Consequently, the computation of the two-loop perturbative correction to this effect constitutes one of the very first multi-loop calculations performed within QED [1].
Quantum corrections mediated through virtual quarks are of special interest. They are closely connected to the total inclusive hadron production cross section at lepton colliders through a dispersion relation [2]. Conversely, it follows from the optical theorem that, up to a simple normalisation factor, the cross section is equal to the imaginary part of the quark contribution to the vacuum polarisation. More precisely, the heavy-quark polarisation function Π and the cross section are related via where the R-ratio for a heavy quark Q is defined as R(s) = σ(e + e − → QQX)/σ 0 with σ 0 = 4πα 2 3s . Starting at four loops in the perturbative expansion of the polarisation function, there is a contribution from flavour-singlet diagrams with massless cuts. 1 These cuts do not correspond to the production of heavy quarks. In the following, we will therefore restrict ourselves to the discussion of the non-singlet polarisation function.
In the limit where the center-of-mass energy is far above both the scale of nonperturbative dynamics and the masses of the quarks, the polarisation function is known at four-loop order [3]. The closely connected Adler function D(s) = −12π 2 s d ds Π(s) is even known at five-loop order for massless quarks [4,5]. The dimensionless polarisation function can only depend on the energy through logarithms, which in turn give rise to the complete imaginary part of the polarisation function. Thus, as per the optical theorem Eq. (1.1), the knowledge of the five-loop Adler function allows a N 4 LO prediction of the total quark production cross section.
However, in the production of heavy quarks the approximation of small quark masses is not always justified. In fact, sufficiently close to the production threshold the full quark mass dependence has to be taken into account. A prominent scenario is the production of top-antitop pairs at the projected first stage of CLIC at a center-of-mass energy of 380 GeV [6]. For the determination of the charm-and the bottom-quark mass it is even the opposite limit of large quark masses (or small center-of-mass energies) that is most relevant. The coefficients in such a low-energy expansion can be identified with (inverse) moments of the heavy-quark production cross section via the aforementioned dispersion relation. These moments in turn are the main ingredient in sum-rule determinations of the quark masses [2,7].
In the kinematic region where the quark mass is non-negligible much less is known about the vacuum polarisation corrections than in the limit of massless quarks. The first major step towards obtaining the three-loop corrections was taken about 20 years ago [8], when expansions in the low-energy, threshold, and high-energy kinematic regions were exploited to construct a Padé-based approximation. Since then, many more terms in the low-energy and high-energy expansions have become available [9][10][11][12], allowing a systematic improvement of the approximation (see e.g. [13]). An alternative approximation procedure based on Mellin-Barnes transforms was explored in [14]. Independently, the cross section corresponding to the imaginary part of the vacuum polarisation was computed numerically in [15,16]. Corrections involving both massive and massless quarks were obtained already much earlier in [17].
At four-loop order, the same approaches were again used for an approximate reconstruction of the heavy-quark corrections to the vacuum polarisation [14,[18][19][20]. These approximations were in turn expanded again in the low-energy limit in order to obtain estimates for higher moments used in sum-rule analyses. In the most precise determinations of the charm-and bottom-quark masses from relativistic sum rules to date [21][22][23][24][25][26][27][28] the exactly computed first three physical moments [29][30][31][32] were considered together with an estimate of the fourth moment.
To summarise, current knowledge of quark-mass corrections to the vacuum polarisation at three-and four-loop order is based to a large degree on approximations. If and in which sense approximations based on the scheme considered in [8,18,19] converge to the true results as more information is added is an open question, which we do not intend to address in this work. Our goal is rather to analyse to which extent approximations based on current knowledge and their heuristic error estimates can be relied on. We aim to ameliorate the dependence on approximations by providing new exact results at three and four loops. At three-loop order we numerically calculate the vacuum polarisation for general kinematics and compare to a new Padé-based approximation constructed from many coefficients in the low-and high-energy expansions as well as to the approximation obtained in [14]. At four loops, we present an analytic result for the fourth term in the low-energy expansion and compare to the various estimates based on the approximations [18][19][20] to the four-loop polarisation function.

Conventions
The quark contribution to the vacuum polarisation is given by the correlator of two vector currents, viz.
where the vector current is j µ =ψγ µ ψ. The polarisation function Π is conventionally renormalised in the on-shell scheme, so that Π(0) ≡ 0. Its perturbative expansion in the strong coupling constant α s for a quark with charge e Q can be written as We set the renormalisation scale µ = m Q , where m Q is the quark mass renormalised in the on-shell scheme [33][34][35][36][37][38][39][40].
In the following we are interested in the three-loop coefficient Π (2) and the four-loop coefficient Π (3) . The kinematic dependence of the polarisation function is described by a single ratio of energy and mass, which we define as z = q 2 /(4m 2 Q ). We consider the general case of complex z, which is needed for example when describing unstable quarks, like the top quark. The perturbative coefficients of the polarisation function are analytic functions of z, apart from a branch cut along the positive real axis. Since we neglect contributions from diagrams with massless cuts, the branch cut starts at the open quark production threshold z = 1. The low-energy expansions of the perturbative coefficients therefore converge for |z| < 1. In close analogy, we write the expansions in the threshold region z → 1 and the high-energy region z → −∞ as The three-loop coefficients C (2) n , D n,m are known up to n = 30 [10][11][12]. At four loops, the coefficients C can be extracted from calculations in a non-relativistic effective theory [41][42][43]; explicit expressions obtained from NNLO results are given in [18,19].

Calculational setup
We generate the diagrams contributing to the polarisation function with QGRAF [44], obtaining 36 diagrams at three loops and 700 diagrams at four loops. For inserting the Feynman rules, evaluating traces, and performing general symbolic manipulations we use FORM [45]. Colour factors are computed with the color [46] package. At four loops, we also perform an expansion around z = 0 up to order z 4 . The resulting scalar integrals are reduced to master integrals by exploiting integration-by-parts identities [47] according to Laporta's algorithm [48] as implemented in Crusher 2 [51].
At four loops, the expansion around z = 0 results in vacuum integrals, and the resulting master integrals are known analytically [52][53][54][55][56][57][58][59][60][61]. At three loops, we derive differential equations [62][63][64] for the master integrals expanded in the dimensional regularisation parameter . We solve the differential equations using the Runge-Kutta-Dormand-Prince [65] method as implemented in the odeint C++ library [66]. As boundary condition we choose values of the integrals at z 0 ≈ 0, which we obtain from the low-energy expansion performed in [11]. Note that we avoid z 0 = 0, since the differential equations exhibit a singularity at this point. For general complex z, we integrate the differential equations along a straight line from z 0 to z. However, there are further singularities along the positive real axis, even below the physical branch cut starting at z = 1. When z is close to the real axis, we therefore perform a contour deformation into the complex plane. In principle, any path that bypasses the singularities is sufficient. In practice, we choose a piecewise linear path from z 0 over Re(z 0 ) + i sgn(Im(z)) Re(z) and Re(z) + i sgn(Im(z)) Re(z) to z.

Three-loop quark contribution to the polarisation function
In the following, we present our new result for the three-loop polarisation function and compare to approximations based on previously known expansion coefficients.

Comparison to Padé-based approximation
We construct Padé-based approximants according to the procedure described in [67]. We briefly summarise the main aspects. First, we use subtraction functions listed in [67] to split Π (2) into two parts, where all known logarithms and poles in the threshold and high-energy expansions Eqs. (2.4), (2.5) have been absorbed into Π log . We then make a Padé ansatz of the form where the variable ω is defined by the relation (4. 3) The approximants [N/0], [N − 1/0] are fixed by requiring that the coefficients match the terms in the Maclaurin series of The degree N corresponds to the number of known coefficients C (2) n , D n,0 , so N = 61. Note that the threshold expansion Eq. (2.4) is only used in the construction of Π (2) log . In particular, terms that are analytic in √ 1 − z are not considered for the approximation. Further approximants are then constructed with the recurrence relations [68] whereη j is the numerator of the approximant in the form of Eq. (4.2) andθ j its denominator. We discard all approximants with poles inside the unit circle, which translate to unphysical poles in the variable z. Instead of constructing new approximants for various fixed numbers n l of massless quark flavours, we decompose and construct separate approximations for the n l -independent coefficients Π . After discarding unphysical approximants as described above we obtain 80 approximants for each Π (2) n 0 l and Π (2) n 1 l . The expressions for the approximants are quite lengthy and provided as ancillary files. Diagonal approximants with n = m are generally expected to perform best, so we select the Padé approximants [n/m] that minimise the distance |n − m| for the following comparison. This corresponds to [30/30] for Π (2) n 0 l and either of [32/28] or [28/32] for Π (2) n 1 l . Since the two latter approximants are numerically essentially indistinguishable, we somewhat arbitrarily select [32/28].
In Figure 1 we compare the Padé-based approximants to the exact result, which we compute as described in Section 3. For the sake of a clear presentation, we restrict ourselves to real values of z, choosing the physical branch on the upper complex half-plane for z > 1. This is implemented in the numerical evaluation by adding a small imaginary part, i.e. by shifting the argument z → (1 + 10 −10 i)z. We solve the differential equations for 198 values of z, which requires about 14 seconds on a single core of an Intel Core i5-4200M processor. It should be noted that the time required for the calculation of a single point increases greatly in the vicinity of singularities. We find excellent agreement over the whole kinematic range, including the region around the Coulomb singularity at z = 1. In fact, the difference Figure 1. Comparison for Π (2) between Padé-based approximation (dotted) and results obtained by numerically solving differential equations with a requested absolute error of 10 −10 (solid lines).
The panel on the left shows the corrections without any light quark flavours, whereas on the right the corrections including a virtual massless quark loop are considered. Note that the differential equations contain spurious singularities for z ∈ {0, 0.25, 0.5}. is typically of the order of the numeric precision requested when solving the differential equation. We conclude that for all practical purposes the approximation is indistinguishable from the true result.
With this degree of accuracy, it is also possible to omit a number of expansion terms in the construction of the approximation while still retaining agreement with the exact result at the level of 10 −10 . For instance, we find that limiting ourselves to coefficients C with n, m < 22 does not lead to a visible increase in the deviation. When omitting further coefficients the accuracy degrades notably in the region above threshold, e.g. to the level of 10 −9 for a [20/20] approximant constructed from coefficients with n < 20, m < 19 and 10 −8 for a [16/15] approximant from expansion terms with n, m < 15.

Comparison to approximation based on Mellin-Barnes transform
In Figure 2 we compare the exact result to the approximation of [14], which is based on the Mellin-Barnes transform. In [14], a flexible number of N * coefficients in the low-energy expansion, all known coefficients in the threshold expansion, and terms up to order z −2 in the high-energy expansion are employed in the construction of the approximation. For the comparison we take into account all N * = 30 low-energy coefficients, but make no attempt at improving the approximation over what was done in [14]. Similarly to section 4.1, we focus on values of z that are close to the real axis. However, we choose a somewhat larger imaginary part by shifting z → (1 + 0.01i)z in both the approximation and the exact result. The reason for this is that the expression for the approximant contains sums of the form ∞ n=1 log n n (±ω) n , which are difficult to evaluate close to the branch cut |ω| = 1. As for the Padé-based approximation the agreement in the low-energy region z < 1 is Im Π (2) Figure 2. Comparison for Π (2) with n l = 3 massless quark flavours between the approximation of [14] (dotted) and results obtained by numerically solving differential equations with a requested absolute error of 10 −10 (solid lines).
remarkably good. Above the threshold, the difference is of the order of 10 −4 , bigger than for the Padé-based approximation. We expect that the inclusion of the complete known high-energy expansion up to z −30 [12] would improve the precision in this region further.

Low-energy expansion at four loops
In the following we compare our new analytic result for C to various estimates. A similar comparison at three-loop order using restricted input in the construction of a Padé-based approximation was already performed in [18], where good agreement between the approximate and exact results for C (2) 4 was found. The low-energy expansion coefficients C n can be decomposed according to their colour structure: 3 hA,n + C F C hF,n ) sing,n . (5.1) As usual, C F and C A denote the eigenvalues of the quadratic Casimir operators in the fundamental and the adjoint representation, respectively. T F is the trace normalisation defined by Tr(T a T b ) = T F δ ab , where T b , T b are generators of the fundamental representation. For QCD, the values of these colour factors are C F = 4/3, C A = 3, T F = 1/2. The number 3 Our notation differs slightly from the one employed in e.g. [32], since the latter does not generalise well to the purely bosonic corrections.
of quark flavours with mass m Q is denoted by n h . The remaining factors in Eq. (5.1) are the dimension of the fundamental representation D F and d F F However, this colour structure only appears in the flavour-singlet contribution. As already mentioned in Section 1, we will therefore not consider the coefficient C (3) sing,n . Since the bosonic contribution for C 3 has only been presented for a SU(3) gauge group in previous works [31,32], we provide the general colour decomposition in Appendix A. Our new result for C Ref. [20] Ref. [ between predictions from [18][19][20] and the exact analytic result for different numbers of light quark flavours.C (3) 4 is the coefficient in the MS scheme. We have refrained from converting the results from [18] and [20] to different schemes. In Table 1 we compare the numerical values for QCD with n l = 3, 4, 5 to the estimates obtained in [18][19][20]. We find excellent agreement, especially for the predictions from [19]. In fact, the true approximation error of [19] seems to be almost an order of magnitude less than estimated.

Conclusion
We have tested the quality of three-and four-loop approximations for the quark contribution to the vacuum polarisation. To this end, we have calculated the three-loop contribution numerically, finding almost perfect agreement with a newly constructed Padé-based approximation and very good agreement with an approximation from [14]. At four loops, we have computed analytically the fourth term in the low-energy expansion, which is also relevant for relativistic sum-rule determinations of the charm-and bottom-quark masses. We find excellent agreement with the Padé-based prediction [19], well within the error estimate. Within their errors, the less precise estimates from [18] and [20] also agree well with the exact result.

A.1 Coefficients in the on-shell scheme
In the following, we show the four-loop coefficients C

A.2 Coefficients in the MS scheme
Renormalising the heavy-quark mass in the MS instead of the on-shell scheme we obtain the low-energy expansion wherem Q now denotes the MS mass [69][70][71][72][73][74] at the scale µ =m Q . The analytic results arē