Moments of heavy quark correlators with two masses: exact mass dependence to three loops

We compute moments of non-diagonal correlators with two massive quarks. Results are obtained where no restriction on the ratio of the masses is assumed. Both analytical and numerical methods are applied in order to evaluate the two-scale master integrals at three loops. We provide explicit results for the latter which are useful for other calculations. As a by-product we obtain results for the electroweak $\rho$ parameter up to three loops which can be applied to a fourth generation of quarks with arbitrary masses.


QCD corrections to current correlators
In recent years many multi-loop results to current correlators became available. In the massless approximation the state-of-the art are four-loop integrals which have been applied to τ and Z-boson decays in order to extract α s [1,2]. On the other hand, in the limit of small external momentum the so-called moments have been obtained also up to four-loop order [3][4][5][6][7][8][9][10][11] and both α s and the charm and bottom quark masses could be extracted by comparing to experimental data or lattice simulations [12][13][14][15][16]. The currents involved in these calculations only couple to one quark flavour.
As far as the non-diagonal correlators, where the external currents couple to two quark flavours with different masses, are concerned results for the moments up to three loops are known for the limit where one mass is zero [17][18][19][20][21][22]. At four-loop order the vector and axial-vector correlators have only been computed for vanishing external momentum (see Ref. [23,24]) in order to obtain four-loop corrections to the electroweak ρ parameter in the Standard Model [23][24][25].
It has been argued in Ref. [26] that it would be desirable to have at hand moments of the heavy-light currents also for general values of the two quark masses. The reason to not only consider physical mass values of the involved quarks is related to lattice simulations where often a variety of different masses are considered in order to be able to extrapolate to the mass values of physical interest. In Ref. [22] this task has been solved with the help of expansions around the equal-mass case and around the limit where one mass is much smaller than the other. In this paper we present exact results up to three loops and thus improve the findings of [22].
We adopt the notation from Ref. [22]. Let us nevertheless present the relevant formulas in order to have a self-contained presentation of the results.
In our calculations we are allowed to use anti-commuting γ 5 since for ψ 1 = ψ 2 only nonsinglet diagrams contribute.
We work in n f -flavour QCD with two massive and n f − 2 massless quarks. In our final results we adopt the MS renormalization scheme both for the parameters (strong coupling constant α s and the quark masses m 1 and m 2 ) and the overall renormalization of the correlator which in the following is indicated by a bar.
In the limit q 2 → 0 the quantityΠ δ (q 2 ) (and analogouslyΠ δ L (q 2 )) can be cast in the form where the dimensionless variables 1 have been introduced. Note that we assume 0 ≤ x ≤ 1. Results for m 2 > m 1 are easily obtained by interchanging m 1 and m 2 . We refer to the coefficientsC δ n (x) also as moments. Their perturbative expansion can be written as where the arguments x and µ (renormalization scale) are suppressed. Note that we have the symmetry relationsC s n (x) =C p n (−x) andC v n (x) =C a n (−x) (see, e.g., Ref. [27]) which we use to cross check our results. A further check is provided by computing the moments of the longitudinal contribution Π v,a L (q 2 ) and compare the result of the moments of Π s,p (q 2 ) with the help of the relationsC For completeness, let us mention that the QED-like normalization of the polarization function with Π δ (0) = 0 can be obtained fromΠ δ (q 2 ) with the help of Sample Feynman diagrams occurring at one-, two-and three-loop order are shown in Fig. 1.

Master integrals for non-diagonal correlators
For the evaluation of the Feynman diagrams a well-tested chain of programs has been used which works hand-in-hand in order to avoid error-prone manual interactions. Since the reduction to master integrals is nowadays pretty standard we restrict ourselves to a brief description up to this step. The Feynman diagrams are generated with QGRAF [28] and the various diagram topologies are identified and transformed to FORM [29,30] notation 1 In the following we do not write m  with the help of q2e and exp [31,32]. In a next step appropriate projectors are applied, the expansion in the external momentum is performed and traces are taken. Afterwards the scalar integrals are mapped to functions which are then passed to FIRE [33] where the reduction to master integrals is performed.
The master integrals which are needed for our calculation are shown in Figs. 2 and 3. At one-and two-loop order there are one and three master integrals, respectively, while at three-loop order the reduction leads to three single-scale (c.f. Fig. 2) and 16 two-scale master (c.f. Fig. 3) integrals. All one-and two-loop integrals and the single-scale threeloop integrals in Fig. 2 are available in the literature (see, e.g., Ref. [34,35]). For the three-loop integrals in Fig. 3 only partial results exist. In this work these results are checked by independent calculations and results for all integrals are presented. Since the problem at hand is symmetric under the exchange of m 1 and m 2 most master integrals come in two variants which are marked by the subscripts "a" and "b".

Differential equation method
Most of the master integrals can be calculated with the method of differential equations (see, e.g., Refs. [36][37][38][39][40][41]). In the following we exemplify the method for the two-loop master integral J (2) 1 of Fig. 3 which is given by where d = 4 − 2ǫ is the space-time dimension. The differential equations are derived by taking the derivative with respect to x and reducing the resulting integrals to master integrals. This leads to d dx J where K (1) (x) denotes the one-loop integral in Fig. 2 with mass xm 1 which is given by with the normalization To solve the differential equation it is necessary to know the initial condition at either x = 0 or x = 1, which corresponds to the single-scale integrals K (2) 1 and K (2) 2 shown in Fig. 2. In this simple example it is possible to solve the differential equation without expanding in ǫ with the result This result agrees with the one given in Ref. [34]. In general it is not possible to obtain a closed solution as in Eq. (12). In those cases one has to expand the original differential equation in ǫ and thus obtains simpler differential equations for the corresponding J (2) 1 Figure 3: Two-scale master integrals where thick, thin and dotted lines denote heavy, light and massless particles. A dot on a thick line indicates that the corresponding propagator is squared and a cross marks propagators which are raised to power minus one. coefficients. In our example one gets from Eq. (9) for the first three coefficients where the J 1,n (x) denotes the coefficient of the n th order in the expansion in ǫ. One observes that the structure of the homogeneous part remains the same whereas the inhomogeneous contribution becomes successively more complicated. While solving the differential equations it is advantageous to use Harmonic Polylogarithms (HPLs) [42], which form a special class of functions suitable for these kind of problems. It is particularly convenient to use the implementation in Mathematica from Ref. [43,44] which allows the application of the method of the variation of constants in order to solve the differential equation in Eq. (13). For our example 2 at hand the result looks as follows which can also be obtained by expanding Eq. (12) and agrees with Ref. [34].
The method of differential equations described above can immediately be applied to the the three-loop integrals J 3 , J 7a and J 7b . J 1 is given in Ref. [41]. In this reference also the result for an integral related to J by integrationby-parts identities is listed. We find agreement including terms of order ǫ.
To our knowledge the results for J J 7b can be obtained from J There are more checks to verify the obtained results. First of all we have checked that the results for the master integrals satisfy the original differential equations. Furthermore, since it is sufficient to use the value of the integral at x = 0 or x = 1 as initial condition for the differential equation the value at the other boundary can be used as cross check.
9b ) obey a system of two coupled differential equations which could not be solved analytically using the method described above. Providing initial conditions a numerical solution is possible, however, the achieved accuracy for the master integrals is not sufficient to compute higher order moments of the current correlators since there are large numerical cancellations between contributions from different master integrals. Thus we decided to apply the Mellin-Barnes method [46,47] which provides for a given value of x a high-precision numerical result with about 30 significant digits. Considering the integral J where the integration contour is chosen to separate the poles from the Gamma functions of the form Γ(z 1 + z 2 + . . .) from the ones originating from Γ(−z 1 + . . .) or Γ(−z 2 + . . .).
We used the package MB [48] to evaluate the Mellin-Barnes representation numerically for 0.1 ≤ x ≤ 1 in steps of 0.005 with high precision. 3 The results can be interpolated to obtain results for all values of x. An important cross check of the numerical approach constitutes the comparison of the exact analytical results for the poles which can be obtained from asymptotic expansions or the Mellin-Barnes representation. They are given by Note that for these integrals the order ǫ 1 terms are not needed for the final results.

J
10a and J (3) 10b The differential equations for J 10a and J 10b only contain inhomogeneous parts involving the integrals J 9b , respectively, and can therefore not be solved analytically. However, in this case the straightforward numerical solution of the differential equations would lead to results which are sufficiently precise for the moments considered in this paper.
A more precise result for J (3) 10a and J (3) 10b which is furthermore simpler to handle can be obtained with the help of expansions around x = 0 and x = 1. In contrast to Ref. [22] where expansions of the whole correlator have been considered and thus the expansion depth was limited to 8 and 9 terms, respectively, in this paper results up to order x 38 and (1 − x) 28 could be obtained for the scalar integrals J Also for these integrals the order ǫ 1 terms do not enter the final results. Thus numerical expressions are only needed for the finite contributions to the moments whereas the cancellation of the poles can be checked analytically.
All available information about the master integrals is provided in the above mentioned Mathematica file TwoMassTadpoles.m [45].

Moments
In this Section we present numerical results for the first four moments of the vector, axialvector, scalar and pseudo-scalar current correlators. We adapt the colour factors to SU (3) and set the number of massless quarks to three. Furthermore, we set the renormalization  For comparison we also show in Figs. 5 to 8 the results obtained in Ref. [22] as dashed lines. In Ref. [22] expansions around x = 0 and x = 1 have been computed which showed good convergence properties up to x ∼ < 0.1 and x ∼ > 0.5, respectively. Inbetween an interpolation has been performed.
In most cases good agreement is found. Small deviations for x ≈ 0.2 . . . 0.3 are observed forC (2),v 4 and, to a lesser extent, also forC

ρ parameter
In this section we discuss an immediate application of C δ −1 of the vector and axial-vector current, namely the QCD corrections to the electroweak ρ parameter which is given by with Σ W (0) and Σ Z (0) are the transverse parts of the W and Z boson propagators defined through where Π µν W/Z are the corresponding polarization functions. The QCD corrections to δρ in the Standard Model (SM), where all quark masses except the top quark mass are set to zero, are known up to four loops [23][24][25][49][50][51][52]. In the following we consider a generic fourth generation of quarks which couples to the W and Z boson in the same way as the top and bottom quarks of the SM. For simplicity we denote the new doublet by (t ′ , b ′ ).
In the case of Σ Z (0) there is no contribution form the vector part. Furthermore one has to consider additional contributions at three-loop order which are only present for currents where ψ 1 = ψ 2 in Eq. (2). Thus we cast the Z-boson self energy in the form where the quantities δC a,(2) −1,db and δC a,(2) −1,sing receive contributions from the double bubble and singlet Feynman diagrams involving two mass scales (see Fig. 9). The double bubble diagrams involving only one massive quark are contained inC a −1,diag together with all other one-scale contributions. Note that the singlet-type contributions involving massless up, down, strange and charm quarks add up to zero.
In Appendix A explicit analytical results are presented forC a −1,diag , δC a,(2) −1,db and δC a,(2) −1,sing . Let us mention that the treatment of γ 5 for the singlet diagrams can be found in Ref. [53].
Before presenting results we transform the MS quark masses m t ′ and m b ′ to the on-shell scheme [54,55] since this is common practice in the context of analyzing electroweak precision data. For M b ′ = 0 we reproduce the SM result [51,52]. The general result can be cast in the form (a) Similarly, the expansion around M b ′ ≈ M t ′ leads to In the file rho.m [45] analytical results up to order X 30 and (1 − X) 20 are provided.
In the above equations n l labels the number of massless quark flavours. In order to reproduce the SM result with massless bottom quark one has n l = 4 and x = 0. The contribution of a fourth generation is obtained for n l = 6 which assumes a massless top quark at three-loop order.
In the SM the higher order corrections to the ρ parameter are quite important. E.g., the three-loop corrections correspond to a shift in the top quark mass of about 2 GeV which is larger than the current experimental uncertainty. In the limit of equal quark masses δρ is proportional to the mass difference and thus the numerical impact of the higher order corrections is reduced as can be seen from Eq. (32).

Summary
In this paper moments of correlators formed by vector, axial-vector, scalar and pseudoscalar currents are considered to three-loop order. The currents couple to quarks with different masses m 1 and m 2 which leads to two-scale vacuum integrals after expanding in the external momentum of the correlators. The moments can be expressed as a linear combination of 12 non-trivial three-loop master integrals which are discussed in details in Section 2. In particular, a Mathematica package TwoMassTadpoles.m is provided [45] which contains analytical results for all but six master integrals. For the latter highprecision numerical results are available. In a further Mathematica file coefhl2.m we also provide results for the first four moments for each correlator.
As a by-product we compute three-loop corrections to the ρ parameter for a generic fourth generation of quarks with masses m t ′ and m b ′ for the up-and down-type quark. In addition to the calculation for the moments of the non-diagonal correlators one has to compute double-bubble and singlet-type diagrams where two different masses can be present in the fermion triangles. For the latter analytical results are provided; the complete results for the ρ parameter up to three loops is available in rho.m [45].
The three-loop resultC a,(2) −1,diag depends on the number of massless flavours n l . Let us finally mention that the Mathematica package rho.m contains replacement rules which express the HPLs occurring in Eqs. (34) and (35) in terms of logarithms and dilogarithms.