Validity of Pad´e approximations in vacuum polarization at three- and four-loop order

The heavy-quark contribution to the polarization function Π at higher perturbative orders is presently only known approximately. We scrutinize 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 nonsinglet heavy-quark vacuum polarization in order to test the prediction for this moment based on the Pad´e approximation.


I. INTRODUCTION
Vacuum polarization 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 multiloop 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 normalization factor, the cross section is equal to the imaginary part of the quark contribution to the vacuum polarization. More precisely, the heavy-quark polarization function Π and the cross section are related via RðsÞ sðs − s 0 Þ ; 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 polarization function, there is a contribution from flavor-singlet diagrams with massless cuts [3]. These cuts do not correspond to the production of heavy quarks. In the following, we will therefore restrict ourselves to the discussion of the nonsinglet polarization 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 polarization function is known at fourloop order [4]. The closely connected Adler function DðsÞ ¼ −12π 2 s d ds ΠðsÞ is even known at five-loop order for massless quarks [5,6]. The dimensionless polarization function can only depend on the energy through logarithms, which in turn give rise to the complete imaginary part of the polarization function. Thus, as per the optical theorem (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 [7]. 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,8].
In the kinematic region where the quark mass is nonnegligible much less is known about the vacuum polarization corrections than in the limit of massless quarks. The first major step towards obtaining the three-loop corrections was taken about 20 years ago [9], 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 highenergy expansions have become available [10][11][12][13], allowing a systematic improvement of the approximation (see e.g. Ref. [14]). An alternative approximation procedure based on Mellin-Barnes transforms was explored in Ref. [15]. Independently, the cross section corresponding to the imaginary part of the vacuum polarization was computed numerically in Refs. [16,17]. Corrections involving both massive and massless quarks were obtained already much earlier in Ref. [18].
At four-loop order, the same approaches were again used for an approximate reconstruction of the heavy-quark corrections to the vacuum polarization [15,[19][20][21]. These approximations were in turn expanded again in the lowenergy 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 [22][23][24][25][26][27][28][29] the exactly computed first three physical moments [30][31][32][33] were considered together with an estimate of the fourth moment.
To summarize, current knowledge of quark-mass corrections to the vacuum polarization 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 Refs. [9,19,20] 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 analyze 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 polarization 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 Ref. [15]. 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 [19][20][21] to the four-loop polarization function.

II. CONVENTIONS
The quark contribution to the vacuum polarization is given by the correlator of two vector currents, viz.
where the vector current is j μ ¼ψγ μ ψ. The polarization function Π is conventionally renormalized 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 renormalization scale μ ¼ m Q , where m Q is the quark mass renormalized in the on-shell scheme [34][35][36][37][38][39][40][41].
In the following we are interested in the three-loop coefficient Π ð2Þ and the four-loop coefficient Π ð3Þ . The kinematic dependence of the polarization 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 polarization 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 jzj < 1.

III. CALCULATIONAL SETUP
We generate the diagrams contributing to the polarization function with QGRAF [45], 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 [46]. Color factors are computed with the COLOR [47] 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 [48] according to Laporta's algorithm [49] as implemented in CRUSHER [50,51].
At four loops, the expansion around z ¼ 0 results in vacuum integrals, and the resulting master integrals are known analytically [54][55][56][57][58][59][60][61][62][63]. At three loops, we derive differential equations [64][65][66] for the master integrals expanded in the dimensional regularization parameter ϵ. We solve the differential equations using the Runge-Kutta-Dormand-Prince [67] method as implemented in the Odeint C++ library [68]. As a boundary condition we choose values of the integrals at z 0 ≈ 0, which we obtain from the low-energy expansion performed in Ref. [12]. 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 Þ þ isgnðImðzÞÞReðzÞ and ReðzÞ þ isgnðImðzÞÞReðzÞ to z.

IV. THREE-LOOP QUARK CONTRIBUTION TO THE POLARIZATION FUNCTION
In the following, we present our new result for the threeloop polarization function and compare to approximations based on previously known expansion coefficients.

A. Comparison to Padé-based approximation
We construct Padé-based approximants according to the procedure described in Ref. [69]. We briefly summarize the main aspects. First, we use subtraction functions listed in Ref. [69] to split Π ð2Þ into two parts, where all known logarithms and poles in the threshold and high-energy expansions (5) and (6) have been absorbed into Π ð2Þ log . We then make a Padé ansatz of the form where the variable ω is defined by the relation 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 ð2Þ n;0 , so N ¼ 61. Note that the threshold expansion (5) is only used in the construction of Π ð2Þ log . In particular, terms that are analytic in ffiffiffiffiffiffiffiffiffiffi 1 − z p are not considered for the approximation.
Further approximants are then constructed with the recurrence relations [70] h whereη j is the numerator of the approximant in the form of Eq. (8) andθ j is 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 flavors, we decompose and construct separate approximations for the n l -independent coefficients Π . Since the two latter approximants are numerically essentially indistinguishable, we somewhat arbitrarily select ½32=28.
In Fig. 1 we compare the Padé-based approximants to the exact result, which we compute as described in Sec. III. 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 is typically of the order of the numeric precision requested when solving the differential equation. We conclude that VALIDITY OF PADé APPROXIMATIONS IN VACUUM … PHYS. REV. D 97, 056016 (2018) 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 ð2Þ n , D ð2Þ m;0 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.

B. Comparison to the approximation based on the Mellin-Barnes transform
In Fig. 2 we compare the exact result to the approximation of Ref. [15], which is based on the Mellin-Barnes transform. In Ref. [15], 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 highenergy expansion are employed in the construction of the FIG. 1. Comparison for Π ð2Þ between the 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 flavors, whereas on the right the corrections including a virtual massless quark loop are considered. Note that the differential equations contain spurious singularities for z ∈ f0; 0.25; 0.5g.

FIG. 2.
Comparison for Π ð2Þ with n l ¼ 3 massless quark flavors between the approximation of Ref. [15] (dotted) and results obtained by numerically solving differential equations with a requested absolute error of 10 −10 (solid lines). 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 Ref. [15]. Similarly to Sec. IVA, 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 P ∞ n¼1 log n n ðAEωÞ n , which are difficult to evaluate close to the branch cut jωj ¼ 1.
As for the Padé-based approximation the agreement in the low-energy region z < 1 is 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 [13] would improve the precision in this region further.

V. LOW-ENERGY EXPANSION AT FOUR LOOPS
In the following we compare our new analytic result for C ð3Þ 4 to various estimates. A similar comparison at threeloop order using restricted input in the construction of a Padé-based approximation was already performed in Ref. [19], where good agreement between the approximate and exact results for C ð2Þ 4 was found. The low-energy expansion coefficients C ð3Þ n can be decomposed according to their color structure [72]: ll;n þ C F T 2 F n l n h C lA;n þ C F C sing;n : 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 normalization 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 color factors are The number of quark flavors with mass m Q is denoted by n h . The remaining factors in Eq. (14) are the dimension of the fundamental representation D F and d FF 33 ¼ ½ 1 2 TrðT a T b T c þ T a T c T b Þ 2 . However, this color structure only appears in the flavor-singlet contribution. As already mentioned in Sec. I, we will therefore not consider the coefficient C sing;n .
Since the bosonic contribution for C 3 has only been presented for a SU(3) gauge group in previous works [32,33], we provide the general color decomposition in the Appendix. Our new result for C where ζ n ¼ P ∞ k¼1 1 k n denote values of the Riemann ζ function, the auxiliary constants c 4 , c 5 are given by c 4 ¼ 24Li 4 1 2 þ log 4 ð2Þ − π 2 log 2 ð2Þ; ð25Þ and Li n ð 1 2 Þ ¼ P ∞ k¼1 1 2 k k n are values of polylogarithm functions. The corresponding results for the coefficients in the MS scheme are given in Appendix A 2. The expressions in both schemes are also available in computer-readable form as Supplemental Material to this article [71].
In Table I we compare the numerical values for QCD with n l ¼ 3, 4, 5 to the estimates obtained in Refs. [19][20][21]. We find excellent agreement, especially for the predictions from Ref. [20]. In fact, the true approximation error of Ref. [20] seems to be almost an order of magnitude less than estimated.

VI. CONCLUSION
We have tested the quality of three-and four-loop approximations for the quark contribution to the vacuum polarization. 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 Ref. [15]. 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 found excellent agreement with the Padé-based prediction [20], well within the error estimate. Within their errors, the less precise estimates from Refs. [19,21] also agree well with the exact result. between predictions from Refs. [19][20][21] and the exact analytic result for different numbers of light quark flavors.C ð3Þ 4 is the coefficient in the MS scheme. We have refrained from converting the results from Refs. [19,21] to different schemes.

Coefficients in the MS scheme
Renormalizing 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 [73][74][75][76][77][78] at the scale μ ¼m Q . The analytic results arē ll;n þ C F T 2 F n l n hC